diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp index f3d2555c9..337768578 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.cpp @@ -12,7 +12,9 @@ #include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/utility/export.h" +#include "storm/utility/solver.h" +#include "storm/exceptions/UnexpectedException.h" #include "storm/exceptions/InvalidOperationException.h" @@ -21,12 +23,12 @@ namespace storm { namespace multiobjective { template - DeterministicParetoExplorer::Point::Point(std::vector const& coordinates) : coordinates(coordinates), paretoOptimal(false) { + DeterministicParetoExplorer::Point::Point(std::vector const& coordinates) : coordinates(coordinates), paretoOptimal(false), onFacet(false) { STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported"); } template - DeterministicParetoExplorer::Point::Point(std::vector&& coordinates) : coordinates(std::move(coordinates)), paretoOptimal(false) { + DeterministicParetoExplorer::Point::Point(std::vector&& coordinates) : coordinates(std::move(coordinates)), paretoOptimal(false), onFacet(false) { STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported"); } @@ -92,6 +94,16 @@ namespace storm { return paretoOptimal; } + template + void DeterministicParetoExplorer::Point::setOnFacet(bool value) { + onFacet = value; + } + + template + bool DeterministicParetoExplorer::Point::liesOnFacet() const { + return onFacet; + } + template std::string DeterministicParetoExplorer::Point::toString(bool convertToDouble) const { std::stringstream out; @@ -147,7 +159,12 @@ namespace storm { STORM_LOG_WARN("Potential precision issues: Found a point that dominates another point which was flagged as pareto optimal. Distance of points is " << std::sqrt(storm::utility::convertNumber(storm::storage::geometry::squaredEuclideanDistance(pointsIt->second.get(), point.get())))); point.setParetoOptimal(true); } - pointsIt = points.erase(pointsIt); + if (pointsIt->second.liesOnFacet()) { + // Do not erase points that lie on a facet + ++pointsIt; + } else { + pointsIt = points.erase(pointsIt); + } break; case Point::DominanceResult::Dominated: // The new point is dominated by another point. @@ -156,6 +173,9 @@ namespace storm { if (point.isParetoOptimal()) { pointsIt->second.setParetoOptimal(); } + if (point.liesOnFacet()) { + pointsIt->second.setOnFacet(); + } return pointsIt->first; } } @@ -234,8 +254,17 @@ namespace storm { } template - void DeterministicParetoExplorer::Facet::addPoint(PointId const& pointId) { + void DeterministicParetoExplorer::Facet::addPoint(PointId const& pointId, Point const& point) { inducedSimplex = nullptr; + GeometryValueType product = storm::utility::vector::dotProduct(getHalfspace().normalVector(), point.get()); + if (product != getHalfspace().offset()) { + if (product < getHalfspace().offset()) { + STORM_LOG_DEBUG("The point on the facet actually has distance " << storm::utility::convertNumber(getHalfspace().euclideanDistance(point.get()))); + } else { + STORM_LOG_DEBUG("Halfspace of facet is shifted by " << storm::utility::convertNumber(getHalfspace().euclideanDistance(point.get())) << " to capture all points that are supposed to lie on the facet."); + halfspace.offset() = product; + } + } paretoPointsOnFacet.push_back(pointId); } @@ -282,6 +311,11 @@ namespace storm { return inducedSimplex; } + template + DeterministicParetoExplorer::FacetAnalysisContext::FacetAnalysisContext(Facet& f) : facet(f) { + // Intentionally left empty + } + template DeterministicParetoExplorer::DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) : model(preprocessorResult.preprocessedModel), objectives(preprocessorResult.objectives) { originalModelInitialState = *preprocessorResult.originalModel.getInitialStates().begin(); @@ -391,7 +425,7 @@ namespace storm { Point p(storm::utility::vector::convertNumericVector(*pIt)); if (last) { last = false; - p.setParetoOptimal(true); + p.setOnFacet(); // Adapt the overapproximation std::vector normalVector(objectives.size(), storm::utility::zero()); normalVector[objIndex] = storm::utility::one(); @@ -406,7 +440,7 @@ namespace storm { Facet f(std::move(h)); for (auto const& p : pointset) { if (f.getHalfspace().isPointOnBoundary(p.second.get())) { - f.addPoint(p.first); + f.addPoint(p.first, p.second); } } STORM_LOG_ASSERT(std::count(f.getHalfspace().normalVector().begin(), f.getHalfspace().normalVector().end(), storm::utility::zero()) + f.getNumberOfPoints() == objectives.size(), "Unexpected number of points on facet."); @@ -441,17 +475,18 @@ namespace storm { if (optimizeAndSplitFacet(env,f)) { return; } - std::set collectedPoints; - if (findAndCheckCachedPoints(env, f, collectedPoints)) { + FacetAnalysisContext context = createAnalysisContext(env, f); + + if (findAndCheckCachedPoints(env, context)) { return; } - if (analyzePointsOnFacet(env, f, collectedPoints)) { + if (analyzePointsOnFacet(env, context)) { return; } - if (analyzePointsInSimplex(env, f, collectedPoints)) { + if (analyzePointsInSimplex(env, context)) { return; } @@ -459,10 +494,47 @@ namespace storm { STORM_LOG_ERROR("Facet " << f.getHalfspace().toString(true) << " could not be analyzed completely."); } + template + typename DeterministicParetoExplorer::FacetAnalysisContext DeterministicParetoExplorer::createAnalysisContext(Environment const& env, Facet& f) { + FacetAnalysisContext res(f); + + res.expressionManager = std::make_shared(); + res.smtSolver = storm::utility::solver::SmtSolverFactory().create(*res.expressionManager); + + Polytope const& inducedPoly = res.facet.getInducedSimplex(pointset); + + res.x = inducedPoly->declareVariables(*res.expressionManager, "x"); + for (auto const& c : inducedPoly->getConstraints(*res.expressionManager, res.x)) { + res.smtSolver->add(c); + } + + res.xMinusEps = inducedPoly->declareVariables(*res.expressionManager, "y"); + for (auto const& c : inducedPoly->getConstraints(*res.expressionManager, res.xMinusEps)) { + res.smtSolver->add(c); + } + + auto eps = res.expressionManager->rational(env.modelchecker().multi().getPrecision()); + storm::expressions::Expression xme; + for (uint64_t i = 0; i < res.x.size(); ++i) { + storm::expressions::Expression subExpr = (res.xMinusEps[i].getExpression() == res.x[i].getExpression() - eps); + if (i == 0) { + xme = subExpr; + } else { + xme = xme && subExpr; + } + } + res.smtSolver->add(xme); + + res.smtSolver->push(); + + return res; + } + template bool DeterministicParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet& f) { // Obtain the correct weight vector auto weightVector = storm::utility::vector::convertNumericVector(f.getHalfspace().normalVector()); + bool weightVectorYieldsParetoOptimalPoint = !storm::utility::vector::hasZeroEntry(weightVector); for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { weightVector[objIndex] *= -storm::utility::one(); @@ -482,7 +554,8 @@ namespace storm { Point p(storm::utility::vector::convertNumericVector(*pIt)); if (last) { last = false; - p.setParetoOptimal(true); + p.setParetoOptimal(weightVectorYieldsParetoOptimalPoint); + p.setOnFacet(); addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), p); optPointId = pointset.addPoint(env, std::move(p)); } else { @@ -508,13 +581,13 @@ namespace storm { if (!storm::utility::vector::hasNegativeEntry(h.normalVector())) { STORM_LOG_ASSERT(h.isPointOnBoundary(optPoint.get()), "Unexpected facet found while splitting."); Facet fNew(std::move(h)); - fNew.addPoint(optPointId.get()); + fNew.addPoint(optPointId.get(), optPoint); auto vertexIt = vertices.begin(); ++vertexIt; for (auto const& pId : f.getPoints()) { assert(pointset.getPoint(pId).get() == *vertexIt); if (fNew.getHalfspace().isPointOnBoundary(*vertexIt)) { - fNew.addPoint(pId); + fNew.addPoint(pId, pointset.getPoint(pId)); } ++vertexIt; } @@ -533,20 +606,49 @@ namespace storm { } template - bool DeterministicParetoExplorer::findAndCheckCachedPoints(Environment const& env, Facet& f, std::set& collectedPoints) { - // TODO + bool DeterministicParetoExplorer::findAndCheckCachedPoints(Environment const& env, FacetAnalysisContext& context) { + Polytope inducedPoly = context.facet.getInducedSimplex(pointset); + pointset.collectPointsInPolytope(context.collectedPoints, inducedPoly); + if (context.collectedPoints.size() != context.facet.getNumberOfPoints()) { + STORM_LOG_ASSERT(context.collectedPoints.size() > context.facet.getNumberOfPoints(), "Did not find all points on the facet"); + // formulate SMT constraints to assert that xMinusEps is not dominated by one of the cached points + for (auto const& pId : context.collectedPoints) { + auto const& coordinates = pointset.getPoint(pId).get(); + storm::expressions::Expression pointAchievesXMinusEps; + for (uint64_t i = 0; i < coordinates.size(); ++i) { + storm::expressions::Expression subExpr = context.xMinusEps[i] <= context.expressionManager->rational(coordinates[i]); + if (i == 0) { + pointAchievesXMinusEps = subExpr; + } else { + pointAchievesXMinusEps = pointAchievesXMinusEps && subExpr; + } + } + context.smtSolver->add(!pointAchievesXMinusEps); + } + + context.smtSolver->push(); + + auto smtCheckResult = context.smtSolver->check(); + if (smtCheckResult == storm::solver::SmtSolver::CheckResult::Unsat) { + // There is no point x that violates the 'precision criterion' for this facet + return true; + } else { + STORM_LOG_THROW(smtCheckResult == storm::solver::SmtSolver::CheckResult::Sat, storm::exceptions::UnexpectedException, "The smt solver did not yield sat or unsat."); + return false; + } + } return false; } template - bool DeterministicParetoExplorer::analyzePointsOnFacet(Environment const& env, Facet& f, std::set& collectedPoints) { + bool DeterministicParetoExplorer::analyzePointsOnFacet(Environment const& env, FacetAnalysisContext& context) { // TODO return false; } template - bool DeterministicParetoExplorer::analyzePointsInSimplex(Environment const& env, Facet& f, std::set& collectedPoints) { + bool DeterministicParetoExplorer::analyzePointsInSimplex(Environment const& env, FacetAnalysisContext& context) { // TODO return false; } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h index aeafc244c..1ac75e161 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h @@ -9,6 +9,8 @@ #include "storm/storage/geometry/Polytope.h" #include "storm/storage/geometry/Halfspace.h" #include "storm/modelchecker/results/CheckResult.h" +#include "storm/storage/expressions/ExpressionManager.h" +#include "storm/solver/SmtSolver.h" namespace storm { @@ -44,12 +46,15 @@ namespace storm { void setParetoOptimal(bool value = true); bool isParetoOptimal() const; + void setOnFacet(bool value = true); + bool liesOnFacet() const; std::string toString(bool convertToDouble = false) const; private: std::vector coordinates; bool paretoOptimal; + bool onFacet; }; @@ -102,7 +107,7 @@ namespace storm { Facet(storm::storage::geometry::Halfspace const& halfspace); Facet(storm::storage::geometry::Halfspace&& halfspace); storm::storage::geometry::Halfspace const& getHalfspace() const; - void addPoint(PointId const& pointId); + void addPoint(PointId const& pointId, Point const& point); std::vector const& getPoints() const; uint64_t getNumberOfPoints() const; @@ -121,6 +126,19 @@ namespace storm { Polytope inducedSimplex; }; + struct FacetAnalysisContext { + FacetAnalysisContext(Facet& f); + + Facet& facet; + std::set collectedPoints; + std::unique_ptr smtSolver; + std::shared_ptr expressionManager; + + // Variables that encode two points that lie in the induced simplex of the analyzed facet + // xMinusEps = (x_1-eps, x_m-eps) + std::vector x, xMinusEps; + }; + DeterministicParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult); @@ -170,6 +188,8 @@ namespace storm { */ void processFacet(Environment const& env, Facet& f); + FacetAnalysisContext createAnalysisContext(Environment const& env, Facet& f); + /*! * Optimizes in the facet direction. If this results in a point that does not lie on the facet, * 1. The new Pareto optimal point is added @@ -182,7 +202,7 @@ namespace storm { * Finds all points that lie within the induced Simplex of the given facet. * Returns true if the facet is sufficiently precise when considering all added points */ - bool findAndCheckCachedPoints(Environment const& env, Facet& f, std::set& collectedPoints); + bool findAndCheckCachedPoints(Environment const& env, FacetAnalysisContext& context); /*! * Finds points that lie on the facet @@ -192,9 +212,9 @@ namespace storm { * use smt to find an eps-box in the simplex that does not contain a point. Add new points until one in the box is found. repeat. * stop when no more points or boxes can be found. */ - bool analyzePointsOnFacet(Environment const& env, Facet& f, std::set& collectedPoints); + bool analyzePointsOnFacet(Environment const& env, FacetAnalysisContext& context); - bool analyzePointsInSimplex(Environment const& env, Facet& f, std::set& collectedPoints); + bool analyzePointsInSimplex(Environment const& env, FacetAnalysisContext& context); Pointset pointset; std::queue unprocessedFacets;