diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index 8704a40c4..17f9ebd1c 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -603,41 +603,92 @@ namespace storm { if (polytopeTree.getChildren().empty()) { // At leaf nodes we need to perform the actual check. - STORM_LOG_TRACE("\tSolving MILP..."); - swCheckAreas.start(); - lpModel->optimize(); - swCheckAreas.stop(); - ++numLpQueries; - STORM_LOG_TRACE("\tDone solving MILP..."); - // STORM_PRINT_AND_LOG("Writing model to file '" << polytopeTree.toId() << ".lp'" << std::endl;); - // lpModel->writeModelToFile(polytopeTree.toId() + ".lp"); - - if (lpModel->isInfeasible()) { - infeasableAreas.push_back(polytopeTree.getPolytope()); - polytopeTree.clear(); - } else { - STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); - swValidate.start(); - Point newPoint = validateCurrentModel(env); - swValidate.stop(); - GeometryValueType offset = storm::utility::convertNumber(lpModel->getObjectiveValue()); - // Get the gap between the found solution and the known bound. - offset += storm::utility::convertNumber(lpModel->getMILPGap(false)); - // we might want to shift the halfspace to guarantee that our point is included. - offset = std::max(offset, storm::utility::vector::dotProduct(currentWeightVector, newPoint)); - auto halfspace = storm::storage::geometry::Halfspace(currentWeightVector, offset).invert(); - infeasableAreas.push_back(polytopeTree.getPolytope()->intersection(halfspace)); - if (infeasableAreas.back()->isEmpty()) { - infeasableAreas.pop_back(); - } - polytopeTree.setMinus(storm::storage::geometry::Polytope::create({halfspace})); - foundPoints.push_back(newPoint); - polytopeTree.substractDownwardClosure(newPoint, eps); - if (!polytopeTree.isEmpty()) { - checkRecursive(env, polytopeTree, eps, foundPoints, infeasableAreas, depth); + // Numerical instabilities might yield a point that is not actually inside the nodeConstraints. + // If the downward closure is disjoint from this tree node, we will not make further progress. + // In this case, we sharpen the violating constraints just a little bit. + // This way, valid solutions might be excluded, so technically, this can yield false negatives. + // However, since we apparently are dealing with numerical algorithms, we can't be sure about correctness anyway. + uint64_t num_sharpen = 0; + auto halfspaces = polytopeTree.getPolytope()->getHalfspaces(); + while (true) { + STORM_LOG_TRACE("\tSolving MILP..."); + swCheckAreas.start(); + lpModel->optimize(); + swCheckAreas.stop(); + ++numLpQueries; + STORM_LOG_TRACE("\tDone solving MILP..."); + + // STORM_PRINT_AND_LOG("Writing model to file '" << polytopeTree.toId() << ".lp'" << std::endl;); + // lpModel->writeModelToFile(polytopeTree.toId() + ".lp"); + + if (lpModel->isInfeasible()) { + infeasableAreas.push_back(polytopeTree.getPolytope()); + polytopeTree.clear(); + break; + } else { + STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); + swValidate.start(); + Point newPoint = validateCurrentModel(env); + swValidate.stop(); + // Check whether this new point yields any progress. + // There is no progress if (due to numerical inaccuracies) the downwardclosure (including points that are epsilon close to it) contained in this polytope. + // We multiply eps by 0.999 so that points that lie on the boundary of polytope and downw. do not count in the intersection. + Point newPointPlusEps = newPoint; + storm::utility::vector::addScaledVector(newPointPlusEps, eps, storm::utility::convertNumber(0.999)); + if (polytopeTree.getPolytope()->contains(newPoint) || !polytopeTree.getPolytope()->intersection(storm::storage::geometry::Polytope::createDownwardClosure({newPointPlusEps}))->isEmpty()) { + GeometryValueType offset = storm::utility::convertNumber(lpModel->getObjectiveValue()); + // Get the gap between the found solution and the known bound. + offset += storm::utility::convertNumber(lpModel->getMILPGap(false)); + // we might want to shift the halfspace to guarantee that our point is included. + offset = std::max(offset, storm::utility::vector::dotProduct(currentWeightVector, newPoint)); + auto halfspace = storm::storage::geometry::Halfspace(currentWeightVector, offset).invert(); + infeasableAreas.push_back(polytopeTree.getPolytope()->intersection(halfspace)); + if (infeasableAreas.back()->isEmpty()) { + infeasableAreas.pop_back(); + } + polytopeTree.setMinus(storm::storage::geometry::Polytope::create({halfspace})); + foundPoints.push_back(newPoint); + polytopeTree.substractDownwardClosure(newPoint, eps); + if (!polytopeTree.isEmpty()) { + checkRecursive(env, polytopeTree, eps, foundPoints, infeasableAreas, depth); + } + break; + } else { + // If we end up here, we have to sharpen the violated constraints for this polytope + for (auto& h : halfspaces) { + GeometryValueType distance = h.distance(newPoint); + // Check if the found point is outside of this halfspace + if (!storm::utility::isZero(distance)) { + // The issue has to be for some normal vector with a negative entry. Otherwise, the intersection with the downward closure wouldn't be empty + bool normalVectorContainsNegative = false; + for (auto const& hi : h.normalVector()) { + if (hi < storm::utility::zero()) { + normalVectorContainsNegative = true; + break; + } + } + if (normalVectorContainsNegative) { + h.offset() -= distance / storm::utility::convertNumber(2); + if (num_sharpen == 0) { + lpModel->push(); + } + lpModel->addConstraint("", h.toExpression(lpModel->getManager(), currentObjectiveVariables)); + ++num_sharpen; + break; + } + } + } + STORM_LOG_TRACE("\tSharpened LP"); + } } } + if (num_sharpen > 0) { + STORM_LOG_WARN("Numerical instabilities detected: LP Solver found an achievable point outside of the search area. The search area had to be sharpened " << num_sharpen << " times."); + // pop sharpened constraints + lpModel->pop(); + } + } else { // Traverse all the non-empty children. for (uint64_t childId = 0; childId < polytopeTree.getChildren().size(); ++childId) { @@ -709,4 +760,4 @@ namespace storm { template class DeterministicSchedsLpChecker, storm::RationalNumber>; } } -} \ No newline at end of file +} diff --git a/src/storm/storage/geometry/Halfspace.h b/src/storm/storage/geometry/Halfspace.h index fcd421983..b8a1aa7fc 100644 --- a/src/storm/storage/geometry/Halfspace.h +++ b/src/storm/storage/geometry/Halfspace.h @@ -5,6 +5,8 @@ #include #include "storm/utility/constants.h" #include "storm/utility/vector.h" +#include "storm/storage/expressions/Expressions.h" +#include "storm/storage/expressions/ExpressionManager.h" namespace storm { namespace storage { @@ -98,6 +100,16 @@ namespace storm { return stream.str(); } + storm::expressions::Expression toExpression(storm::expressions::ExpressionManager const& manager, std::vector const& variables) { + STORM_LOG_ASSERT(variables.size() == normalVector().size(), "Dimension missmatch."); + STORM_LOG_ASSERT(normalVector().size() != 0, "Invalid dimension."); + storm::expressions::Expression lhs = manager.rational(normalVector()[0]) * variables[0].getExpression(); + for (uint64_t dim = 1; dim < normalVector().size(); ++dim) { + lhs = lhs + manager.rational(normalVector()[dim]) * variables[dim].getExpression(); + } + return lhs <= manager.rational(offset()); + } + std::vector const& normalVector() const { return mNormalVector; }