sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

deterministic schedulers: Refactored code for lp-based checker.

main
Tim Quatmann 7 years ago
parent
886bead969
commit
ee090b630e
8 changed files with 1596 additions and 2 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 253
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp
  2. 71
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h
  3. 159
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp
  4. 44
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h
  5. 699
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp
  6. 226
      src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h
  7. 4
      src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp
  8. 142
      src/storm/storage/geometry/PolytopeTree.h

253
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp
View File

@ -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 <typename ModelType, typename GeometryValueType>
DeterministicSchedsLpChecker<ModelType, GeometryValueType>::DeterministicSchedsLpChecker(ModelType const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) : model(model) {
swInit.start();
initializeObjectiveHelper(objectives);
initializeLpModel();
swInit.stop();
}
template <typename ModelType, typename GeometryValueType>
DeterministicSchedsLpChecker<ModelType, GeometryValueType>::~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 <typename ModelType, typename GeometryValueType>
void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::setCurrentWeightVector(std::vector<GeometryValueType> 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<ValueType>(weightVector[objIndex])));
lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == initialStateResults[objIndex]);
}
lpModel->update();
swLpBuild.stop();
}
template <typename ModelType, typename GeometryValueType>
std::vector<GeometryValueType> DeterministicSchedsLpChecker<ModelType, GeometryValueType>::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<GeometryValueType>(lpModel->getContinuousValue(objVar)));
}
swCheckVertices.stop(); swCheck.stop();
return newPoint;
}
template <typename ModelType, typename GeometryValueType>
std::pair<std::vector<std::vector<GeometryValueType>>, std::vector<std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>>>> DeterministicSchedsLpChecker<ModelType, GeometryValueType>::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& 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<Point> foundPoints;
std::vector<Polytope> infeasableAreas;
checkRecursive(polytopeTree, eps, foundPoints, infeasableAreas);
swCheck.stop();
std::cout << " done!" << std::endl;
return {foundPoints, infeasableAreas};
}
template <typename ModelType, typename GeometryValueType>
void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::initializeObjectiveHelper(std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives) {
objectiveHelper.reserve(objectives.size());
for (auto const& obj : objectives) {
objectiveHelper.emplace_back(model, obj);
}
}
template <typename ModelType, typename GeometryValueType>
void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::initializeLpModel() {
uint64_t numStates = model.getNumberOfStates();
lpModel = storm::utility::solver::getLpSolver<ValueType>("model");
lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize);
initialStateResults.clear();
auto one = lpModel->getConstant(storm::utility::one<ValueType>());
// Create choice variables and assert that at least one choice is taken at each state.
std::vector<storm::expressions::Expression> 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<storm::expressions::Expression> 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<storm::expressions::Expression> 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 <typename ModelType, typename GeometryValueType>
void DeterministicSchedsLpChecker<ModelType, GeometryValueType>::checkRecursive(storm::storage::geometry::PolytopeTree <GeometryValueType>& polytopeTree, GeometryValueType const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& 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<GeometryValueType>(lpModel->getContinuousValue(objVar)));
}
auto halfspace = storm::storage::geometry::Halfspace<GeometryValueType>(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<GeometryValueType>::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::models::sparse::Mdp<double>, storm::RationalNumber>;
template class DeterministicSchedsLpChecker<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
template class DeterministicSchedsLpChecker<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
template class DeterministicSchedsLpChecker<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
}
}
}

71
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h
View File

@ -0,0 +1,71 @@
#pragma once
#include <vector>
#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 <typename ModelType, typename GeometryValueType>
class DeterministicSchedsLpChecker {
public:
typedef typename ModelType::ValueType ValueType;
typedef typename std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> Polytope;
typedef typename std::vector<GeometryValueType> Point;
DeterministicSchedsLpChecker(ModelType const& model, std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
~DeterministicSchedsLpChecker();
/*!
* Specifies the current direction.
*/
void setCurrentWeightVector(std::vector<GeometryValueType> 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<Point>, std::vector<Polytope>> check(storm::Environment const& env, storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps);
private:
void initializeObjectiveHelper(std::vector<storm::modelchecker::multiobjective::Objective<ValueType>> const& objectives);
void initializeLpModel();
void checkRecursive(storm::storage::geometry::PolytopeTree<GeometryValueType>& polytopeTree, GeometryValueType const& eps, std::vector<Point>& foundPoints, std::vector<Polytope>& infeasableAreas);
ModelType const& model;
std::vector<DeterministicSchedsObjectiveHelper<ModelType>> objectiveHelper;
std::unique_ptr<storm::solver::LpSolver<ValueType>> lpModel;
std::vector<storm::expressions::Expression> initialStateResults;
std::vector<storm::expressions::Variable> currentObjectiveVariables;
std::vector<GeometryValueType> 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;
};
}
}
}

159
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp
View File

@ -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 <typename ModelType>
DeterministicSchedsObjectiveHelper<ModelType>::DeterministicSchedsObjectiveHelper(ModelType const& model, storm::modelchecker::multiobjective::Objective<ValueType> const& objective) : model(model), objective(objective) {
// Intentionally left empty
}
template <typename ModelType>
storm::storage::BitVector evaluatePropositionalFormula(ModelType const& model, storm::logic::Formula const& formula) {
storm::modelchecker::SparsePropositionalModelChecker<ModelType> 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 <typename ModelType>
std::map<uint64_t, typename ModelType::ValueType> const& DeterministicSchedsObjectiveHelper<ModelType>::getSchedulerIndependentStateValues() const {
if (!schedulerIndependentStateValues) {
auto const& formula = *objective.formula;
std::map<uint64_t, ValueType> 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<ValueType>();
}
}
{
storm::storage::BitVector prob0States = storm::utility::graph::performProb0A(backwardTransitions, phiStates, psiStates);
for (auto const& prob0State : prob0States) {
result[prob0State] = storm::utility::zero<ValueType>();
}
}
} 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<ValueType>();
}
} 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<ValueType>();
}
} else {
STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported.");
}
schedulerIndependentStateValues = std::move(result);
}
return schedulerIndependentStateValues.get();
}
template <typename ModelType>
std::map<uint64_t, typename ModelType::ValueType> const& DeterministicSchedsObjectiveHelper<ModelType>::getChoiceValueOffsets() const {
if (!choiceValueOffsets) {
auto const& formula = *objective.formula;
auto const& subformula = formula.getSubformula();
std::map<uint64_t, ValueType> 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<ValueType> 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<ValueType> const* rates = nullptr;
storm::storage::BitVector const* ms = nullptr;
if (model.isOfType(storm::models::ModelType::MarkovAutomaton)) {
auto ma = model.template as<storm::models::sparse::MarkovAutomaton<ValueType>>();
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<ValueType>();
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>
typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper<ModelType>::getUpperValueBoundAtState(uint64_t state) const{
return objective.upperResultBound.get();
}
template <typename ModelType>
typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper<ModelType>::getLowerValueBoundAtState(uint64_t state) const{
return objective.lowerResultBound.get();
}
template class DeterministicSchedsObjectiveHelper<storm::models::sparse::Mdp<double>>;
template class DeterministicSchedsObjectiveHelper<storm::models::sparse::Mdp<storm::RationalNumber>>;
template class DeterministicSchedsObjectiveHelper<storm::models::sparse::MarkovAutomaton<double>>;
template class DeterministicSchedsObjectiveHelper<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>>;
}
}
}

44
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h
View File

@ -0,0 +1,44 @@
#pragma once
#include <map>
#include <boost/optional.hpp>
#include "storm/modelchecker/multiobjective/Objective.h"
namespace storm {
namespace modelchecker {
namespace multiobjective {
template <typename ModelType>
class DeterministicSchedsObjectiveHelper {
public:
typedef typename ModelType::ValueType ValueType;
DeterministicSchedsObjectiveHelper(ModelType const& model, Objective<ValueType> const& objective);
/*!
* Returns states and values for states that are independent of the scheduler.
*/
std::map<uint64_t, ValueType> 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<uint64_t, ValueType> const& getChoiceValueOffsets() const;
ValueType const& getUpperValueBoundAtState(uint64_t state) const;
ValueType const& getLowerValueBoundAtState(uint64_t state) const;
private:
mutable boost::optional<std::map<uint64_t, ValueType>> schedulerIndependentStateValues;
mutable boost::optional<std::map<uint64_t, ValueType>> choiceValueOffsets;
ModelType const& model;
Objective<ValueType> const& objective;
};
}
}
}

699
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp
View File

@ -0,0 +1,699 @@
#include <sstream>
#include <algorithm>
#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 <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::Point(std::vector<GeometryValueType> const& coordinates) : coordinates(coordinates), paretoOptimal(false), onFacet(false) {
STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported");
}
template <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::Point(std::vector<GeometryValueType>&& coordinates) : coordinates(std::move(coordinates)), paretoOptimal(false), onFacet(false) {
STORM_LOG_ASSERT(!this->coordinates.empty(), "Points with dimension 0 are not supported");
}
template <class SparseModelType, typename GeometryValueType>
std::vector<GeometryValueType>& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::get() {
return coordinates;
}
template <class SparseModelType, typename GeometryValueType>
std::vector<GeometryValueType> const& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::get() const {
return coordinates;
}
template <class SparseModelType, typename GeometryValueType>
uint64_t DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::dimension() const {
STORM_LOG_ASSERT(!coordinates.empty(), "Points with dimension 0 are not supported");
return coordinates.size();
}
template <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::DominanceResult DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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 <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::setParetoOptimal(bool value) {
paretoOptimal = value;
}
template <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::isParetoOptimal() const {
return paretoOptimal;
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::setOnFacet(bool value) {
onFacet = value;
}
template <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point::liesOnFacet() const {
return onFacet;
}
template <class SparseModelType, typename GeometryValueType>
std::string DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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<double>(pi);
} else {
out << pi;
}
}
return out.str();
}
// template <class SparseModelType, typename GeometryValueType>
// bool operator<(typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point const& lhs, typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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 <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::Pointset() : currId(1) {
// Intentionally left empty
}
template <class SparseModelType, typename GeometryValueType>
boost::optional<typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::PointId> DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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<double>(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 <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Point const& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::getPoint(PointId const& id) const {
return points.at(id);
}
template <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::iterator_type DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::begin() const {
return points.begin();
}
template <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::iterator_type DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::end() const {
return points.end();
}
template <class SparseModelType, typename GeometryValueType>
uint64_t DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::size() const {
return points.size();
}
template <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Polytope DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::downwardClosure() const {
std::vector<std::vector<GeometryValueType>> pointsAsVector;
pointsAsVector.reserve(size());
for (auto const& p : points) {
pointsAsVector.push_back(p.second.get());
}
return storm::storage::geometry::Polytope<GeometryValueType>::createDownwardClosure(std::move(pointsAsVector));
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Pointset::collectPointsInPolytope(std::set<PointId>& collectedPoints, Polytope const& polytope) {
for (auto const& p : points) {
if (polytope->contains(p.second.get())) {
collectedPoints.insert(p.first);
}
}
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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 <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::Facet(storm::storage::geometry::Halfspace<GeometryValueType> const& halfspace) : halfspace(halfspace) {
// Intentionally left empty
}
template <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::Facet(storm::storage::geometry::Halfspace<GeometryValueType>&& halfspace) : halfspace(std::move(halfspace)) {
// Intentionally left empty
}
template <class SparseModelType, typename GeometryValueType>
storm::storage::geometry::Halfspace<GeometryValueType> const& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::getHalfspace() const {
return halfspace;
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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<double>(getHalfspace().euclideanDistance(point.get())));
} else {
STORM_LOG_DEBUG("Halfspace of facet is shifted by " << storm::utility::convertNumber<double>(getHalfspace().euclideanDistance(point.get())) << " to capture all points that are supposed to lie on the facet.");
halfspace.offset() = product;
}
}
paretoPointsOnFacet.push_back(pointId);
}
template <class SparseModelType, typename GeometryValueType>
std::vector<typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::PointId> const& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::getPoints() const {
return paretoPointsOnFacet;
}
template <class SparseModelType, typename GeometryValueType>
uint64_t DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::getNumberOfPoints() const {
return paretoPointsOnFacet.size();
}
template <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Polytope const& DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::Facet::getInducedSimplex(Pointset const& pointset, std::vector<GeometryValueType> const& referenceCoordinates) {
if (!inducedSimplex) {
std::vector<std::vector<GeometryValueType>> 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<GeometryValueType>::create(vertices);
} else {
inducedSimplex = storm::storage::geometry::Polytope<GeometryValueType>::createSelectiveDownwardClosure(vertices, dimensionsForDownwardClosure);
}
}
return inducedSimplex;
}
template <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::FacetAnalysisContext::FacetAnalysisContext(Facet& f) : facet(f) {
// Intentionally left empty
}
template <class SparseModelType, typename GeometryValueType>
DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::DeterministicSchedsParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult<SparseModelType>& preprocessorResult) : model(preprocessorResult.preprocessedModel), objectives(preprocessorResult.objectives) {
originalModelInitialState = *preprocessorResult.originalModel.getInitialStates().begin();
lpChecker = std::make_shared<DeterministicSchedsLpChecker<SparseModelType, GeometryValueType>>(*model, objectives);
}
template <class SparseModelType, typename GeometryValueType>
std::unique_ptr<CheckResult> DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::check(Environment const& env) {
clean();
initializeFacets(env);
while (!unprocessedFacets.empty()) {
Facet f = std::move(unprocessedFacets.front());
unprocessedFacets.pop();
processFacet(env, f);
}
std::vector<std::vector<ModelValueType>>paretoPoints;
paretoPoints.reserve(pointset.size());
for (auto const& p : pointset) {
paretoPoints.push_back(storm::utility::vector::convertNumericVector<ModelValueType>(transformObjectiveValuesToOriginal(objectives, p.second.get())));
}
return std::make_unique<storm::modelchecker::ExplicitParetoCurveCheckResult<ModelValueType>>(originalModelInitialState, std::move(paretoPoints),
nullptr, nullptr);
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::clean() {
pointset = Pointset();
unprocessedFacets = std::queue<Facet>();
overApproximation = storm::storage::geometry::Polytope<GeometryValueType>::createUniversalPolytope();
unachievableAreas.clear();
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::addHalfspaceToOverApproximation(Environment const& env, std::vector<GeometryValueType> 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<double>(xi);
}
std::cout << "];[" << storm::utility::convertNumber<double>(offset) << "]" << std::endl;
}
storm::storage::geometry::Halfspace<GeometryValueType> overApproxHalfspace(normalVector, offset);
overApproximation = overApproximation->intersection(overApproxHalfspace);
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::addUnachievableArea(Environment const& env, Polytope const& area) {
if (env.modelchecker().multi().isPrintResultsSet()) {
std::vector<std::vector<GeometryValueType>> 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<double>(vi);
}
std::cout << "]";
}
std::cout << std::endl;
}
unachievableAreas.push_back(area);
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::initializeFacets(Environment const& env) {
for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
std::vector<GeometryValueType> weightVector(objectives.size(), storm::utility::zero<ModelValueType>());
if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
weightVector[objIndex] = -storm::utility::one<GeometryValueType>();
} else {
weightVector[objIndex] = storm::utility::one<GeometryValueType>();
}
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<ModelValueType>();
}
Point p(storm::utility::vector::convertNumericVector<GeometryValueType>(point));
p.setOnFacet();
// Adapt the overapproximation
std::vector<GeometryValueType> normalVector(objectives.size(), storm::utility::zero<GeometryValueType>());
normalVector[objIndex] = storm::utility::one<GeometryValueType>();
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<GeometryValueType>()) + f.getNumberOfPoints() == objectives.size(), "Unexpected number of points on facet.");
if (!checkFacetPrecision(env, f)) {
unprocessedFacets.push(std::move(f));
}
}
}
template <class SparseModelType, typename GeometryValueType>
std::vector<GeometryValueType> DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::getReferenceCoordinates() const {
std::vector<GeometryValueType> 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<GeometryValueType>(value));
}
return result;
}
template <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::checkFacetPrecision(Environment const& env, Facet& f) {
// TODO:
return false;
/*
auto const& inducedSimplex = f.getInducedSimplex(pointset);
GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
// get a polytope that contains exactly the points y, such that y+eps is in the induced simplex
std::vector<GeometryValueType> 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 <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::checkFacetPrecision(Environment const& env, Facet& f, std::set<PointId> const& collectedSimplexPoints) {
assert(false);
return false;
}
template <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::processFacet(Environment const& env, Facet& f) {
lpChecker->setCurrentWeightVector(f.getHalfspace().normalVector());
if (optimizeAndSplitFacet(env,f)) {
return;
}
GeometryValueType eps = storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
eps += eps; // The unknown area (box) can actually have size 2*eps
storm::storage::geometry::PolytopeTree<GeometryValueType> 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 <class SparseModelType, typename GeometryValueType>
typename DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::FacetAnalysisContext DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::createAnalysisContext(Environment const& env, Facet& f) {
FacetAnalysisContext res(f);
/*
res.expressionManager = std::make_shared<storm::expressions::ExpressionManager>();
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 <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::optimizeAndSplitFacet(Environment const& env, Facet& f) {
// Obtain the correct weight vector
auto weightVector = storm::utility::vector::convertNumericVector<ModelValueType>(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<ModelValueType>();
}
}
// Invoke optimization and insert the explored points
boost::optional<PointId> 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<ModelValueType>();
}
}
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<std::vector<GeometryValueType>> 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<GeometryValueType>::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 <class SparseModelType, typename GeometryValueType>
bool DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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 <class SparseModelType, typename GeometryValueType>
void DeterministicSchedsParetoExplorer<SparseModelType, GeometryValueType>::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<GeometryValueType> boundaries(std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()), std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()));
std::vector<std::vector<GeometryValueType>> 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<std::string> columnHeaders = {"x", "y"};
std::vector<std::vector<double>> 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<double>(v));
}
storm::utility::exportDataToCSVFile<double, std::string>(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<double>(v));
}
storm::utility::exportDataToCSVFile<double, std::string>(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<double>(v));
}
storm::utility::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathParetoPoints().get(), pointsForPlotting, columnHeaders);
}
};
*/
}
template class DeterministicSchedsParetoExplorer<storm::models::sparse::Mdp<double>, storm::RationalNumber>;
template class DeterministicSchedsParetoExplorer<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
template class DeterministicSchedsParetoExplorer<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
template class DeterministicSchedsParetoExplorer<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
}
}
}

226
src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h
View File

@ -0,0 +1,226 @@
#pragma once
#include <memory>
#include <queue>
#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 SparseModelType, typename GeometryValueType>
class DeterministicSchedsParetoExplorer {
public:
typedef uint64_t PointId;
typedef typename std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> Polytope;
typedef typename SparseModelType::ValueType ModelValueType;
class Point {
public:
Point(std::vector<GeometryValueType> const& coordinates);
Point(std::vector<GeometryValueType>&& coordinates);
std::vector<GeometryValueType> const& get() const;
std::vector<GeometryValueType>& 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<GeometryValueType> coordinates;
bool paretoOptimal;
bool onFacet;
};
class Pointset {
public:
typedef typename std::map<PointId, Point>::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<PointId> 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<PointId>& collectedPoints, Polytope const& polytope);
void printToStream(std::ostream& out, bool includeIDs = true, bool convertToDouble = false);
private:
std::map<PointId, Point> points;
PointId currId;
};
class Facet {
public:
Facet(storm::storage::geometry::Halfspace<GeometryValueType> const& halfspace);
Facet(storm::storage::geometry::Halfspace<GeometryValueType>&& halfspace);
storm::storage::geometry::Halfspace<GeometryValueType> const& getHalfspace() const;
void addPoint(PointId const& pointId, Point const& point);
std::vector<PointId> 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<GeometryValueType> const& referenceCoordinates);
private:
storm::storage::geometry::Halfspace<GeometryValueType> halfspace;
std::vector<PointId> paretoPointsOnFacet;
Polytope inducedSimplex;
};
struct FacetAnalysisContext {
FacetAnalysisContext(Facet& f);
Facet& facet;
std::set<PointId> collectedPoints;
std::unique_ptr<storm::solver::SmtSolver> smtSolver;
std::shared_ptr<storm::expressions::ExpressionManager> 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<storm::expressions::Variable> x, xMinusEps;
};
DeterministicSchedsParetoExplorer(preprocessing::SparseMultiObjectivePreprocessorResult<SparseModelType>& preprocessorResult);
virtual std::unique_ptr<CheckResult> 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<GeometryValueType> 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<GeometryValueType> 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<PointId> 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<Facet> unprocessedFacets;
Polytope overApproximation;
std::vector<Polytope> unachievableAreas;
std::shared_ptr<DeterministicSchedsLpChecker<SparseModelType, GeometryValueType>> lpChecker;
std::shared_ptr<SparseModelType> const& model;
uint64_t originalModelInitialState;
std::vector<Objective<ModelValueType>> const& objectives;
};
}
}
}

4
src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp
View File

@ -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<SparseModelType>::QueryType::Pareto, storm::exceptions::NotImplementedException, "Currently, only Pareto queries with scheduler restrictions are implemented.");
auto explorer = DeterministicParetoExplorer<SparseModelType, storm::RationalNumber>(preprocessorResult);
auto explorer = DeterministicSchedsParetoExplorer<SparseModelType, storm::RationalNumber>(preprocessorResult);
result = explorer.check(env);
if (env.modelchecker().multi().isExportPlotSet()) {
explorer.exportPlotOfCurrentApproximation(env);

142
src/storm/storage/geometry/PolytopeTree.h
View File

@ -0,0 +1,142 @@
#pragma once
#include <vector>
#include <memory>
#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 <typename ValueType>
class PolytopeTree {
public:
PolytopeTree(std::shared_ptr<Polytope<ValueType>> const& polytope = nullptr) : polytope(polytope) {
// Intentionally left empty
}
/*!
* Substracts the given rhs from this polytope.
*/
void setMinus(std::shared_ptr<Polytope<ValueType>> 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<PolytopeTree<ValueType>> newChildren;
std::shared_ptr<Polytope<ValueType>> 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<ValueType> const& point, ValueType const& offset = storm::utility::zero<ValueType>()) {
if (storm::utility::isZero(offset)) {
setMinus(Polytope<ValueType>::createDownwardClosure({point}));
} else {
std::vector<ValueType> pointPrime;
pointPrime.reserve(point.size());
for (auto const& coordinate : point) {
pointPrime.push_back(coordinate + offset);
}
setMinus(Polytope<ValueType>::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<Polytope<ValueType>>& getPolytope() {
return polytope;
}
/*!
* Gets the children at this node.
*/
std::vector<PolytopeTree>& 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<double>(vi) << ",";
}
s << "]\t";
}
s << std::endl;
return s.str();
}
private:
std::shared_ptr<Polytope<ValueType>> polytope;
std::vector<PolytopeTree<ValueType>> children;
};
}
}
}
Write Preview
Loading…
Cancel
Save