diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index da5815cc4..fac4052a5 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -50,23 +50,33 @@ namespace storm { } template - std::vector DeterministicSchedsLpChecker::check(storm::Environment const& env) { + boost::optional> DeterministicSchedsLpChecker::check(storm::Environment const& env, Polytope area) { STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); STORM_LOG_TRACE("Checking a vertex..."); - swCheck.start(); swLpSolve.start(); swCheckVertices.start(); + swCheck.start(); + swLpBuild.start(); + lpModel->push(); + auto areaConstraints = area->getConstraints(lpModel->getManager(), currentObjectiveVariables); + for (auto const& c : areaConstraints) { + lpModel->addConstraint("", c); + } + lpModel->update(); + swLpBuild.stop(); + swLpSolve.start(); swCheckVertices.start(); lpModel->optimize(); swCheckVertices.stop(); swLpSolve.stop(); STORM_LOG_TRACE("\t Done checking a vertex..."); - STORM_LOG_ASSERT(!lpModel->isInfeasible(), "LP result is infeasable."); - STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); - - Point newPoint; - for (auto const& objVar : currentObjectiveVariables) { - newPoint.push_back(storm::utility::convertNumber(lpModel->getContinuousValue(objVar))); + boost::optional result; + if (!lpModel->isInfeasible()) { + STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); + result = Point(); + for (auto const& objVar : currentObjectiveVariables) { + result->push_back(storm::utility::convertNumber(lpModel->getContinuousValue(objVar))); + } } - + lpModel->pop(); swCheck.stop(); - return newPoint; + return result; } template diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h index a1352e3d1..adf3079fd 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h @@ -36,7 +36,7 @@ namespace storm { * Optimizes in the currently given direction * @return some optimal point found in that direction. */ - Point check(storm::Environment const& env); + boost::optional check(storm::Environment const& env, Polytope area); /*! * Optimizes in the currently given direction, recursively checks for points in the given area. diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp index 034fbe6fe..2c7ca7ec5 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp @@ -310,7 +310,6 @@ namespace storm { template std::unique_ptr DeterministicSchedsParetoExplorer::check(Environment const& env) { - clean(); initializeFacets(env); while (!unprocessedFacets.empty()) { @@ -324,8 +323,7 @@ namespace storm { for (auto const& p : pointset) { paretoPoints.push_back(storm::utility::vector::convertNumericVector(transformObjectiveValuesToOriginal(objectives, p.second.get()))); } - return std::make_unique>(originalModelInitialState, std::move(paretoPoints), - nullptr, nullptr); + return std::make_unique>(originalModelInitialState, std::move(paretoPoints), nullptr, nullptr); } template @@ -337,8 +335,7 @@ namespace storm { } template - void DeterministicSchedsParetoExplorer::addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace) { - GeometryValueType offset = storm::utility::vector::dotProduct(normalVector, pointOnHalfspace.get()); + void DeterministicSchedsParetoExplorer::addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, GeometryValueType const& offset) { if (env.modelchecker().multi().isPrintResultsSet()) { std::cout << "## unachievable halfspace: ["; bool first = true; @@ -391,29 +388,46 @@ namespace storm { } template - void DeterministicSchedsParetoExplorer::initializeFacets(Environment const& env) { + typename DeterministicSchedsParetoExplorer::Polytope DeterministicSchedsParetoExplorer::negateMinObjectives(Polytope const& polytope) const { + std::vector zeroRow(objectives.size(), storm::utility::zero()); + std::vector> transformationMatrix(objectives.size(), zeroRow); for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { - std::vector weightVector(objectives.size(), storm::utility::zero()); if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { - weightVector[objIndex] = -storm::utility::one(); + transformationMatrix[objIndex][objIndex] = -storm::utility::one(); } else { - weightVector[objIndex] = storm::utility::one(); + transformationMatrix[objIndex][objIndex] = storm::utility::one(); } - lpChecker->setCurrentWeightVector(weightVector); - auto point = lpChecker->check(env); - for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { - if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { - point[objIndex] *= -storm::utility::one(); - } - Point p(storm::utility::vector::convertNumericVector(point)); - p.setOnFacet(); - // Adapt the overapproximation - std::vector normalVector(objectives.size(), storm::utility::zero()); - normalVector[objIndex] = storm::utility::one(); - addHalfspaceToOverApproximation(env, normalVector, p); - pointset.addPoint(env, std::move(p)); + } + return polytope->affineTransformation(transformationMatrix, zeroRow); + } + + template + void DeterministicSchedsParetoExplorer::negateMinObjectives(std::vector& vector) const { + for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { + vector[objIndex] *= -storm::utility::one(); } } + } + + template + void DeterministicSchedsParetoExplorer::initializeFacets(Environment const& env) { + for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { + std::vector weightVector(objectives.size(), storm::utility::zero()); + weightVector[objIndex] = storm::utility::one(); + negateMinObjectives(weightVector); + lpChecker->setCurrentWeightVector(weightVector); + auto point = lpChecker->check(env, negateMinObjectives(this->overApproximation)); + STORM_LOG_THROW(point.is_initialized(), storm::exceptions::UnexpectedException, "Unable to find a point in the current overapproximation."); + negateMinObjectives(weightVector); + negateMinObjectives(point.get()); + Point p(point.get()); + p.setOnFacet(); + // Adapt the overapproximation + GeometryValueType offset = storm::utility::vector::dotProduct(weightVector, p.get()); + addHalfspaceToOverApproximation(env, weightVector, offset); + pointset.addPoint(env, std::move(p)); + } auto initialHalfspaces = pointset.downwardClosure()->getHalfspaces(); for (auto& h : initialHalfspaces) { @@ -424,9 +438,8 @@ namespace storm { } } 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."); - if (!checkFacetPrecision(env, f)) { - unprocessedFacets.push(std::move(f)); - } + + unprocessedFacets.push(std::move(f)); } } @@ -434,40 +447,18 @@ namespace storm { std::vector DeterministicSchedsParetoExplorer::getReferenceCoordinates() const { std::vector result; for (auto const& obj : objectives) { - ModelValueType value = storm::solver::minimize(obj.formula->getOptimalityType()) ? obj.upperResultBound.get() : obj.lowerResultBound.get(); + // TODO: use objectiveHelper here. + ModelValueType value = storm::solver::minimize(obj.formula->getOptimalityType()) ? -obj.upperResultBound.get() : obj.lowerResultBound.get(); result.push_back(storm::utility::convertNumber(value)); } return result; } - - template - bool DeterministicSchedsParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f) { - // TODO: - return false; - /* - auto const& inducedSimplex = f.getInducedSimplex(pointset); - - GeometryValueType eps = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); - // get a polytope that contains exactly the points y, such that y+eps is in the induced simplex - std::vector offsetVector(objectives.size(), -eps); - auto shiftedSimplex = inducedSimplex->shift(offsetVector); - - // If the intersection of both polytopes is empty, it means that there can not be a point y in the simplex - // such that y-eps is also in the simplex, i.e., the facet is already precise enough. - return inducedSimplex->intersection(shiftedSimplex)->isEmpty(); - */ - } - - template - bool DeterministicSchedsParetoExplorer::checkFacetPrecision(Environment const& env, Facet& f, std::set const& collectedSimplexPoints) { - assert(false); - return false; - } - template void DeterministicSchedsParetoExplorer::processFacet(Environment const& env, Facet& f) { - lpChecker->setCurrentWeightVector(f.getHalfspace().normalVector()); + std::vector weightVector = f.getHalfspace().normalVector(); + negateMinObjectives(weightVector); + lpChecker->setCurrentWeightVector(weightVector); if (optimizeAndSplitFacet(env,f)) { return; @@ -485,74 +476,32 @@ namespace storm { if (!polytopeTree.isEmpty()) { auto res = lpChecker->check(env, polytopeTree, eps); for (auto const& infeasableArea : res.second) { - addUnachievableArea(env, infeasableArea); + addUnachievableArea(env, negateMinObjectives(infeasableArea)); } for (auto& achievablePoint : res.first) { + negateMinObjectives(achievablePoint); pointset.addPoint(env, Point(std::move(achievablePoint))); } } - - } - - template - typename DeterministicSchedsParetoExplorer::FacetAnalysisContext DeterministicSchedsParetoExplorer::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); - */ - return res; } template bool DeterministicSchedsParetoExplorer::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 (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { - if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { - weightVector[objIndex] *= -storm::utility::one(); - } - } // Invoke optimization and insert the explored points boost::optional optPointId; - auto point = lpChecker->check(env); - for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { - if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) { - point[objIndex] *= -storm::utility::one(); - } + auto currentArea = negateMinObjectives(overApproximation->intersection(f.getHalfspace().invert())); + auto point = lpChecker->check(env, currentArea); + if (point.is_initialized()) { + negateMinObjectives(point.get()); + Point p(point.get()); + p.setOnFacet(); + GeometryValueType offset = storm::utility::vector::dotProduct(f.getHalfspace().normalVector(), p.get()); + addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), offset); + optPointId = pointset.addPoint(env, std::move(p)); + } else { + addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), f.getHalfspace().offset()); } - Point p(point); - p.setParetoOptimal(weightVectorYieldsParetoOptimalPoint); - p.setOnFacet(); - addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), p); - optPointId = pointset.addPoint(env, std::move(p)); // Potentially generate new facets if (optPointId) { @@ -584,9 +533,7 @@ namespace storm { ++vertexIt; } assert(vertexIt == vertices.end()); - if (!checkFacetPrecision(env, fNew)) { - unprocessedFacets.push(std::move(fNew)); - } + unprocessedFacets.push(std::move(fNew)); } } return true; @@ -596,36 +543,7 @@ namespace storm { return false; } } - - template - bool DeterministicSchedsParetoExplorer::addNewSimplexPoint(FacetAnalysisContext& context, PointId const& pointId, bool performCheck) { - auto const& coordinates = pointset.getPoint(pointId).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); - if (performCheck) { - auto smtCheckResult = context.smtSolver->check(); - if (smtCheckResult == storm::solver::SmtSolver::CheckResult::Unsat) { - // For all points x, there is a cached point that dominates or is equal to (x-eps). - // (we have a constraint pointAchievesXminusEps that does not not hold (double negation) - return true; - } else { - STORM_LOG_THROW(smtCheckResult == storm::solver::SmtSolver::CheckResult::Sat, storm::exceptions::UnexpectedException, "The smt solver did not yield sat or unsat."); - // there is a point x such that (x-eps) is not dominated by or equal to a cached point. - return false; - } - } else { - return false; - } - } - + template void DeterministicSchedsParetoExplorer::exportPlotOfCurrentApproximation(Environment const& env) { /* @@ -689,6 +607,7 @@ namespace storm { }; */ } + template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h index 0fe5dc408..c283d9aa2 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h @@ -154,7 +154,7 @@ namespace storm { /*! * Intersects the overapproximation with the given halfspace */ - void addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace); + void addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, GeometryValueType const& offset); /*! * Adds a polytope which consists of unachievable points @@ -173,26 +173,10 @@ namespace storm { std::vector getReferenceCoordinates() const; /*! - * Checks the precision of the given Facet and returns true, if no further processing of the facet is necessary - */ - bool checkFacetPrecision(Environment const& env, Facet& f); - - /*! - * Checks the precision of the given Facet and returns true, if no further processing of the facet is necessary. - * Also takes the given points within the simplex of the facet into account - */ - bool checkFacetPrecision(Environment const& env, Facet& f, std::set const& collectedSimplexPoints); - - /*! Processes the given facet as follows: - * 1. Optimize in the facet direction. Potentially, this adds new, unprocessed facets - * 2. Find points that have already been collected so far such that they lie in the induced simplex of the facet. - * 3. Find more points that lie on the facet - * 4. Find all points that lie in the induced simplex or prove that there are none + * Processes the given facet */ 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 @@ -201,14 +185,8 @@ namespace storm { */ bool optimizeAndSplitFacet(Environment const& env, Facet& f); - /*! - * Adds a new point that lies within the induced simplex of the given facet to the analysis context. - * @param context the analysis context - * @param pointId the id of the given point. - * @param performCheck if true, it is checked whether the facet is sufficiently precise now. If false, no check is performed. - * @return true iff performCheck is true and the facet is sufficiently precise. - */ - bool addNewSimplexPoint(FacetAnalysisContext& context, PointId const& pointId, bool performCheck); + Polytope negateMinObjectives(Polytope const& polytope) const; + void negateMinObjectives(std::vector& vector) const; Pointset pointset; std::queue unprocessedFacets;