From ee090b630e90c8507d47843e492c7870ede33471 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Mon, 29 Apr 2019 07:21:22 +0200 Subject: [PATCH] deterministic schedulers: Refactored code for lp-based checker. --- .../DeterministicSchedsLpChecker.cpp | 253 +++++++ .../DeterministicSchedsLpChecker.h | 71 ++ .../DeterministicSchedsObjectiveHelper.cpp | 159 ++++ .../DeterministicSchedsObjectiveHelper.h | 44 ++ .../DeterministicSchedsParetoExplorer.cpp | 699 ++++++++++++++++++ .../DeterministicSchedsParetoExplorer.h | 226 ++++++ .../multiObjectiveModelChecking.cpp | 4 +- src/storm/storage/geometry/PolytopeTree.h | 142 ++++ 8 files changed, 1596 insertions(+), 2 deletions(-) create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h create mode 100644 src/storm/storage/geometry/PolytopeTree.h diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp new file mode 100644 index 000000000..f621a78b6 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -0,0 +1,253 @@ +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h" +#include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/utility/solver.h" + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + DeterministicSchedsLpChecker::DeterministicSchedsLpChecker(ModelType const& model, std::vector> const& objectives) : model(model) { + swInit.start(); + initializeObjectiveHelper(objectives); + initializeLpModel(); + swInit.stop(); + } + + template + DeterministicSchedsLpChecker::~DeterministicSchedsLpChecker() { + std::cout << "Deterministic Scheds LP CHECKER STATISTICS: " << std::endl; + std::cout << "\t" << swInit << " seconds for initialization" << std::endl; + std::cout << "\t" << swCheck << " seconds for checking, including" << std::endl; + std::cout << "\t\t" << swLpBuild << " seconds for LP building" << std::endl; + std::cout << "\t\t" << swLpSolve << " seconds for LP solving" << std::endl; + std::cout << "\t\t" << swCheckVertices << " seconds for checking the vertices of the convex hull." << std::endl; + std::cout << "\t" << swAux << " seconds for aux stuff" << std::endl; + } + + template + void DeterministicSchedsLpChecker::setCurrentWeightVector(std::vector const& weightVector) { + STORM_LOG_ASSERT(!weightVector.empty(), "Setting an empty weight vector is not supported."); + swLpBuild.start(); + if (!currentWeightVector.empty()) { + // Pop information of the current weight vector. + lpModel->pop(); + lpModel->update(); + currentObjectiveVariables.clear(); + } + currentWeightVector = weightVector; + + lpModel->push(); + // set up objective function for the given weight vector + for (uint64_t objIndex = 0; objIndex < initialStateResults.size(); ++objIndex) { + currentObjectiveVariables.push_back(lpModel->addUnboundedContinuousVariable("w_" + std::to_string(objIndex), storm::utility::convertNumber(weightVector[objIndex]))); + lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == initialStateResults[objIndex]); + } + lpModel->update(); + swLpBuild.stop(); + } + + template + std::vector DeterministicSchedsLpChecker::check(storm::Environment const& env) { + STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); + STORM_LOG_TRACE("Checking a vertex..."); + swCheck.start(); swCheckVertices.start(); swLpSolve.start(); + lpModel->optimize(); + swLpSolve.stop(); + 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))); + } + + swCheckVertices.stop(); swCheck.stop(); + return newPoint; + } + + template + std::pair>, std::vector>>> DeterministicSchedsLpChecker::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps) { + std::cout << "Checking " << polytopeTree.toString() << std::endl << "\t"; + swCheck.start(); + STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); + if (polytopeTree.isEmpty()) { + return {{}, {}}; + } + + std::vector foundPoints; + std::vector infeasableAreas; + checkRecursive(polytopeTree, eps, foundPoints, infeasableAreas); + + swCheck.stop(); + std::cout << " done!" << std::endl; + return {foundPoints, infeasableAreas}; + } + + template + void DeterministicSchedsLpChecker::initializeObjectiveHelper(std::vector> const& objectives) { + objectiveHelper.reserve(objectives.size()); + for (auto const& obj : objectives) { + objectiveHelper.emplace_back(model, obj); + } + } + + template + void DeterministicSchedsLpChecker::initializeLpModel() { + uint64_t numStates = model.getNumberOfStates(); + lpModel = storm::utility::solver::getLpSolver("model"); + lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + initialStateResults.clear(); + + auto one = lpModel->getConstant(storm::utility::one()); + + // Create choice variables and assert that at least one choice is taken at each state. + std::vector choiceVars; + choiceVars.reserve(model.getNumberOfChoices()); + for (uint64_t state = 0; state < numStates; ++state) { + uint64_t numChoices = model.getNumberOfChoices(state); + if (numChoices == 1) { + choiceVars.emplace_back(); + } else { + std::vector localChoices; + for (uint64_t choice = 0; choice < numChoices; ++choice) { + localChoices.push_back(lpModel->addBoundedIntegerVariable("c" + std::to_string(state) + "_" + std::to_string(choice), 0, 1).getExpression()); + } + lpModel->addConstraint("", storm::expressions::sum(localChoices).reduceNesting() >= one); + choiceVars.insert(choiceVars.end(), localChoices.begin(), localChoices.end()); + } + } + + for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { + auto const& schedulerIndependentStates = objectiveHelper[objIndex].getSchedulerIndependentStateValues(); + // Create state variables + std::vector stateVars; + stateVars.reserve(numStates); + for (uint64_t state = 0; state < numStates; ++state) { + auto valIt = schedulerIndependentStates.find(state); + if (valIt == schedulerIndependentStates.end()) { + stateVars.push_back(lpModel->addBoundedContinuousVariable("x" + std::to_string(objIndex) + "_" + std::to_string(state), objectiveHelper[objIndex].getLowerValueBoundAtState(state), objectiveHelper[objIndex].getUpperValueBoundAtState(state)).getExpression()); + } else { + stateVars.push_back(lpModel->getConstant(valIt->second)); + } + if (state == *model.getInitialStates().begin()) { + initialStateResults.push_back(stateVars.back()); + } + } + + // Create and assert choice values + auto const& choiceValueOffsets = objectiveHelper[objIndex].getChoiceValueOffsets(); + for (uint64_t state = 0; state < numStates; ++state) { + if (schedulerIndependentStates.find(state) != schedulerIndependentStates.end()) { + continue; + } + uint64_t numChoices = model.getNumberOfChoices(state); + uint64_t choiceOffset = model.getTransitionMatrix().getRowGroupIndices()[state]; + for (uint64_t choice = 0; choice < numChoices; ++choice) { + storm::expressions::Expression choiceValue; + auto valIt = choiceValueOffsets.find(choiceOffset + choice); + if (valIt != choiceValueOffsets.end()) { + choiceValue = lpModel->getConstant(valIt->second); + } + for (auto const& transition : model.getTransitionMatrix().getRow(state, choice)) { + storm::expressions::Expression transitionValue = lpModel->getConstant(transition.getValue()) * stateVars[transition.getColumn()]; + if (choiceValue.isInitialized()) { + choiceValue = choiceValue + transitionValue; + } else { + choiceValue = transitionValue; + } + } + choiceValue = choiceValue.simplify().reduceNesting(); + if (numChoices == 1) { + lpModel->addConstraint("", stateVars[state] == choiceValue); + } else { + uint64_t globalChoiceIndex = model.getTransitionMatrix().getRowGroupIndices()[state] + choice; + storm::expressions::Expression maxDiff = lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(state) - objectiveHelper[objIndex].getLowerValueBoundAtState(state)) * (one - choiceVars[globalChoiceIndex]); + lpModel->addConstraint("", stateVars[state] - choiceValue <= maxDiff); + lpModel->addConstraint("", choiceValue - stateVars[state] <= maxDiff); + } + } + } + } + lpModel->update(); + } + + template + void DeterministicSchedsLpChecker::checkRecursive(storm::storage::geometry::PolytopeTree & polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas) { + std::cout << "."; + std::cout.flush(); + STORM_LOG_ASSERT(!polytopeTree.isEmpty(), "Tree node is empty"); + STORM_LOG_ASSERT(!polytopeTree.getPolytope()->isEmpty(), "Tree node is empty."); + STORM_LOG_TRACE("Checking " << polytopeTree.toString()); + + swLpBuild.start(); + lpModel->push(); + // Assert the constraints of the current polytope + auto nodeConstraints = polytopeTree.getPolytope()->getConstraints(lpModel->getManager(), currentObjectiveVariables); + for (auto const& constr : nodeConstraints) { + lpModel->addConstraint("", constr); + } + lpModel->update(); + swLpBuild.stop(); + + if (polytopeTree.getChildren().empty()) { + // At leaf nodes we need to perform the actual check. + swLpSolve.start(); + STORM_LOG_TRACE("\tSolving MILP..."); + lpModel->optimize(); + STORM_LOG_TRACE("\tDone solving MILP..."); + swLpSolve.stop(); + + if (lpModel->isInfeasible()) { + infeasableAreas.push_back(polytopeTree.getPolytope()); + polytopeTree.clear(); + } else { + 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))); + } + auto halfspace = storm::storage::geometry::Halfspace(currentWeightVector, storm::utility::vector::dotProduct(currentWeightVector, newPoint)).invert(); + infeasableAreas.push_back(polytopeTree.getPolytope()->intersection(halfspace)); + if (infeasableAreas.back()->isEmpty()) { + infeasableAreas.pop_back(); + } + swAux.start(); + polytopeTree.setMinus(storm::storage::geometry::Polytope::create({halfspace})); + foundPoints.push_back(newPoint); + polytopeTree.substractDownwardClosure(newPoint, eps); + swAux.stop(); + if (!polytopeTree.isEmpty()) { + checkRecursive(polytopeTree, eps, foundPoints, infeasableAreas); + } + } + } else { + // Traverse all the children. + for (uint64_t childId = 0; childId < polytopeTree.getChildren().size(); ++childId) { + uint64_t newPointIndex = foundPoints.size(); + checkRecursive(polytopeTree.getChildren()[childId], eps, foundPoints, infeasableAreas); + STORM_LOG_ASSERT(polytopeTree.getChildren()[childId].isEmpty(), "expected empty children."); + // Make the new points known to the right siblings + for (; newPointIndex < foundPoints.size(); ++newPointIndex) { + for (uint64_t siblingId = childId + 1; siblingId < polytopeTree.getChildren().size(); ++siblingId) { + polytopeTree.getChildren()[siblingId].substractDownwardClosure(foundPoints[newPointIndex], eps); + } + } + } + // All children are empty now, so this node becomes empty. + polytopeTree.clear(); + } + swLpBuild.start(); + lpModel->pop(); + swLpBuild.stop(); + } + + template class DeterministicSchedsLpChecker, storm::RationalNumber>; + template class DeterministicSchedsLpChecker, storm::RationalNumber>; + template class DeterministicSchedsLpChecker, storm::RationalNumber>; + template class DeterministicSchedsLpChecker, storm::RationalNumber>; + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h new file mode 100644 index 000000000..3a88b0940 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h @@ -0,0 +1,71 @@ +#pragma once + +#include + +#include "storm/modelchecker/multiobjective/Objective.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h" +#include "storm/storage/geometry/Polytope.h" +#include "storm/storage/geometry/PolytopeTree.h" +#include "storm/solver/LpSolver.h" +#include "storm/utility/Stopwatch.h" + +namespace storm { + + class Environment; + + namespace modelchecker { + namespace multiobjective { + + template + class DeterministicSchedsLpChecker { + public: + + typedef typename ModelType::ValueType ValueType; + typedef typename std::shared_ptr> Polytope; + typedef typename std::vector Point; + + DeterministicSchedsLpChecker(ModelType const& model, std::vector> const& objectives); + ~DeterministicSchedsLpChecker(); + + /*! + * Specifies the current direction. + */ + void setCurrentWeightVector(std::vector const& weightVector); + + /*! + * Optimizes in the currently given direction + * @return some optimal point found in that direction. + */ + Point check(storm::Environment const& env); + + /*! + * Optimizes in the currently given direction, recursively checks for points in the given area. + * @return all pareto optimal points in the area given by polytopeTree as well as a set of area in which no solution lies (the points might be achievable via some point outside of this area, though) + */ + std::pair, std::vector> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps); + + private: + void initializeObjectiveHelper(std::vector> const& objectives); + void initializeLpModel(); + + void checkRecursive(storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas); + + ModelType const& model; + std::vector> objectiveHelper; + + std::unique_ptr> lpModel; + std::vector initialStateResults; + std::vector currentObjectiveVariables; + std::vector currentWeightVector; + + storm::utility::Stopwatch swInit; + storm::utility::Stopwatch swCheck; + storm::utility::Stopwatch swCheckVertices; + storm::utility::Stopwatch swLpSolve; + storm::utility::Stopwatch swLpBuild; + storm::utility::Stopwatch swAux; + }; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp new file mode 100644 index 000000000..67bc9d24d --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp @@ -0,0 +1,159 @@ +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h" + +#include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/storage/BitVector.h" +#include "storm/utility/graph.h" +#include "storm/utility/FilteredRewardModel.h" +#include "storm/logic/Formulas.h" + +#include "storm/exceptions/UnexpectedException.h" +namespace storm { + namespace modelchecker { + namespace multiobjective { + + + template + DeterministicSchedsObjectiveHelper::DeterministicSchedsObjectiveHelper(ModelType const& model, storm::modelchecker::multiobjective::Objective const& objective) : model(model), objective(objective) { + // Intentionally left empty + } + + template + storm::storage::BitVector evaluatePropositionalFormula(ModelType const& model, storm::logic::Formula const& formula) { + storm::modelchecker::SparsePropositionalModelChecker mc(model); + auto checkResult = mc.check(formula); + STORM_LOG_THROW(checkResult && checkResult->isExplicitQualitativeCheckResult(), storm::exceptions::UnexpectedException, "Unexpected type of check result for subformula " << formula << "."); + return checkResult->asExplicitQualitativeCheckResult().getTruthValuesVector(); + } + + template + std::map const& DeterministicSchedsObjectiveHelper::getSchedulerIndependentStateValues() const { + if (!schedulerIndependentStateValues) { + auto const& formula = *objective.formula; + std::map result; + if (formula.isProbabilityOperatorFormula() && formula.getSubformula().isUntilFormula()) { + storm::storage::BitVector phiStates = evaluatePropositionalFormula(model, formula.getSubformula().asUntilFormula().getLeftSubformula()); + storm::storage::BitVector psiStates = evaluatePropositionalFormula(model, formula.getSubformula().asUntilFormula().getRightSubformula()); + auto backwardTransitions = model.getBackwardTransitions(); + { + storm::storage::BitVector prob1States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), backwardTransitions, phiStates, psiStates); + for (auto const& prob1State : prob1States) { + result[prob1State] = storm::utility::one(); + } + } + { + storm::storage::BitVector prob0States = storm::utility::graph::performProb0A(backwardTransitions, phiStates, psiStates); + for (auto const& prob0State : prob0States) { + result[prob0State] = storm::utility::zero(); + } + } + } else if (formula.getSubformula().isEventuallyFormula() && (formula.isRewardOperatorFormula() || formula.isTimeOperatorFormula())) { + storm::storage::BitVector rew0States = evaluatePropositionalFormula(model, formula.getSubformula().asEventuallyFormula().getSubformula()); + if (formula.isRewardOperatorFormula()) { + auto const& baseRewardModel = formula.asRewardOperatorFormula().hasRewardModelName() ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()) : model.getUniqueRewardModel(); + auto rewardModel = storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), formula.getSubformula().asEventuallyFormula()); + storm::storage::BitVector statesWithoutReward = rewardModel.get().getStatesWithZeroReward(model.getTransitionMatrix()); + rew0States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), statesWithoutReward, rew0States); + } + for (auto const& rew0State : rew0States) { + result[rew0State] = storm::utility::zero(); + } + } else if (formula.isRewardOperatorFormula() && formula.getSubformula().isTotalRewardFormula()) { + auto const& baseRewardModel = formula.asRewardOperatorFormula().hasRewardModelName() ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()) : model.getUniqueRewardModel(); + auto rewardModel = storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), formula.getSubformula().asTotalRewardFormula()); + storm::storage::BitVector statesWithoutReward = rewardModel.get().getStatesWithZeroReward(model.getTransitionMatrix()); + storm::storage::BitVector rew0States = storm::utility::graph::performProbGreater0E(model.getBackwardTransitions(), statesWithoutReward, ~statesWithoutReward); + rew0States.complement(); + for (auto const& rew0State : rew0States) { + result[rew0State] = storm::utility::zero(); + } + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported."); + } + schedulerIndependentStateValues = std::move(result); + } + return schedulerIndependentStateValues.get(); + } + + template + std::map const& DeterministicSchedsObjectiveHelper::getChoiceValueOffsets() const { + if (!choiceValueOffsets) { + auto const& formula = *objective.formula; + auto const& subformula = formula.getSubformula(); + std::map result; + if (formula.isProbabilityOperatorFormula() && subformula.isUntilFormula()) { + // In this case, there is nothing to be done. + } else if (formula.isRewardOperatorFormula() && (subformula.isTotalRewardFormula() || subformula.isEventuallyFormula())) { + auto const& baseRewardModel = formula.asRewardOperatorFormula().hasRewardModelName() ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()) : model.getUniqueRewardModel(); + auto rewardModel = subformula.isEventuallyFormula() ? storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), subformula.asEventuallyFormula()) : storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), subformula.asTotalRewardFormula()); + std::vector choiceBasedRewards = rewardModel.get().getTotalRewardVector(model.getTransitionMatrix()); + // Set entries for all non-zero reward choices at states whose value is not already known. + // This relies on the fact that for goal states in reachability reward formulas, getSchedulerIndependentStateValues()[state] is set to zero. + auto const& rowGroupIndices = model.getTransitionMatrix().getRowGroupIndices(); + auto const& stateValues = getSchedulerIndependentStateValues(); + for (uint64_t state = 0; state < model.getNumberOfStates(); ++state) { + if (stateValues.find(state) == stateValues.end()) { + for (uint64_t choice = rowGroupIndices[state]; choice < rowGroupIndices[state + 1]; ++choice) { + if (!storm::utility::isZero(choiceBasedRewards[choice])) { + result[choice] = choiceBasedRewards[choice]; + } + } + } + } + } else if (formula.isTimeOperatorFormula() && subformula.isEventuallyFormula()) { + auto const& rowGroupIndices = model.getTransitionMatrix().getRowGroupIndices(); + auto const& stateValues = getSchedulerIndependentStateValues(); + std::vector const* rates = nullptr; + storm::storage::BitVector const* ms = nullptr; + if (model.isOfType(storm::models::ModelType::MarkovAutomaton)) { + auto ma = model.template as>(); + rates = &ma->getExitRates(); + ms = &ma->getMarkovianStates(); + } + if (model.isOfType(storm::models::ModelType::Mdp)) { + // Set all choice offsets to one, except for the ones at states in scheduerIndependentStateValues. + for (uint64_t state = 0; state < model.getNumberOfStates(); ++state) { + if (stateValues.find(state) == stateValues.end()) { + ValueType value = storm::utility::one(); + if (rates) { + if (ms->get(state)) { + value /= (*rates)[state]; + } else { + // Nothing to be done for probabilistic states + continue; + } + } + for (uint64_t choice = rowGroupIndices[state]; choice < rowGroupIndices[state + 1]; ++choice) { + result[choice] = value; + } + } + } + } + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported."); + } + choiceValueOffsets = std::move(result); + } + return choiceValueOffsets.get(); + } + + template + typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getUpperValueBoundAtState(uint64_t state) const{ + return objective.upperResultBound.get(); + } + + template + typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getLowerValueBoundAtState(uint64_t state) const{ + return objective.lowerResultBound.get(); + } + + template class DeterministicSchedsObjectiveHelper>; + template class DeterministicSchedsObjectiveHelper>; + template class DeterministicSchedsObjectiveHelper>; + template class DeterministicSchedsObjectiveHelper>; + } + } +} diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h new file mode 100644 index 000000000..5bb778c87 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h @@ -0,0 +1,44 @@ +#pragma once + +#include +#include + +#include "storm/modelchecker/multiobjective/Objective.h" + + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + class DeterministicSchedsObjectiveHelper { + public: + + typedef typename ModelType::ValueType ValueType; + DeterministicSchedsObjectiveHelper(ModelType const& model, Objective const& objective); + + /*! + * Returns states and values for states that are independent of the scheduler. + */ + std::map const& getSchedulerIndependentStateValues() const; + + /*! + * Returns offsets of each choice value (e.g., the reward) if non-zero. + * This does not include choices of states with independent state values + */ + std::map const& getChoiceValueOffsets() const; + + ValueType const& getUpperValueBoundAtState(uint64_t state) const; + ValueType const& getLowerValueBoundAtState(uint64_t state) const; + private: + + mutable boost::optional> schedulerIndependentStateValues; + mutable boost::optional> choiceValueOffsets; + + ModelType const& model; + Objective const& objective; + + }; + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp new file mode 100644 index 000000000..1b5662e99 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp @@ -0,0 +1,699 @@ +#include +#include + + +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h" +#include "storm/storage/geometry/coordinates.h" +#include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" +#include "storm/modelchecker/multiobjective/MultiObjectivePostprocessing.h" +#include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h" +#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" + + +namespace storm { + namespace modelchecker { + namespace multiobjective { + + template + DeterministicSchedsParetoExplorer::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 + DeterministicSchedsParetoExplorer::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"); + } + + template + std::vector& DeterministicSchedsParetoExplorer::Point::get() { + return coordinates; + } + + template + std::vector const& DeterministicSchedsParetoExplorer::Point::get() const { + return coordinates; + } + + template + uint64_t DeterministicSchedsParetoExplorer::Point::dimension() const { + STORM_LOG_ASSERT(!coordinates.empty(), "Points with dimension 0 are not supported"); + return coordinates.size(); + } + + template + typename DeterministicSchedsParetoExplorer::Point::DominanceResult DeterministicSchedsParetoExplorer::Point::getDominance(Point const& other) const { + STORM_LOG_ASSERT(this->dimension() == other.dimension(), "Non-Equal dimensions of points: [" << this->toString() << "] vs. [" << other.toString() << "]"); + auto thisIt = this->get().begin(); + auto otherIt = other.get().begin(); + auto thisItE = this->get().end(); + + // Find the first entry where the points differ + while (*thisIt == *otherIt) { + ++thisIt; + ++otherIt; + if (thisIt == thisItE) { + return DominanceResult::Equal; + } + } + + if (*thisIt > *otherIt) { + // *this might dominate other + for (++thisIt, ++otherIt; thisIt != thisItE; ++thisIt, ++otherIt) { + if (*thisIt < *otherIt) { + return DominanceResult::Incomparable; + } + } + return DominanceResult::Dominates; + } else { + assert(*thisIt < *otherIt); + // *this might be dominated by other + for (++thisIt, ++otherIt; thisIt != thisItE; ++thisIt, ++otherIt) { + if (*thisIt > *otherIt) { + return DominanceResult::Incomparable; + } + } + return DominanceResult::Dominated; + } + } + + template + void DeterministicSchedsParetoExplorer::Point::setParetoOptimal(bool value) { + paretoOptimal = value; + } + + template + bool DeterministicSchedsParetoExplorer::Point::isParetoOptimal() const { + return paretoOptimal; + } + + template + void DeterministicSchedsParetoExplorer::Point::setOnFacet(bool value) { + onFacet = value; + } + + template + bool DeterministicSchedsParetoExplorer::Point::liesOnFacet() const { + return onFacet; + } + + template + std::string DeterministicSchedsParetoExplorer::Point::toString(bool convertToDouble) const { + std::stringstream out; + bool first = true; + for (auto const& pi : this->get()) { + if (first) { + first = false; + } else { + out << ", "; + } + if (convertToDouble) { + out << storm::utility::convertNumber(pi); + } else { + out << pi; + } + } + return out.str(); + } + + // template + // bool operator<(typename DeterministicSchedsParetoExplorer::Point const& lhs, typename DeterministicSchedsParetoExplorer::Point const& rhs) { + // STORM_LOG_ASSERT(lhs.dimension() == rhs.dimension(), "Non-Equal dimensions of points: " << lhs << " vs. " << rhs); + // for (uint64_t i = 0; i < lhs.dimension(); ++i) { + // if (lhs.get()[i] < rhs.get()[i]) { + // return true; + // } else if (lhs.get()[i] != rhs.get()[i]) { + // return false; + // } + // } + // return false; + // } + + template + DeterministicSchedsParetoExplorer::Pointset::Pointset() : currId(1) { + // Intentionally left empty + } + + template + boost::optional::PointId> DeterministicSchedsParetoExplorer::Pointset::addPoint(Environment const& env, Point&& point) { + + // Find dominated and dominating points + auto pointsIt = points.begin(); + while (pointsIt != points.end()) { + switch (point.getDominance(pointsIt->second)) { + case Point::DominanceResult::Incomparable: + // Nothing to be done for this point + ++pointsIt; + break; + case Point::DominanceResult::Dominates: + // Found a point in the set that is dominated by the new point, so we erase it + if (pointsIt->second.isParetoOptimal()) { + 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); + } + 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. + return boost::none; + case Point::DominanceResult::Equal: + if (point.isParetoOptimal()) { + pointsIt->second.setParetoOptimal(); + } + if (point.liesOnFacet()) { + pointsIt->second.setOnFacet(); + } + return pointsIt->first; + } + } + + if (env.modelchecker().multi().isPrintResultsSet()) { + std::cout << "## achievable point: [" << point.toString(true) << "]" << std::endl; + } + + points.emplace_hint(points.end(), currId, std::move(point)); + return currId++; + } + + template + typename DeterministicSchedsParetoExplorer::Point const& DeterministicSchedsParetoExplorer::Pointset::getPoint(PointId const& id) const { + return points.at(id); + } + + template + typename DeterministicSchedsParetoExplorer::Pointset::iterator_type DeterministicSchedsParetoExplorer::Pointset::begin() const { + return points.begin(); + } + + template + typename DeterministicSchedsParetoExplorer::Pointset::iterator_type DeterministicSchedsParetoExplorer::Pointset::end() const { + return points.end(); + } + + template + uint64_t DeterministicSchedsParetoExplorer::Pointset::size() const { + return points.size(); + } + + template + typename DeterministicSchedsParetoExplorer::Polytope DeterministicSchedsParetoExplorer::Pointset::downwardClosure() const { + std::vector> pointsAsVector; + pointsAsVector.reserve(size()); + for (auto const& p : points) { + pointsAsVector.push_back(p.second.get()); + } + return storm::storage::geometry::Polytope::createDownwardClosure(std::move(pointsAsVector)); + } + + template + void DeterministicSchedsParetoExplorer::Pointset::collectPointsInPolytope(std::set& collectedPoints, Polytope const& polytope) { + for (auto const& p : points) { + if (polytope->contains(p.second.get())) { + collectedPoints.insert(p.first); + } + } + } + + template + void DeterministicSchedsParetoExplorer::Pointset::printToStream(std::ostream& out, bool includeIDs, bool convertToDouble) { + for (auto const& p : this->points) { + if (includeIDs) { + out << p.first << ": [" << p.second.toString(convertToDouble) << "]" << std::endl; + } else { + out << p.second.toString(convertToDouble) << std::endl; + } + } + } + + template + DeterministicSchedsParetoExplorer::Facet::Facet(storm::storage::geometry::Halfspace const& halfspace) : halfspace(halfspace) { + // Intentionally left empty + } + + template + DeterministicSchedsParetoExplorer::Facet::Facet(storm::storage::geometry::Halfspace&& halfspace) : halfspace(std::move(halfspace)) { + // Intentionally left empty + } + + template + storm::storage::geometry::Halfspace const& DeterministicSchedsParetoExplorer::Facet::getHalfspace() const { + return halfspace; + } + + template + void DeterministicSchedsParetoExplorer::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); + } + + template + std::vector::PointId> const& DeterministicSchedsParetoExplorer::Facet::getPoints() const { + return paretoPointsOnFacet; + } + + template + uint64_t DeterministicSchedsParetoExplorer::Facet::getNumberOfPoints() const { + return paretoPointsOnFacet.size(); + } + + + template + typename DeterministicSchedsParetoExplorer::Polytope const& DeterministicSchedsParetoExplorer::Facet::getInducedSimplex(Pointset const& pointset, std::vector const& referenceCoordinates) { + if (!inducedSimplex) { + std::vector> vertices = {referenceCoordinates}; + for (auto const& pId : paretoPointsOnFacet) { + vertices.push_back(pointset.getPoint(pId).get()); + } + // This facet might lie at the 'border', which means that the downward closure has to be taken in some directions + storm::storage::BitVector dimensionsForDownwardClosure = storm::utility::vector::filterZero(this->halfspace.normalVector()); + STORM_LOG_ASSERT(dimensionsForDownwardClosure.getNumberOfSetBits() + vertices.size() >= halfspace.normalVector().size() + 1, "The number of points on the facet is insufficient"); + if (dimensionsForDownwardClosure.empty()) { + inducedSimplex = storm::storage::geometry::Polytope::create(vertices); + } else { + inducedSimplex = storm::storage::geometry::Polytope::createSelectiveDownwardClosure(vertices, dimensionsForDownwardClosure); + } + } + return inducedSimplex; + } + + template + DeterministicSchedsParetoExplorer::FacetAnalysisContext::FacetAnalysisContext(Facet& f) : facet(f) { + // Intentionally left empty + } + + template + DeterministicSchedsParetoExplorer::DeterministicSchedsParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) : model(preprocessorResult.preprocessedModel), objectives(preprocessorResult.objectives) { + originalModelInitialState = *preprocessorResult.originalModel.getInitialStates().begin(); + lpChecker = std::make_shared>(*model, objectives); + } + + template + std::unique_ptr DeterministicSchedsParetoExplorer::check(Environment const& env) { + + clean(); + initializeFacets(env); + while (!unprocessedFacets.empty()) { + Facet f = std::move(unprocessedFacets.front()); + unprocessedFacets.pop(); + processFacet(env, f); + } + + std::vector>paretoPoints; + paretoPoints.reserve(pointset.size()); + 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); + } + + template + void DeterministicSchedsParetoExplorer::clean() { + pointset = Pointset(); + unprocessedFacets = std::queue(); + overApproximation = storm::storage::geometry::Polytope::createUniversalPolytope(); + unachievableAreas.clear(); + } + + template + void DeterministicSchedsParetoExplorer::addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace) { + GeometryValueType offset = storm::utility::vector::dotProduct(normalVector, pointOnHalfspace.get()); + if (env.modelchecker().multi().isPrintResultsSet()) { + std::cout << "## unachievable halfspace: ["; + bool first = true; + for (auto const& xi : normalVector) { + if (first) { + first = false; + } else { + std::cout << ","; + } + std::cout << storm::utility::convertNumber(xi); + } + std::cout << "];[" << storm::utility::convertNumber(offset) << "]" << std::endl; + } + storm::storage::geometry::Halfspace overApproxHalfspace(normalVector, offset); + overApproximation = overApproximation->intersection(overApproxHalfspace); + } + + template + void DeterministicSchedsParetoExplorer::addUnachievableArea(Environment const& env, Polytope const& area) { + if (env.modelchecker().multi().isPrintResultsSet()) { + std::vector> vertices; + if (objectives.size() == 2) { + vertices = area->getVerticesInClockwiseOrder(); + } else { + vertices = area->getVertices(); + } + std::cout << "## unachievable polytope: "; + bool firstVertex = true; + for (auto const& v : vertices) { + if (firstVertex) { + firstVertex = false; + } else { + std::cout << ";"; + } + std::cout << "["; + bool firstEntry = true; + for (auto const& vi : v) { + if (firstEntry) { + firstEntry = false; + } else { + std::cout << ","; + } + std::cout << storm::utility::convertNumber(vi); + } + std::cout << "]"; + } + std::cout << std::endl; + } + unachievableAreas.push_back(area); + } + + template + void DeterministicSchedsParetoExplorer::initializeFacets(Environment const& env) { + 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(); + } else { + weightVector[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)); + } + } + + auto initialHalfspaces = pointset.downwardClosure()->getHalfspaces(); + for (auto& h : initialHalfspaces) { + Facet f(std::move(h)); + for (auto const& p : pointset) { + if (f.getHalfspace().isPointOnBoundary(p.second.get())) { + 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."); + if (!checkFacetPrecision(env, f)) { + unprocessedFacets.push(std::move(f)); + } + } + } + + template + 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(); + 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()); + + if (optimizeAndSplitFacet(env,f)) { + return; + } + + GeometryValueType eps = storm::utility::convertNumber(env.modelchecker().multi().getPrecision()); + eps += eps; // The unknown area (box) can actually have size 2*eps + storm::storage::geometry::PolytopeTree polytopeTree(f.getInducedSimplex(pointset, getReferenceCoordinates())); + for (auto const& point : pointset) { + polytopeTree.substractDownwardClosure(point.second.get(), eps); + if (polytopeTree.isEmpty()) { + break; + } + } + if (!polytopeTree.isEmpty()) { + auto res = lpChecker->check(env, polytopeTree, eps); + for (auto const& infeasableArea : res.second) { + addUnachievableArea(env, infeasableArea); + } + for (auto& achievablePoint : res.first) { + 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(); + } + } + 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) { + auto const& optPoint = pointset.getPoint(*optPointId); + // TODO: this check might suffer from numerical errors. Check how much this would hurt us. + if (f.getHalfspace().contains(optPoint.get())) { + // The point is contained in the halfspace which means that no more splitting is possible. + return false; + } else { + // Found a new Pareto optimal point -> generate new facets + std::vector> vertices; + vertices.push_back(optPoint.get()); + for (auto const& pId : f.getPoints()) { + vertices.push_back(pointset.getPoint(pId).get()); + } + auto newHalfspaceCandidates = storm::storage::geometry::Polytope::createSelectiveDownwardClosure(vertices, storm::utility::vector::filterZero(f.getHalfspace().normalVector()))->getHalfspaces(); + for (auto& h : newHalfspaceCandidates) { + 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(), 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, pointset.getPoint(pId)); + } + ++vertexIt; + } + assert(vertexIt == vertices.end()); + if (!checkFacetPrecision(env, fNew)) { + unprocessedFacets.push(std::move(fNew)); + } + } + } + return true; + } + } else { + // If the 'optimal point' was dominated by an existing point, we can not split the facet any further. + 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) { + /* + STORM_LOG_ERROR_COND(objectives.size()==2, "Exporting plot requested but this is only implemented for the two-dimensional case."); + + auto transformedUnderApprox = transformPolytopeToOriginalModel(underApproximation); + auto transformedOverApprox = transformPolytopeToOriginalModel(overApproximation); + + // Get pareto points as well as a hyperrectangle that is used to guarantee that the resulting polytopes are bounded. + storm::storage::geometry::Hyperrectangle boundaries(std::vector(objectives.size(), storm::utility::zero()), std::vector(objectives.size(), storm::utility::zero())); + std::vector> paretoPoints; + paretoPoints.reserve(refinementSteps.size()); + for(auto const& step : refinementSteps) { + paretoPoints.push_back(transformPointToOriginalModel(step.lowerBoundPoint)); + boundaries.enlarge(paretoPoints.back()); + } + auto underApproxVertices = transformedUnderApprox->getVertices(); + for(auto const& v : underApproxVertices) { + boundaries.enlarge(v); + } + auto overApproxVertices = transformedOverApprox->getVertices(); + for(auto const& v : overApproxVertices) { + boundaries.enlarge(v); + } + + //Further enlarge the boundaries a little + storm::utility::vector::scaleVectorInPlace(boundaries.lowerBounds(), GeometryValueType(15) / GeometryValueType(10)); + storm::utility::vector::scaleVectorInPlace(boundaries.upperBounds(), GeometryValueType(15) / GeometryValueType(10)); + + auto boundariesAsPolytope = boundaries.asPolytope(); + std::vector columnHeaders = {"x", "y"}; + + std::vector> pointsForPlotting; + if (env.modelchecker().multi().getPlotPathUnderApproximation()) { + underApproxVertices = transformedUnderApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder(); + pointsForPlotting.reserve(underApproxVertices.size()); + for(auto const& v : underApproxVertices) { + pointsForPlotting.push_back(storm::utility::vector::convertNumericVector(v)); + } + storm::utility::exportDataToCSVFile(env.modelchecker().multi().getPlotPathUnderApproximation().get(), pointsForPlotting, columnHeaders); + } + + if (env.modelchecker().multi().getPlotPathOverApproximation()) { + pointsForPlotting.clear(); + overApproxVertices = transformedOverApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder(); + pointsForPlotting.reserve(overApproxVertices.size()); + for(auto const& v : overApproxVertices) { + pointsForPlotting.push_back(storm::utility::vector::convertNumericVector(v)); + } + storm::utility::exportDataToCSVFile(env.modelchecker().multi().getPlotPathOverApproximation().get(), pointsForPlotting, columnHeaders); + } + + if (env.modelchecker().multi().getPlotPathParetoPoints()) { + pointsForPlotting.clear(); + pointsForPlotting.reserve(paretoPoints.size()); + for(auto const& v : paretoPoints) { + pointsForPlotting.push_back(storm::utility::vector::convertNumericVector(v)); + } + storm::utility::exportDataToCSVFile(env.modelchecker().multi().getPlotPathParetoPoints().get(), pointsForPlotting, columnHeaders); + } + }; + */ + } + + template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; + template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; + template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; + template class DeterministicSchedsParetoExplorer, storm::RationalNumber>; + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h new file mode 100644 index 000000000..bbf074b5c --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h @@ -0,0 +1,226 @@ +#pragma once + +#include +#include + +#include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h" + +#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 { + + class Environment; + + namespace modelchecker { + namespace multiobjective { + + template + class DeterministicSchedsParetoExplorer { + public: + typedef uint64_t PointId; + typedef typename std::shared_ptr> Polytope; + typedef typename SparseModelType::ValueType ModelValueType; + + class Point { + public: + Point(std::vector const& coordinates); + Point(std::vector&& coordinates); + + std::vector const& get() const; + std::vector& get(); + + uint64_t dimension() const; + + enum class DominanceResult { + Incomparable, + Dominates, + Dominated, + Equal + }; + DominanceResult getDominance(Point const& other) const; + + 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; + }; + + + class Pointset { + public: + + typedef typename std::map::const_iterator iterator_type; + + + Pointset(); + + /*! + * If the given point is not dominated by another point in the set, it is added + * to the set and its ID is returned. + * If the point is dominated by another point, boost::none is returned. + * Erases all points in the set, that are dominated by the given point. + * If the same point is already contained in the set, its id is returned + */ + boost::optional addPoint(Environment const& env, Point&& point); + + /*! + * Returns the point with the given ID + */ + Point const& getPoint(PointId const& id) const; + + iterator_type begin() const; + iterator_type end() const; + + /*! + * Returns the number of points currently contained in the set + */ + uint64_t size() const; + + /*! + * Returns the downward closure of the contained points. + */ + Polytope downwardClosure() const; + + void collectPointsInPolytope(std::set& collectedPoints, Polytope const& polytope); + + void printToStream(std::ostream& out, bool includeIDs = true, bool convertToDouble = false); + + private: + std::map points; + PointId currId; + }; + + class Facet { + public: + 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, Point const& point); + std::vector const& getPoints() const; + uint64_t getNumberOfPoints() const; + + /*! + * Creates a polytope that captures all points that lie 'under' the facet + */ + Polytope const& getInducedSimplex(Pointset const& pointset, std::vector const& referenceCoordinates); + + + + private: + storm::storage::geometry::Halfspace halfspace; + std::vector paretoPointsOnFacet; + 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; + }; + + + DeterministicSchedsParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult); + + virtual std::unique_ptr check(Environment const& env); + + void exportPlotOfCurrentApproximation(Environment const& env); + + private: + + /*! + * Cleans up all cached results from a previous check call + */ + void clean(); + + /*! + * Intersects the overapproximation with the given halfspace + */ + void addHalfspaceToOverApproximation(Environment const& env, std::vector const& normalVector, Point const& pointOnHalfspace); + + /*! + * Adds a polytope which consists of unachievable points + */ + void addUnachievableArea(Environment const& env, Polytope const& area); + + /*! + * Builds the initial facets by optimizing the objectives individually. + * Adds the facets that need further processing to unprocessedFacets + */ + void initializeFacets(Environment const& env); + + /*! + * Gets reference coordinates used to subdividing the downwardclosure + */ + 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 + */ + 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 + * 2. New facets are generated and (if not already precise enough) added to unprocessedFacets + * 3. true is returned + */ + 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); + + Pointset pointset; + std::queue unprocessedFacets; + Polytope overApproximation; + std::vector unachievableAreas; + + std::shared_ptr> lpChecker; + std::shared_ptr const& model; + uint64_t originalModelInitialState; + std::vector> const& objectives; + }; + + } + } +} \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp b/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp index 0c15d9a81..68e9fcd5e 100644 --- a/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp +++ b/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp @@ -10,7 +10,7 @@ #include "storm/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.h" #include "storm/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.h" #include "storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.h" -#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicParetoExplorer.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h" #include "storm/settings/SettingsManager.h" #include "storm/settings/modules/CoreSettings.h" #include "storm/utility/Stopwatch.h" @@ -53,7 +53,7 @@ namespace storm { { if (env.modelchecker().multi().isSchedulerRestrictionSet()) { STORM_LOG_THROW(preprocessorResult.queryType == preprocessing::SparseMultiObjectivePreprocessorResult::QueryType::Pareto, storm::exceptions::NotImplementedException, "Currently, only Pareto queries with scheduler restrictions are implemented."); - auto explorer = DeterministicParetoExplorer(preprocessorResult); + auto explorer = DeterministicSchedsParetoExplorer(preprocessorResult); result = explorer.check(env); if (env.modelchecker().multi().isExportPlotSet()) { explorer.exportPlotOfCurrentApproximation(env); diff --git a/src/storm/storage/geometry/PolytopeTree.h b/src/storm/storage/geometry/PolytopeTree.h new file mode 100644 index 000000000..bf5306752 --- /dev/null +++ b/src/storm/storage/geometry/PolytopeTree.h @@ -0,0 +1,142 @@ +#pragma once + +#include +#include +#include "storm/storage/geometry/Polytope.h" + +namespace storm { + namespace storage { + namespace geometry { + + /*! + * Represents a set of points in Euclidean space. + * The set is defined as the union of the polytopes at the leafs of the tree. + * The polytope at inner nodes should always be the convex union of its children. + * The sets described by the children of a node are disjoint. + * A child is always non-empty, i.e., isEmpty() should only hold for the root node. + */ + template + class PolytopeTree { + + public: + PolytopeTree(std::shared_ptr> const& polytope = nullptr) : polytope(polytope) { + // Intentionally left empty + } + + /*! + * Substracts the given rhs from this polytope. + */ + void setMinus(std::shared_ptr> const& rhs) { + // This operation only has an effect if the intersection of this and rhs is non-empty. + if (!isEmpty() && !polytope->intersection(rhs)->isEmpty()) { + if (children.empty()) { + // This is a leaf node. + // Apply splitting. + auto newChildren = polytope->setMinus(rhs); + if (newChildren.empty()) { + // Delete this node. + polytope = nullptr; + } else if (newChildren.size() == 1) { + // Replace this node with its only child + polytope = newChildren.front()->clean(); + } else { + // Add the new children to this node. There is no need to traverse them. + for (auto& c : newChildren) { + children.push_back(c->clean()); + } + } + } else { + // This is an inner node. Traverse the children and set this to the convex union of its children. + std::vector> newChildren; + std::shared_ptr> newPolytope = nullptr; + for (auto& c : children) { + c.setMinus(rhs); + if (c.polytope != nullptr) { + newChildren.push_back(c); + if (newPolytope) { + newPolytope->convexUnion(c.polytope); + } else { + newPolytope = c.polytope; + } + } + } + polytope = newPolytope; // nullptr, if no children left + children = std::move(newChildren); + } + } + } + + /*! + * Substracts the downward closure of the given point from this set. + * @param point the given point + * @param offset + */ + void substractDownwardClosure(std::vector const& point, ValueType const& offset = storm::utility::zero()) { + if (storm::utility::isZero(offset)) { + setMinus(Polytope::createDownwardClosure({point})); + } else { + std::vector pointPrime; + pointPrime.reserve(point.size()); + for (auto const& coordinate : point) { + pointPrime.push_back(coordinate + offset); + } + setMinus(Polytope::createDownwardClosure({pointPrime})); + } + } + + /*! + * Returns true if this is the empty set. + */ + bool isEmpty() const { + return polytope == nullptr; + } + + /*! + * Clears all contents of this set, making it the empty set. + */ + void clear() { + children.clear(); + polytope = nullptr; + } + + /*! + * Gets the polytope at this node + */ + std::shared_ptr>& getPolytope() { + return polytope; + } + + /*! + * Gets the children at this node. + */ + std::vector& getChildren() { + return children; + } + + /*! + * Returns a string representation of this node (for debugging purposes) + */ + std::string toString() { + if (isEmpty()) { + return "Empty PolytopeTree"; + } + std::stringstream s; + s << "PolytopeTree node with " << getChildren().size() << " children: " << getPolytope()->toString(true) << std::endl << "Vertices: "; + auto vertices = getPolytope()->getVertices(); + for (auto const& v : vertices) { + s << "["; + for (auto const& vi : v) { + s << storm::utility::convertNumber(vi) << ","; + } + s << "]\t"; + } + s << std::endl; + return s.str(); + } + private: + std::shared_ptr> polytope; + std::vector> children; + }; + } + } +} \ No newline at end of file