From f8754c0f50d309c5f5b020ed5dc59333180d0e11 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 5 Nov 2019 16:11:18 +0100 Subject: [PATCH] LPSolvers: Allowing premature termination by specifying a mip gap. Fixes for incremental solving with Z3Lpsolver. --- src/storm/solver/GlpkLpSolver.cpp | 68 +++++++++++++++++++++++++---- src/storm/solver/GlpkLpSolver.h | 15 +++++++ src/storm/solver/GurobiLpSolver.h | 16 ++----- src/storm/solver/LpSolver.h | 14 ++++++ src/storm/solver/Z3LpSolver.cpp | 71 ++++++++++++++++++++++++++----- src/storm/solver/Z3LpSolver.h | 16 ++++--- 6 files changed, 163 insertions(+), 37 deletions(-) diff --git a/src/storm/solver/GlpkLpSolver.cpp b/src/storm/solver/GlpkLpSolver.cpp index 62518fba6..76007fc91 100644 --- a/src/storm/solver/GlpkLpSolver.cpp +++ b/src/storm/solver/GlpkLpSolver.cpp @@ -20,7 +20,6 @@ #include "storm/settings/modules/GlpkSettings.h" - namespace storm { namespace solver { @@ -37,6 +36,11 @@ namespace storm { // Set whether the glpk output shall be printed to the command line. glp_term_out(storm::settings::getModule().isDebugSet() || storm::settings::getModule().isOutputSet() ? GLP_ON : GLP_OFF); + // Set the maximal allowed MILP gap to its default value + glp_iocp* defaultParameters = new glp_iocp(); + glp_init_iocp(defaultParameters); + this->maxMILPGap = defaultParameters->mip_gap; + this->maxMILPGapRelative = true; } template @@ -213,6 +217,24 @@ namespace storm { this->currentModelHasBeenOptimized = false; } + // Method used within the MIP solver to terminate early + void callback(glp_tree* t, void* info) { + auto& mipgap = *static_cast*>(info); + double actualRelativeGap = glp_ios_mip_gap(t); + double factor = storm::utility::one(); + if (!mipgap.second) { + // Compute absolute gap + factor = storm::utility::abs(glp_mip_obj_val(glp_ios_get_prob(t))) + DBL_EPSILON; + assert(factor >= 0.0); + } + if (actualRelativeGap * factor <= mipgap.first) { + // Terminate early + mipgap.first = actualRelativeGap; + mipgap.second = true; // The gap is relative. + glp_ios_terminate(t); + } + } + template void GlpkLpSolver::optimize() const { // First, reset the flags. @@ -228,9 +250,22 @@ namespace storm { glp_init_iocp(parameters); parameters->presolve = GLP_ON; parameters->tol_int = storm::settings::getModule().getIntegerTolerance(); + + // Check whether we allow sub-optimal solutions via a non-zero MIP gap. + // parameters->mip_gap = this->maxMILPGap; (only works for relative values. Also, we need to obtain the actual gap anyway. + std::pair mipgap(this->maxMILPGap, this->maxMILPGapRelative); + if (!storm::utility::isZero(this->maxMILPGap)) { + parameters->cb_func = &callback; + parameters->cb_info = &mipgap; + } + + // Invoke mip solving error = glp_intopt(this->lp, parameters); delete parameters; + // mipgap.first has been set to the achieved mipgap (either within the callback function or because it has been set to this->maxMILPGap) + this->actualRelativeMILPGap = mipgap.first; + // In case the error is caused by an infeasible problem, we do not want to view this as an error and // reset the error code. if (error == GLP_ENOPFS) { @@ -239,6 +274,9 @@ namespace storm { } else if (error == GLP_ENODFS) { this->isUnboundedFlag = true; error = 0; + } else if (error == GLP_ESTOP) { + // Early termination due to achieved MIP Gap. That's fine. + error = 0; } else if (error == GLP_EBOUND) { throw storm::exceptions::InvalidStateException() << "The bounds of some variables are illegal. Note that glpk only accepts integer bounds for integer variables."; } @@ -282,13 +320,7 @@ namespace storm { return false; } - int status = 0; - if (this->modelContainsIntegerVariables) { - status = glp_mip_status(this->lp); - } else { - status = glp_get_status(this->lp); - } - return status == GLP_OPT; + return !isInfeasible() && !isUnbounded(); } template @@ -438,7 +470,25 @@ namespace storm { } } } - + + template + void GlpkLpSolver::setMaximalMILPGap(ValueType const& gap, bool relative) { + if (relative) { + this->maxMILPGap = storm::utility::convertNumber(gap); + } else { + this->maxMILPGap = storm::utility::convertNumber(gap); + } + } + + template + ValueType GlpkLpSolver::getMILPGap(bool relative) const { + STORM_LOG_ASSERT(this->isOptimal(), "Asked for the MILP gap although there is no (bounded) solution."); + if (relative) { + return this->actualRelativeMILPGap; + } else { + return storm::utility::abs(this->actualRelativeMILPGap * getObjectiveValue()); + } + } template class GlpkLpSolver; template class GlpkLpSolver; diff --git a/src/storm/solver/GlpkLpSolver.h b/src/storm/solver/GlpkLpSolver.h index be530034b..633a550ec 100644 --- a/src/storm/solver/GlpkLpSolver.h +++ b/src/storm/solver/GlpkLpSolver.h @@ -95,6 +95,9 @@ namespace storm { virtual void push() override; virtual void pop() override; + + virtual void setMaximalMILPGap(ValueType const& gap, bool relative) override; + virtual ValueType getMILPGap(bool relative) const override; private: /*! @@ -122,6 +125,10 @@ namespace storm { mutable bool isInfeasibleFlag; mutable bool isUnboundedFlag; + mutable double maxMILPGap; + mutable bool maxMILPGapRelative; + mutable double actualRelativeMILPGap; + struct IncrementalLevel { std::vector variables; int firstConstraintIndex; @@ -240,6 +247,14 @@ namespace storm { virtual void pop() override { throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without support for glpk. Yet, a method was called that requires this support. Please choose a version of support with glpk support."; } + + virtual ValueType getMILPGap(bool relative) const override { + throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without support for glpk. Yet, a method was called that requires this support. Please choose a version of support with glpk support."; + } + + virtual ValueType getMILPGap(bool relative) const override { + throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without support for glpk. Yet, a method was called that requires this support. Please choose a version of support with glpk support."; + } }; #endif diff --git a/src/storm/solver/GurobiLpSolver.h b/src/storm/solver/GurobiLpSolver.h index 5947ad600..fe6cac2f8 100644 --- a/src/storm/solver/GurobiLpSolver.h +++ b/src/storm/solver/GurobiLpSolver.h @@ -107,6 +107,9 @@ namespace storm { virtual void push() override; virtual void pop() override; + virtual void setMaximalMILPGap(ValueType const& gap, bool relative) override; + virtual ValueType getMILPGap(bool relative) const override; + // Methods to retrieve values of sub-optimal solutions found along the way. void setMaximalSolutionCount(uint64_t value); // How many solutions will be stored (at max) uint64_t getSolutionCount() const; // How many solutions have been found @@ -115,18 +118,7 @@ namespace storm { bool getBinaryValue(storm::expressions::Variable const& name, uint64_t const& solutionIndex) const; ValueType getObjectiveValue(uint64_t const& solutionIndex) const; - /*! - * Specifies the maximum difference between lower- and upper objective bounds that triggers termination. - * That means a solution is considered optimal if - * upperBound - lowerBound < gap * (relative ? |upperBound| : 1). - * Only relevant for programs with integer/boolean variables. - */ - void setMaximalMILPGap(ValueType const& gap, bool relative); - - /*! - * Returns the obtained gap after a call to optimize() - */ - ValueType getMILPGap(bool relative) const; + private: /*! diff --git a/src/storm/solver/LpSolver.h b/src/storm/solver/LpSolver.h index f397352fe..e92691921 100644 --- a/src/storm/solver/LpSolver.h +++ b/src/storm/solver/LpSolver.h @@ -294,6 +294,20 @@ namespace storm { * that were added after the last call to push(). */ virtual void pop() = 0; + + /*! + * Specifies the maximum difference between lower- and upper objective bounds that triggers termination. + * That means a solution is considered optimal if + * upperBound - lowerBound < gap * (relative ? |upperBound| : 1). + * Only relevant for programs with integer/boolean variables. + */ + virtual void setMaximalMILPGap(ValueType const& gap, bool relative) = 0; + + /*! + * Returns the obtained gap after a call to optimize() + */ + virtual ValueType getMILPGap(bool relative) const = 0; + protected: // The manager responsible for the variables. diff --git a/src/storm/solver/Z3LpSolver.cpp b/src/storm/solver/Z3LpSolver.cpp index a1a5fd967..a769dff47 100644 --- a/src/storm/solver/Z3LpSolver.cpp +++ b/src/storm/solver/Z3LpSolver.cpp @@ -9,6 +9,7 @@ #include "storm/utility/macros.h" #include "storm/utility/constants.h" +#include "storm/utility/file.h" #include "storm/storage/expressions/Expression.h" #include "storm/storage/expressions/ExpressionManager.h" @@ -31,7 +32,6 @@ namespace storm { context = std::unique_ptr(new z3::context(config)); solver = std::unique_ptr(new z3::optimize(*context)); expressionAdapter = std::unique_ptr(new storm::adapters::Z3ExpressionAdapter(*this->manager, *context)); - optimizationFunction = this->getManager().rational(storm::utility::zero()); } template @@ -71,7 +71,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getRationalType()); } solver->add(expressionAdapter->translateExpression((newVariable.getExpression() >= this->manager->rational(lowerBound)) && (newVariable.getExpression() <= this->manager->rational(upperBound)))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -84,7 +86,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getRationalType()); } solver->add(expressionAdapter->translateExpression(newVariable.getExpression() >= this->manager->rational(lowerBound))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -97,7 +101,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getRationalType()); } solver->add(expressionAdapter->translateExpression(newVariable.getExpression() <= this->manager->rational(upperBound))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -109,7 +115,9 @@ namespace storm { } else { newVariable = this->manager->declareVariable(name, this->manager->getRationalType()); } - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -122,7 +130,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getIntegerType()); } solver->add(expressionAdapter->translateExpression((newVariable.getExpression() >= this->manager->rational(lowerBound)) && (newVariable.getExpression() <= this->manager->rational(upperBound)))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -135,7 +145,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getIntegerType()); } solver->add(expressionAdapter->translateExpression(newVariable.getExpression() >= this->manager->rational(lowerBound))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -148,7 +160,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getIntegerType()); } solver->add(expressionAdapter->translateExpression(newVariable.getExpression() <= this->manager->rational(upperBound))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -160,7 +174,9 @@ namespace storm { } else { newVariable = this->manager->declareVariable(name, this->manager->getIntegerType()); } - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -173,7 +189,9 @@ namespace storm { newVariable = this->manager->declareVariable(name, this->manager->getIntegerType()); } solver->add(expressionAdapter->translateExpression((newVariable.getExpression() >= this->manager->rational(storm::utility::zero())) && (newVariable.getExpression() <= this->manager->rational(storm::utility::one())))); - optimizationFunction = optimizationFunction + this->manager->rational(objectiveFunctionCoefficient) * newVariable; + if (!storm::utility::isZero(objectiveFunctionCoefficient)) { + optimizationSummands.push_back(this->manager->rational(objectiveFunctionCoefficient) * newVariable); + } return newVariable; } @@ -193,6 +211,7 @@ namespace storm { solver->push(); // Solve the optimization problem depending on the optimization direction + storm::expressions::Expression optimizationFunction = storm::expressions::sum(optimizationSummands); z3::optimize::handle optFuncHandle = this->getOptimizationDirection() == OptimizationDirection::Minimize ? solver->minimize(expressionAdapter->translateExpression(optimizationFunction)) : solver->maximize(expressionAdapter->translateExpression(optimizationFunction)); z3::check_result chkRes = solver->check(); STORM_LOG_THROW(chkRes != z3::unknown, storm::exceptions::InvalidStateException, "Unable to solve LP problem with Z3: Check result is unknown."); @@ -300,20 +319,40 @@ namespace storm { template void Z3LpSolver::writeModelToFile(std::string const& filename) const { + std::ofstream stream; + storm::utility::openFile(filename, stream); + stream << Z3_optimize_to_string(*context, *solver); + storm::utility::closeFile(stream); STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Exporting LP Problems to a file is not implemented for z3."); } template void Z3LpSolver::push() { + incrementaOptimizationSummandIndicators.push_back(optimizationSummands.size()); solver->push(); } template void Z3LpSolver::pop() { + STORM_LOG_ASSERT(!incrementaOptimizationSummandIndicators.empty(), "Tried to pop() without push()ing first."); solver->pop(); + // Delete summands of the optimization function that have been added since the last call to push() + optimizationSummands.resize(incrementaOptimizationSummandIndicators.back()); + incrementaOptimizationSummandIndicators.pop_back(); isIncremental = true; } - + + template + void Z3LpSolver::setMaximalMILPGap(ValueType const&, bool) { + // Since the solver is always exact, setting a gap has no effect. + // Intentionally left empty. + } + + template + ValueType Z3LpSolver::getMILPGap(bool relative) const { + // Since the solver is precise, the milp gap is always zero. + return storm::utility::zero(); + } #else template Z3LpSolver::Z3LpSolver(std::string const&, OptimizationDirection const&) { @@ -453,6 +492,16 @@ namespace storm { void Z3LpSolver::pop() { throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without Z3 or the version of Z3 does not support optimization. Yet, a method was called that requires this support."; } + + template + void Z3LpSolver::setMaximalMILPGap(ValueType const& gap, bool relative) { + throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without Z3 or the version of Z3 does not support optimization. Yet, a method was called that requires this support."; + } + + template + ValueType Z3LpSolver::getMILPGap(bool relative) const { + throw storm::exceptions::NotImplementedException() << "This version of storm was compiled without Z3 or the version of Z3 does not support optimization. Yet, a method was called that requires this support."; + } #endif template class Z3LpSolver; diff --git a/src/storm/solver/Z3LpSolver.h b/src/storm/solver/Z3LpSolver.h index ece5951f0..800d8fe79 100644 --- a/src/storm/solver/Z3LpSolver.h +++ b/src/storm/solver/Z3LpSolver.h @@ -99,6 +99,9 @@ namespace storm { virtual void push() override; virtual void pop() override; + virtual void setMaximalMILPGap(ValueType const& gap, bool relative) override; + virtual ValueType getMILPGap(bool relative) const override; + private: virtual storm::expressions::Expression getValue(storm::expressions::Variable const& variable) const; @@ -117,14 +120,17 @@ namespace storm { mutable std::unique_ptr lastCheckObjectiveValue; mutable std::unique_ptr lastCheckModel; - // Stores whether this solver is used in an incremental way (with push() and pop()) - bool isIncremental; - // An expression adapter that is used for translating the expression into Z3's format. std::unique_ptr expressionAdapter; - // The function that is to be optimized - storm::expressions::Expression optimizationFunction; + // The function that is to be optimized (we interpret this as a sum) + std::vector optimizationSummands; + + // Stores whether this solver is used in an incremental way (with push() and pop()) + bool isIncremental; + + // Stores the number of optimization summands at each call of push(). + std::vector incrementaOptimizationSummandIndicators; #endif };