From d3052c725b33e19336eba877a7ee9a149f1a0baf Mon Sep 17 00:00:00 2001
From: Tim Quatmann <tim.quatmann@cs.rwth-aachen.de>
Date: Mon, 6 May 2019 10:55:36 +0200
Subject: [PATCH] DeterministicSchedsLpChecker: Make use of added Gurobi
 features: Specify a goal-gap and make use of solutions found during solving.

---
 .../DeterministicSchedsLpChecker.cpp          | 46 +++++++++++++++++--
 .../DeterministicSchedsLpChecker.h            |  2 +
 2 files changed, 45 insertions(+), 3 deletions(-)

diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp
index 66eabf8db..842580ef6 100644
--- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp
+++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp
@@ -92,6 +92,17 @@ namespace storm {
                     return {{}, {}};
                 }
                 
+                if (gurobiLpModel) {
+                    // For gurobi, it is possible to specify a gap between the obtained lower/upper objective bounds.
+                    GeometryValueType milpGap = storm::utility::zero<GeometryValueType>();
+                    for (auto const& wi : currentWeightVector) {
+                        milpGap += wi;
+                    }
+                    milpGap *= eps;
+                    gurobiLpModel->setMaximalMILPGap(storm::utility::convertNumber<ValueType>(milpGap), false);
+                    gurobiLpModel->update();
+                }
+                
                 std::vector<Point> foundPoints;
                 std::vector<Polytope> infeasableAreas;
                 checkRecursive(polytopeTree, eps, foundPoints, infeasableAreas);
@@ -105,6 +116,8 @@ namespace storm {
             void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::initializeLpModel(Environment const& env) {
                 uint64_t numStates = model.getNumberOfStates();
                 lpModel = storm::utility::solver::getLpSolver<ValueType>("model");
+                gurobiLpModel = dynamic_cast<storm::solver::GurobiLpSolver<ValueType>*>(lpModel.get());
+                
                 lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize);
                 initialStateResults.clear();
                 
@@ -176,6 +189,7 @@ namespace storm {
                                 storm::expressions::Expression upperValueBoundAtChoice = lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state));
                                 if (objectiveHelper[objIndex].minimizing()) {
                                     lpModel->addConstraint("", choiceValVar + (upperValueBoundAtChoice * (one - choiceVars[globalChoiceIndex])) >= choiceValue);
+                                    // Optional: lpModel->addConstraint("", choiceValVar <= choiceValue);
                                 } else {
                                     lpModel->addConstraint("", choiceValVar <= choiceValue);
                                     lpModel->addConstraint("", choiceValVar <= upperValueBoundAtChoice * choiceVars[globalChoiceIndex]);
@@ -238,15 +252,41 @@ namespace storm {
                         for (auto const& objVar : currentObjectiveVariables) {
                             newPoint.push_back(storm::utility::convertNumber<GeometryValueType>(lpModel->getContinuousValue(objVar)));
                         }
-                        auto halfspace = storm::storage::geometry::Halfspace<GeometryValueType>(currentWeightVector, storm::utility::vector::dotProduct(currentWeightVector, newPoint)).invert();
+                        std::vector<Point> newPoints = {newPoint};
+                        if (gurobiLpModel) {
+                            // gurobi might have other good solutions.
+                            for (uint64_t solutionIndex = 0; solutionIndex < gurobiLpModel->getSolutionCount(); ++ solutionIndex) {
+                                Point p;
+                                bool considerP = false;
+                                for (uint64_t objIndex = 0; objIndex < currentObjectiveVariables.size(); ++objIndex) {
+                                    p.push_back(storm::utility::convertNumber<GeometryValueType>(gurobiLpModel->getContinuousValue(currentObjectiveVariables[objIndex], solutionIndex)));
+                                    if (p.back() > newPoint[objIndex] + eps / storm::utility::convertNumber<GeometryValueType>(2)) {
+                                        // The other solution dominates the newPoint in this dimension and is also not too close to the newPoint.
+                                        considerP = true;
+                                    }
+                                }
+                                if (considerP) {
+                                    newPoints.push_back(std::move(p));
+                                }
+                            }
+                            std::cout << "found " << (newPoints.size() - 1) << " useful points" << std::endl;
+                        }
+                        GeometryValueType offset = storm::utility::convertNumber<GeometryValueType>(lpModel->getObjectiveValue());
+                        if (gurobiLpModel) {
+                            // Gurobi gives us the gap between the found solution and the known bound.
+                            offset += storm::utility::convertNumber<GeometryValueType>(gurobiLpModel->getMILPGap(false));
+                        }
+                        auto halfspace = storm::storage::geometry::Halfspace<GeometryValueType>(currentWeightVector, offset).invert();
                         infeasableAreas.push_back(polytopeTree.getPolytope()->intersection(halfspace));
                         if (infeasableAreas.back()->isEmpty()) {
                             infeasableAreas.pop_back();
                         }
                         swAux.start();
                         polytopeTree.setMinus(storm::storage::geometry::Polytope<GeometryValueType>::create({halfspace}));
-                        foundPoints.push_back(newPoint);
-                        polytopeTree.substractDownwardClosure(newPoint, eps);
+                        for (auto const& p : newPoints) {
+                            foundPoints.push_back(p);
+                            polytopeTree.substractDownwardClosure(p, eps);
+                        }
                         swAux.stop();
                         if (!polytopeTree.isEmpty()) {
                             checkRecursive(polytopeTree, eps, foundPoints, infeasableAreas);
diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h
index d6a75a194..e22b6d7a8 100644
--- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h
+++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h
@@ -7,6 +7,7 @@
 #include "storm/storage/geometry/Polytope.h"
 #include "storm/storage/geometry/PolytopeTree.h"
 #include "storm/solver/LpSolver.h"
+#include "storm/solver/GurobiLpSolver.h"
 #include "storm/utility/Stopwatch.h"
 
 namespace storm {
@@ -53,6 +54,7 @@ namespace storm {
                 std::vector<DeterministicSchedsObjectiveHelper<ModelType>> const& objectiveHelper;
 
                 std::unique_ptr<storm::solver::LpSolver<ValueType>> lpModel;
+                storm::solver::GurobiLpSolver<ValueType>* gurobiLpModel;
                 std::vector<storm::expressions::Expression> initialStateResults;
                 std::vector<storm::expressions::Variable> currentObjectiveVariables;
                 std::vector<GeometryValueType> currentWeightVector;