LPSolvers: Allowing premature termination by specifying a mip gap. Fixes for incremental solving with Z3Lpsolver.

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<storm::settings::modules::DebugSettings>().isDebugSet() || storm::settings::getModule<storm::settings::modules::GlpkSettings>().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<typename ValueType>
@@ -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<std::pair<double, bool>*>(info);
+            double actualRelativeGap = glp_ios_mip_gap(t);
+            double factor = storm::utility::one<double>();
+            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<typename ValueType>
         void GlpkLpSolver<ValueType>::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<storm::settings::modules::GlpkSettings>().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<double, bool> 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<typename ValueType>
@@ -439,6 +471,24 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void GlpkLpSolver<ValueType>::setMaximalMILPGap(ValueType const& gap, bool relative) {
+            if (relative) {
+                this->maxMILPGap = storm::utility::convertNumber<double>(gap);
+            } else {
+                this->maxMILPGap = storm::utility::convertNumber<double>(gap);
+            }
+        }
+        
+        template<typename ValueType>
+        ValueType GlpkLpSolver<ValueType>::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<ValueType>(this->actualRelativeMILPGap * getObjectiveValue());
+            }
+        }
         
         template class GlpkLpSolver<double>;
         template class GlpkLpSolver<storm::RationalNumber>;

src/storm/solver/GlpkLpSolver.h

@@ -96,6 +96,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:
             /*!
              * Adds a variable with the given name, type, lower and upper bound and objective function coefficient.
@@ -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<storm::expressions::Variable> variables;
                 int firstConstraintIndex;
@@ -241,6 +248,14 @@ namespace storm {
                 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
     }

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:
             /*!

src/storm/solver/LpSolver.h

@@ -295,6 +295,20 @@ namespace storm {
              */
 			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.
             std::shared_ptr<storm::expressions::ExpressionManager> manager;

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<z3::context>(new z3::context(config));
             solver = std::unique_ptr<z3::optimize>(new z3::optimize(*context));
             expressionAdapter = std::unique_ptr<storm::adapters::Z3ExpressionAdapter>(new storm::adapters::Z3ExpressionAdapter(*this->manager, *context));
-            optimizationFunction = this->getManager().rational(storm::utility::zero<storm::RationalNumber>());
         }
         
         template<typename ValueType>
@@ -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<storm::RationalNumber>())) && (newVariable.getExpression() <= this->manager->rational(storm::utility::one<storm::RationalNumber>()))));
-            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<typename ValueType>
         void Z3LpSolver<ValueType>::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<typename ValueType>
         void Z3LpSolver<ValueType>::push()  {
+            incrementaOptimizationSummandIndicators.push_back(optimizationSummands.size());
             solver->push();
         }
         
         template<typename ValueType>
         void Z3LpSolver<ValueType>::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<typename ValueType>
+        void Z3LpSolver<ValueType>::setMaximalMILPGap(ValueType const&, bool) {
+            // Since the solver is always exact, setting a gap has no effect.
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        ValueType Z3LpSolver<ValueType>::getMILPGap(bool relative) const {
+            // Since the solver is precise, the milp gap is always zero.
+            return storm::utility::zero<ValueType>();
+        }
 #else
         template<typename ValueType>
         Z3LpSolver<ValueType>::Z3LpSolver(std::string const&, OptimizationDirection const&) {
@@ -453,6 +492,16 @@ namespace storm {
         void Z3LpSolver<ValueType>::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<typename ValueType>
+        void Z3LpSolver<ValueType>::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<typename ValueType>
+        ValueType Z3LpSolver<ValueType>::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<double>;

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<z3::expr> lastCheckObjectiveValue;
             mutable std::unique_ptr<z3::model> 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<storm::adapters::Z3ExpressionAdapter> 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<storm::expressions::Expression> 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<uint64_t> incrementaOptimizationSummandIndicators;
 
 #endif
         };