From fdace26bee6794a802658bb30969d1a283c868f3 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Thu, 14 Dec 2023 14:43:54 +0100 Subject: [PATCH] Revised LP Encoding for Multi-objective verification under simple strategies (#468) * Update findGUROBI.cmake to find Gurobi 11 * Revised implementation of multi-objective Achievability under simple strategies * Added support for achievability queries --- resources/cmake/find_modules/FindGUROBI.cmake | 5 +- .../examples/testfiles/mdp/multiobj_infrew.nm | 19 + .../MultiObjectiveModelCheckerEnvironment.cpp | 28 + .../MultiObjectiveModelCheckerEnvironment.h | 12 + ...eterministicSchedsAchievabilityChecker.cpp | 90 ++ .../DeterministicSchedsAchievabilityChecker.h | 46 + .../DeterministicSchedsLpChecker.cpp | 946 ++++++++------- .../DeterministicSchedsLpChecker.h | 18 +- .../DeterministicSchedsObjectiveHelper.cpp | 1012 ++++++++++------- .../DeterministicSchedsObjectiveHelper.h | 94 +- .../DeterministicSchedsParetoExplorer.cpp | 106 +- .../DeterministicSchedsParetoExplorer.h | 4 + .../VisitingTimesHelper.cpp | 151 +++ .../deterministicScheds/VisitingTimesHelper.h | 48 + .../multiObjectiveModelChecking.cpp | 19 +- .../modules/MultiObjectiveSettings.cpp | 47 +- .../settings/modules/MultiObjectiveSettings.h | 25 + ...ultiObjectiveSchedRestModelCheckerTest.cpp | 810 ++++++------- 18 files changed, 2007 insertions(+), 1473 deletions(-) create mode 100644 resources/examples/testfiles/mdp/multiobj_infrew.nm create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.cpp create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.cpp create mode 100644 src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h diff --git a/resources/cmake/find_modules/FindGUROBI.cmake b/resources/cmake/find_modules/FindGUROBI.cmake index e97101030c..b1ba090fa7 100644 --- a/resources/cmake/find_modules/FindGUROBI.cmake +++ b/resources/cmake/find_modules/FindGUROBI.cmake @@ -18,6 +18,7 @@ else (GUROBI_INCLUDE_DIR) find_path(GUROBI_INCLUDE_DIR NAMES gurobi_c++.h PATHS "$ENV{GUROBI_HOME}/include" + "/Library/gurobi1100/macos_universal2/include" "/Library/gurobi1000/macos_universal2/include" "/Library/gurobi951/macos_universal2/include" "/Library/gurobi950/mac64/include" @@ -40,6 +41,7 @@ find_path(GUROBI_INCLUDE_DIR find_library( GUROBI_LIBRARY NAMES gurobi + gurobi110 gurobi100 gurobi95 gurobi91 @@ -58,6 +60,7 @@ find_library( GUROBI_LIBRARY gurobi46 gurobi45 PATHS "$ENV{GUROBI_HOME}/lib" + "/Library/gurobi1100/macos_universal2/lib" "/Library/gurobi1000/macos_universal2/lib" "/Library/gurobi951/macos_universal2/lib" "/Library/gurobi950/mac64/lib" @@ -81,7 +84,7 @@ find_library( GUROBI_LIBRARY find_library( GUROBI_CXX_LIBRARY NAMES gurobi_c++ PATHS "$ENV{GUROBI_HOME}/lib" - "/Library/gurobi1000/macos_universal2/lib" + "/Library/gurobi1100/macos_universal2/lib" "/Library/gurobi951/macos_universal2/lib" "/Library/gurobi950/mac64/lib" "/Library/gurobi912/mac64/lib" diff --git a/resources/examples/testfiles/mdp/multiobj_infrew.nm b/resources/examples/testfiles/mdp/multiobj_infrew.nm new file mode 100644 index 0000000000..4c4dc3c7c0 --- /dev/null +++ b/resources/examples/testfiles/mdp/multiobj_infrew.nm @@ -0,0 +1,19 @@ +mdp +module main + + x : [0..2]; + [a] x=0 -> (x'=1); + [b] x=0 -> (x'=2); + [c] x=1 -> (x'=1); + [d] x=1 -> (x'=2); + [e] x=2 -> (x'=1); +endmodule + +rewards "one" + [b] true : 10; + [c] true : 1; +endrewards + +rewards "two" + [e] true : 1; +endrewards \ No newline at end of file diff --git a/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.cpp b/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.cpp index 908e7ad8cc..dd8ff75333 100644 --- a/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.cpp +++ b/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.cpp @@ -33,6 +33,12 @@ MultiObjectiveModelCheckerEnvironment::MultiObjectiveModelCheckerEnvironment() { } else if (multiobjectiveSettings.isFlowEncodingSet()) { encodingType = EncodingType::Flow; } + STORM_LOG_ASSERT(multiobjectiveSettings.isBsccDetectionViaOrderConstraintsSet() || multiobjectiveSettings.isBsccDetectionViaFlowConstraintsSet(), + "unexpected settings"); + bsccOrderEncoding = multiobjectiveSettings.isBsccDetectionViaOrderConstraintsSet(); + STORM_LOG_ASSERT(multiobjectiveSettings.isIndicatorConstraintsSet() || multiobjectiveSettings.isBigMConstraintsSet(), "unexpected settings"); + indicatorConstraints = multiobjectiveSettings.isIndicatorConstraintsSet(); + redundantBsccConstraints = multiobjectiveSettings.isRedundantBsccConstraintsSet(); if (multiobjectiveSettings.isMaxStepsSet()) { maxSteps = multiobjectiveSettings.getMaxSteps(); @@ -121,6 +127,28 @@ void MultiObjectiveModelCheckerEnvironment::setEncodingType(EncodingType const& encodingType = value; } +bool MultiObjectiveModelCheckerEnvironment::getUseIndicatorConstraints() const { + return indicatorConstraints; +} +void MultiObjectiveModelCheckerEnvironment::setUseIndicatorConstraints(bool value) { + indicatorConstraints = value; +} + +bool MultiObjectiveModelCheckerEnvironment::getUseBsccOrderEncoding() const { + return bsccOrderEncoding; +} +void MultiObjectiveModelCheckerEnvironment::setUseBsccOrderEncoding(bool value) { + bsccOrderEncoding = value; +} + +bool MultiObjectiveModelCheckerEnvironment::getUseRedundantBsccConstraints() const { + return redundantBsccConstraints; +} + +void MultiObjectiveModelCheckerEnvironment::setUseRedundantBsccConstraints(bool value) { + redundantBsccConstraints = value; +} + bool MultiObjectiveModelCheckerEnvironment::isMaxStepsSet() const { return maxSteps.is_initialized(); } diff --git a/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h b/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h index 3f697a5ed3..9a87aac958 100644 --- a/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h +++ b/src/storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h @@ -47,6 +47,15 @@ class MultiObjectiveModelCheckerEnvironment { EncodingType const& getEncodingType() const; void setEncodingType(EncodingType const& value); + bool getUseIndicatorConstraints() const; + void setUseIndicatorConstraints(bool value); + + bool getUseBsccOrderEncoding() const; + void setUseBsccOrderEncoding(bool value); + + bool getUseRedundantBsccConstraints() const; + void setUseRedundantBsccConstraints(bool value); + bool isMaxStepsSet() const; uint64_t const& getMaxSteps() const; void setMaxSteps(uint64_t const& value); @@ -69,6 +78,9 @@ class MultiObjectiveModelCheckerEnvironment { storm::RationalNumber precision; PrecisionType precisionType; EncodingType encodingType; + bool indicatorConstraints; + bool bsccOrderEncoding; + bool redundantBsccConstraints; boost::optional maxSteps; boost::optional schedulerRestriction; bool printResults; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.cpp new file mode 100644 index 0000000000..cedca9edfe --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.cpp @@ -0,0 +1,90 @@ +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h" + +#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h" +#include "storm/modelchecker/multiobjective/preprocessing/SparseMultiObjectivePreprocessorResult.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/StandardRewardModel.h" + +namespace storm::modelchecker::multiobjective { + +template +DeterministicSchedsAchievabilityChecker::DeterministicSchedsAchievabilityChecker( + preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult) + : model(preprocessorResult.preprocessedModel), originalModelInitialState(*preprocessorResult.originalModel.getInitialStates().begin()) { + for (auto const& obj : preprocessorResult.objectives) { + objectiveHelper.emplace_back(*model, obj); + if (!objectiveHelper.back().hasThreshold()) { + STORM_LOG_ASSERT(!optimizingObjectiveIndex.has_value(), + "Unexpected input: got a (quantitative) achievability query but there are more than one objectives without a threshold."); + optimizingObjectiveIndex = objectiveHelper.size() - 1; + } + } + lpChecker = std::make_shared>(*model, objectiveHelper); +} + +template +std::unique_ptr DeterministicSchedsAchievabilityChecker::check(Environment const& env) { + using Point = typename std::vector; + uint64_t const numObj = objectiveHelper.size(); + // Create a polytope that contains all the points that would certify achievability of the point induced by the objective threshold + std::vector> thresholdHalfspaces; + std::vector objVector; + for (uint64_t objIndex = 0; objIndex < numObj; ++objIndex) { + auto& obj = objectiveHelper[objIndex]; + obj.computeLowerUpperBounds(env); + if (!optimizingObjectiveIndex.has_value() || *optimizingObjectiveIndex != objIndex) { + objVector.assign(numObj, storm::utility::zero()); + objVector[objIndex] = -storm::utility::one(); + thresholdHalfspaces.emplace_back(objVector, storm::utility::convertNumber(-obj.getThreshold())); + } + } + auto thresholdPolytope = storm::storage::geometry::Polytope::create(std::move(thresholdHalfspaces)); + + // Set objective weights (none if this is a qualitative achievability query) + objVector.assign(numObj, storm::utility::zero()); + auto eps = objVector; + if (optimizingObjectiveIndex.has_value()) { + objVector[*optimizingObjectiveIndex] = storm::utility::one(); + eps[*optimizingObjectiveIndex] = env.modelchecker().multi().getPrecision() * storm::utility::convertNumber(2u); + if (env.modelchecker().multi().getPrecisionType() == MultiObjectiveModelCheckerEnvironment::PrecisionType::RelativeToDiff) { + eps[*optimizingObjectiveIndex] *= storm::utility::convertNumber( + objectiveHelper[*optimizingObjectiveIndex].getUpperValueBoundAtState(originalModelInitialState) - + objectiveHelper[*optimizingObjectiveIndex].getLowerValueBoundAtState(originalModelInitialState)); + } + STORM_LOG_INFO("Computing quantitative achievability value with precision " << eps[*optimizingObjectiveIndex] << "."); + } + lpChecker->setCurrentWeightVector(env, objVector); + + // Search achieving point + auto achievingPoint = lpChecker->check(env, thresholdPolytope, eps); + bool const isAchievable = achievingPoint.has_value(); + if (isAchievable) { + STORM_LOG_INFO( + "Found achievable point: " << storm::utility::vector::toString(achievingPoint->first) << " ( approx. " + << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(achievingPoint->first)) + << " )."); + if (optimizingObjectiveIndex.has_value()) { + using ValueType = typename SparseModelType::ValueType; + // Average between obtained lower- and upper bounds + auto result = + storm::utility::convertNumber(achievingPoint->first[*optimizingObjectiveIndex] + achievingPoint->second); + result /= storm::utility::convertNumber(2u); + if (objectiveHelper[*optimizingObjectiveIndex].minimizing()) { + result *= -storm::utility::one(); + } + return std::make_unique>(originalModelInitialState, result); + } + } + return std::make_unique(originalModelInitialState, isAchievable); +} + +template class DeterministicSchedsAchievabilityChecker, storm::RationalNumber>; +template class DeterministicSchedsAchievabilityChecker, storm::RationalNumber>; +template class DeterministicSchedsAchievabilityChecker, storm::RationalNumber>; +template class DeterministicSchedsAchievabilityChecker, storm::RationalNumber>; +} // namespace storm::modelchecker::multiobjective diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h new file mode 100644 index 0000000000..54026a21c7 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include + +namespace storm { +class Environment; + +namespace modelchecker { + +class CheckResult; + +namespace multiobjective { +template +class DeterministicSchedsLpChecker; + +template +class DeterministicSchedsObjectiveHelper; + +namespace preprocessing { +template +class SparseMultiObjectivePreprocessorResult; +} + +template +class DeterministicSchedsAchievabilityChecker { + public: + typedef typename SparseModelType::ValueType ModelValueType; + + DeterministicSchedsAchievabilityChecker(preprocessing::SparseMultiObjectivePreprocessorResult& preprocessorResult); + + virtual std::unique_ptr check(Environment const& env); + + private: + std::shared_ptr> lpChecker; + std::vector> objectiveHelper; + + std::shared_ptr model; + uint64_t const originalModelInitialState; + std::optional optimizingObjectiveIndex; +}; + +} // namespace multiobjective +} // namespace modelchecker +} // namespace storm \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index d2a9eea316..7a3a3ddf13 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -1,27 +1,21 @@ #include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h" +#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h" #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" - #include "storm/models/sparse/MarkovAutomaton.h" #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/StandardRewardModel.h" #include "storm/settings/SettingsManager.h" #include "storm/settings/modules/CoreSettings.h" #include "storm/storage/MaximalEndComponentDecomposition.h" -#include "storm/storage/Scheduler.h" #include "storm/storage/SparseMatrix.h" -#include "storm/utility/graph.h" #include "storm/utility/solver.h" -#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" - -#include -#include #include "storm/exceptions/InvalidOperationException.h" +#include "storm/exceptions/UnexpectedException.h" -namespace storm { -namespace modelchecker { -namespace multiobjective { +namespace storm::modelchecker::multiobjective { template DeterministicSchedsLpChecker::DeterministicSchedsLpChecker( @@ -71,19 +65,15 @@ void DeterministicSchedsLpChecker::setCurrentWeigh for (uint64_t objIndex = 0; objIndex < initialStateResults.size(); ++objIndex) { currentObjectiveVariables.push_back( lpModel->addUnboundedContinuousVariable("w_" + std::to_string(objIndex), storm::utility::convertNumber(weightVector[objIndex]))); - if (objectiveHelper[objIndex].minimizing() && flowEncoding) { - lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == -initialStateResults[objIndex]); - } else { - lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == initialStateResults[objIndex]); - } + lpModel->addConstraint("", currentObjectiveVariables.back().getExpression() == initialStateResults[objIndex]); } lpModel->update(); swAll.stop(); } template -boost::optional> DeterministicSchedsLpChecker::check(storm::Environment const& env, - Polytope overapproximation) { +std::optional, GeometryValueType>> DeterministicSchedsLpChecker::check( + storm::Environment const& env, Polytope overapproximation, Point const& eps) { swAll.start(); initialize(env); STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); @@ -93,19 +83,31 @@ boost::optional> DeterministicSchedsLpCheckeraddConstraint("", c); } + + if (!eps.empty()) { + STORM_LOG_ASSERT(currentWeightVector.size() == eps.size(), "Eps vector has unexpected size."); + // Specify the allowed gap between the obtained lower/upper objective bounds. + GeometryValueType milpGap = storm::utility::vector::dotProduct(currentWeightVector, eps); + lpModel->setMaximalMILPGap(storm::utility::convertNumber(milpGap), false); + } lpModel->update(); swCheckWeightVectors.start(); lpModel->optimize(); swCheckWeightVectors.stop(); ++numLpQueries; - // STORM_PRINT_AND_LOG("Writing model to file '" << std::to_string(numLpQueries) << ".lp'\n";); - // lpModel->writeModelToFile(std::to_string(numLpQueries) + ".lp"); - boost::optional result; + // STORM_PRINT_AND_LOG("Writing model to file '" << std::to_string(numLpQueries) << ".lp'\n";); + // lpModel->writeModelToFile(std::to_string(numLpQueries) + ".lp"); + std::optional> result; if (!lpModel->isInfeasible()) { STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); swValidate.start(); - result = validateCurrentModel(env); + auto resultPoint = validateCurrentModel(env); swValidate.stop(); + auto resultValue = storm::utility::vector::dotProduct(resultPoint, currentWeightVector); + if (!eps.empty()) { + resultValue += storm::utility::convertNumber(lpModel->getMILPGap(false)); + } + result = std::make_pair(resultPoint, resultValue); } lpModel->pop(); STORM_LOG_TRACE("\t Done checking a vertex..."); @@ -140,233 +142,445 @@ DeterministicSchedsLpChecker::check(storm::Environ } template -std::map getSubEndComponents(storm::storage::SparseMatrix const& mecTransitions) { - auto backwardTransitions = mecTransitions.transpose(true); - std::map unprocessed, processed; - storm::storage::BitVector allStates(mecTransitions.getRowGroupCount(), true); - storm::storage::BitVector allChoices(mecTransitions.getRowCount(), true); - unprocessed[allStates] = allChoices; - while (!unprocessed.empty()) { - auto currentIt = unprocessed.begin(); - storm::storage::BitVector currentStates = currentIt->first; - storm::storage::BitVector currentChoices = currentIt->second; - unprocessed.erase(currentIt); - - bool hasSubEc = false; - for (auto removedState : currentStates) { - storm::storage::BitVector subset = currentStates; - subset.set(removedState, false); - storm::storage::MaximalEndComponentDecomposition subMecs(mecTransitions, backwardTransitions, subset); - for (auto const& subMec : subMecs) { - hasSubEc = true; - // Convert to bitvector - storm::storage::BitVector newEcStates(currentStates.size(), false), newEcChoices(currentChoices.size(), false); - for (auto const& stateChoices : subMec) { - newEcStates.set(stateChoices.first, true); - for (auto const& choice : stateChoices.second) { - newEcChoices.set(choice, true); - } +auto createChoiceVariables(storm::solver::LpSolver& lpModel, storm::storage::SparseMatrix const& matrix) { + std::vector choiceVariables; + choiceVariables.reserve(matrix.getRowCount()); + for (uint64_t state = 0; state < matrix.getRowGroupCount(); ++state) { + auto choices = matrix.getRowGroupIndices(state); + if (choices.size() == 1) { + choiceVariables.push_back(lpModel.getConstant(storm::utility::one())); // Unique choice; no variable necessary + } else { + std::vector localChoices; + for (auto const choice : choices) { + localChoices.push_back(lpModel.addBinaryVariable("a" + std::to_string(choice))); + choiceVariables.push_back(localChoices.back()); + } + lpModel.update(); + lpModel.addConstraint("", storm::expressions::sum(localChoices) == lpModel.getConstant(1)); + } + } + return choiceVariables; +} + +template +std::vector classicConstraints(storm::solver::LpSolver& lpModel, bool const& indicatorConstraints, + storm::storage::SparseMatrix const& matrix, uint64_t initialState, uint64_t objIndex, + HelperType const& objectiveHelper, + std::vector const& choiceVariables) { + // Create variables + std::vector objectiveValueVariables(matrix.getRowGroupCount()); + for (auto const& state : objectiveHelper.getMaybeStates()) { + if (indicatorConstraints) { + objectiveValueVariables[state] = lpModel.addContinuousVariable("x_" + std::to_string(objIndex) + "_" + std::to_string(state)); + } else { + objectiveValueVariables[state] = + lpModel.addBoundedContinuousVariable("x_" + std::to_string(objIndex) + "_" + std::to_string(state), + objectiveHelper.getLowerValueBoundAtState(state), objectiveHelper.getUpperValueBoundAtState(state)); + } + } + std::vector reachVars; + if (objectiveHelper.getInfinityCase() == HelperType::InfinityCase::HasNegativeInfinite) { + reachVars.assign(matrix.getRowGroupCount(), {}); + for (auto const& state : objectiveHelper.getRewMinusInfEStates()) { + reachVars[state] = lpModel.addBinaryVariable("c_" + std::to_string(objIndex) + "_" + std::to_string(state)); + } + STORM_LOG_ASSERT(objectiveHelper.getRewMinusInfEStates().get(initialState), ""); + lpModel.update(); + lpModel.addConstraint("", reachVars[initialState] == lpModel.getConstant(storm::utility::one())); + } + lpModel.update(); + for (auto const& state : objectiveHelper.getMaybeStates()) { + bool const requireReachConstraints = + objectiveHelper.getInfinityCase() == HelperType::InfinityCase::HasNegativeInfinite && objectiveHelper.getRewMinusInfEStates().get(state); + for (auto choice : matrix.getRowGroupIndices(state)) { + auto const& choiceVarAsExpression = choiceVariables.at(choice); + STORM_LOG_ASSERT(choiceVarAsExpression.isVariable() || + (!choiceVarAsExpression.containsVariables() && storm::utility::isOne(choiceVarAsExpression.evaluateAsRational())), + "Unexpected kind of choice variable: " << choiceVarAsExpression); + std::vector summands; + if (!indicatorConstraints && choiceVarAsExpression.isVariable()) { + summands.push_back((lpModel.getConstant(storm::utility::one()) - choiceVarAsExpression) * + lpModel.getConstant(objectiveHelper.getUpperValueBoundAtState(state) - objectiveHelper.getLowerValueBoundAtState(state))); + } + if (auto findRes = objectiveHelper.getChoiceRewards().find(choice); findRes != objectiveHelper.getChoiceRewards().end()) { + auto rewExpr = lpModel.getConstant(findRes->second); + if (requireReachConstraints) { + summands.push_back(reachVars[state] * rewExpr); + } else { + summands.push_back(rewExpr); } - if (processed.count(newEcStates) == 0) { - unprocessed.emplace(std::move(newEcStates), std::move(newEcChoices)); + } + for (auto const& succ : matrix.getRow(choice)) { + if (objectiveHelper.getMaybeStates().get(succ.getColumn())) { + summands.push_back(lpModel.getConstant(succ.getValue()) * objectiveValueVariables.at(succ.getColumn())); + } + if (requireReachConstraints && objectiveHelper.getRewMinusInfEStates().get(succ.getColumn())) { + lpModel.addConstraint( + "", reachVars[state] <= reachVars[succ.getColumn()] + lpModel.getConstant(storm::utility::one()) - choiceVarAsExpression); } } + if (summands.empty()) { + summands.push_back(lpModel.getConstant(storm::utility::zero())); + } + if (indicatorConstraints && choiceVarAsExpression.isVariable()) { + auto choiceVar = choiceVarAsExpression.getBaseExpression().asVariableExpression().getVariable(); + lpModel.addIndicatorConstraint("", choiceVar, true, objectiveValueVariables.at(state) <= storm::expressions::sum(summands)); + } else { + lpModel.addConstraint("", objectiveValueVariables.at(state) <= storm::expressions::sum(summands)); + } } - processed.emplace(std::move(currentStates), std::move(currentChoices)); } + return objectiveValueVariables; +} - return processed; +/// Computes the set of problematic MECS with the objective indices that induced them +/// An EC is problematic if its contained in the "maybestates" of an objective and only considers zero-reward choices. +template +auto computeProblematicMecs(storm::storage::SparseMatrix const& matrix, storm::storage::SparseMatrix const& backwardTransitions, + std::vector const& objectiveHelper) { + std::vector>> problMecs; + for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { + auto const& obj = objectiveHelper[objIndex]; + storm::storage::MaximalEndComponentDecomposition objMecs(matrix, backwardTransitions, obj.getMaybeStates(), + obj.getRelevantZeroRewardChoices()); + for (auto& newMec : objMecs) { + bool found = false; + for (auto& problMec : problMecs) { + if (problMec.first == newMec) { + problMec.second.push_back(objIndex); + found = true; + break; + } + } + if (!found) { + problMecs.emplace_back(std::move(newMec), std::vector({objIndex})); + } + } + } + STORM_LOG_DEBUG("Found " << problMecs.size() << " problematic ECs." << std::endl); + return problMecs; } -template -std::map processEc(storm::storage::MaximalEndComponent const& ec, - storm::storage::SparseMatrix const& transitions, std::string const& varNameSuffix, - std::vector const& choiceVars, - storm::solver::LpSolver& lpModel) { - std::map ecStateVars, ecChoiceVars, ecFlowChoiceVars; - // Compute an upper bound on the expected number of visits of the states in this ec. - // First get a lower bound l on the probability of a path that leaves this MEC. 1-l is an upper bound on Pr_s(X F s). - // The desired upper bound is thus 1/(1-(1-l)) = 1/l. See Baier et al., CAV'17 - - // To compute l, we multiply the smallest transition probabilities occurring at each state and MEC-Choice - // as well as the smallest 'exit' probability - // Observe that the actual transition probabilities do not matter for this encoding. - // Hence, we assume that all distributions are uniform to achieve a better (numerically more stable) bound - ValueType lpath = storm::utility::one(); - for (auto const& stateChoices : ec) { - uint64_t maxEntryCount = transitions.getColumnCount(); - // Choices that leave the EC are not considered in the in[s] below. Hence, also do not need to consider them here. - for (auto const& choice : stateChoices.second) { - maxEntryCount = std::max(maxEntryCount, transitions.getRow(choice).getNumberOfEntries()); - } - lpath *= storm::utility::one() / storm::utility::convertNumber(maxEntryCount); +template +auto problematicMecConstraintsExpVisits(storm::solver::LpSolver& lpModel, bool const& indicatorConstraints, bool const& redundantConstraints, + storm::storage::SparseMatrix const& matrix, storm::storage::SparseMatrix const& backwardChoices, + uint64_t mecIndex, storm::storage::MaximalEndComponent const& problematicMec, + std::vector const& relevantObjectiveIndices, + std::vector> const& objectiveValueVariables, + std::vector const& choiceVariables, + UpperBoundsGetterType const& objectiveStateUpperBoundGetter) { + storm::expressions::Expression visitsUpperBound; + if (!indicatorConstraints) { + visitsUpperBound = lpModel.getConstant(VisitingTimesHelper::computeMecVisitsUpperBound(problematicMec, matrix, true)); } - ValueType expVisitsUpperBound = storm::utility::one() / lpath; - STORM_LOG_WARN_COND(expVisitsUpperBound <= storm::utility::convertNumber(1000.0), - "Large upper bound for expected visiting times: " << expVisitsUpperBound); - // create variables - for (auto const& stateChoices : ec) { - ecStateVars.emplace(stateChoices.first, lpModel - .addBoundedIntegerVariable("e" + std::to_string(stateChoices.first) + varNameSuffix, - storm::utility::zero(), storm::utility::one()) - .getExpression()); - for (auto const& choice : stateChoices.second) { - ecChoiceVars.emplace(choice, lpModel - .addBoundedIntegerVariable("ec" + std::to_string(choice) + varNameSuffix, storm::utility::zero(), - storm::utility::one()) - .getExpression()); - ecFlowChoiceVars.emplace( - choice, - lpModel.addBoundedContinuousVariable("f" + std::to_string(choice) + varNameSuffix, storm::utility::zero(), expVisitsUpperBound) - .getExpression()); + + // Create variables and basic lower/upper bounds + storm::storage::BitVector mecChoices(matrix.getRowCount(), false); + std::map expVisitsVars; // z^C_{s,act} + std::map botVars; // z^C_{s,bot} + std::map bsccIndicatorVariables; // b^C_{s} + for (auto const& stateChoices : problematicMec) { + auto const state = stateChoices.first; + auto bsccIndicatorVar = lpModel.addBinaryVariable("b_" + std::to_string(mecIndex) + "_" + std::to_string(state)); + bsccIndicatorVariables.emplace(state, bsccIndicatorVar.getExpression()); + std::string visitsVarPref = "z_" + std::to_string(mecIndex) + "_"; + auto stateBotVisitsVar = + lpModel.addLowerBoundedContinuousVariable(visitsVarPref + std::to_string(state) + "bot", storm::utility::zero()).getExpression(); + botVars.emplace(state, stateBotVisitsVar); + lpModel.update(); + if (indicatorConstraints) { + lpModel.addIndicatorConstraint("", bsccIndicatorVar, false, stateBotVisitsVar <= lpModel.getConstant(storm::utility::zero())); + } else { + lpModel.addConstraint("", stateBotVisitsVar <= bsccIndicatorVar.getExpression() * visitsUpperBound); + } + for (auto choice : matrix.getRowGroupIndices(state)) { + auto stateActionVisitsVar = + lpModel.addLowerBoundedContinuousVariable(visitsVarPref + std::to_string(choice), storm::utility::zero()).getExpression(); + lpModel.update(); + if (indicatorConstraints) { + if (auto const& a = choiceVariables[choice]; a.isVariable()) { + auto aVar = a.getBaseExpression().asVariableExpression().getVariable(); + lpModel.addIndicatorConstraint("", aVar, false, stateActionVisitsVar <= lpModel.getConstant(storm::utility::zero())); + } + } else { + lpModel.addConstraint("", stateActionVisitsVar <= choiceVariables[choice] * visitsUpperBound); + } + expVisitsVars.emplace(choice, stateActionVisitsVar); + } + for (auto const& ecChoice : stateChoices.second) { + mecChoices.set(ecChoice, true); + } + for (auto const& objIndex : relevantObjectiveIndices) { + if (indicatorConstraints) { + lpModel.addIndicatorConstraint("", bsccIndicatorVar, true, + objectiveValueVariables[objIndex][state] <= lpModel.getConstant(storm::utility::zero())); + } else { + auto const upperBnd = lpModel.getConstant(objectiveStateUpperBoundGetter(objIndex, state)); + lpModel.addConstraint("", objectiveValueVariables[objIndex][state] <= upperBnd - upperBnd * bsccIndicatorVar.getExpression()); + } } } - // create constraints - std::map> ins, outs, ecIns; - for (auto const& stateChoices : ec) { - std::vector ecChoiceVarsAtState; - std::vector out; - for (auto const& choice : stateChoices.second) { - if (choiceVars[choice].isInitialized()) { - lpModel.addConstraint("", ecChoiceVars[choice] <= choiceVars[choice]); - lpModel.addConstraint("", ecFlowChoiceVars[choice] <= lpModel.getConstant(expVisitsUpperBound) * choiceVars[choice]); - } - ecChoiceVarsAtState.push_back(ecChoiceVars[choice]); - out.push_back(ecFlowChoiceVars[choice]); - ValueType fakeProbability = - storm::utility::one() / storm::utility::convertNumber(transitions.getRow(choice).getNumberOfEntries()); - for (auto const& transition : transitions.getRow(choice)) { - lpModel.addConstraint("", ecChoiceVars[choice] <= ecStateVars[transition.getColumn()]); - ins[transition.getColumn()].push_back(lpModel.getConstant(fakeProbability) * ecFlowChoiceVars[choice]); - ecIns[transition.getColumn()].push_back(ecChoiceVars[choice]); + // Create visits constraints + auto const initProb = lpModel.getConstant(storm::utility::one() / storm::utility::convertNumber(problematicMec.size())); + std::vector outVisitsSummands; + for (auto const& stateChoices : problematicMec) { + auto const state = stateChoices.first; + auto const& choices = stateChoices.second; + auto const& stateBotVar = botVars.at(state); + outVisitsSummands.push_back(stateBotVar); + std::vector stateVisitsSummands; + stateVisitsSummands.push_back(stateBotVar); + for (auto choice : matrix.getRowGroupIndices(state)) { + auto const& choiceVisitsVar = expVisitsVars.at(choice); + stateVisitsSummands.push_back(choiceVisitsVar); + if (choices.count(choice) != 0) { + if (redundantConstraints) { + for (auto const& postElem : matrix.getRow(choice)) { + if (storm::utility::isZero(postElem.getValue())) { + continue; + } + auto succ = postElem.getColumn(); + lpModel.addConstraint("", bsccIndicatorVariables.at(state) + choiceVariables.at(choice) <= + lpModel.getConstant(storm::utility::one()) + bsccIndicatorVariables.at(succ)); + } + } + } else { + outVisitsSummands.push_back(choiceVisitsVar); } } - lpModel.addConstraint("", ecStateVars[stateChoices.first] == storm::expressions::sum(ecChoiceVarsAtState)); - out.push_back(lpModel.getConstant(expVisitsUpperBound) * ecStateVars[stateChoices.first]); - // Iterate over choices that leave the ec - for (uint64_t choice = transitions.getRowGroupIndices()[stateChoices.first]; choice < transitions.getRowGroupIndices()[stateChoices.first + 1]; - ++choice) { - if (stateChoices.second.find(choice) == stateChoices.second.end()) { - assert(choiceVars[choice].isInitialized()); - out.push_back(lpModel.getConstant(expVisitsUpperBound) * choiceVars[choice]); + for (auto const& preEntry : backwardChoices.getRow(state)) { + uint64_t const preChoice = preEntry.getColumn(); + if (mecChoices.get(preChoice)) { + ValueType preProb = + storm::utility::one() / storm::utility::convertNumber(matrix.getRow(preChoice).getNumberOfEntries()); + stateVisitsSummands.push_back(lpModel.getConstant(-preProb) * expVisitsVars.at(preChoice)); } } - outs.emplace(stateChoices.first, out); + lpModel.addConstraint("", storm::expressions::sum(stateVisitsSummands) == initProb); } - uint64_t numStates = std::distance(ec.begin(), ec.end()); - for (auto const& stateChoices : ec) { - auto in = ins.find(stateChoices.first); - STORM_LOG_ASSERT(in != ins.end(), "ec state does not seem to have an incoming transition."); - // Assume a uniform initial distribution - in->second.push_back(lpModel.getConstant(storm::utility::one() / storm::utility::convertNumber(numStates))); - auto out = outs.find(stateChoices.first); - STORM_LOG_ASSERT(out != outs.end(), "out flow of ec state was not set."); - lpModel.addConstraint("", storm::expressions::sum(in->second) <= storm::expressions::sum(out->second)); - auto ecIn = ecIns.find(stateChoices.first); - STORM_LOG_ASSERT(ecIn != ecIns.end(), "ec state does not seem to have an incoming transition."); - lpModel.addConstraint("", storm::expressions::sum(ecIn->second) >= ecStateVars[stateChoices.first]); + lpModel.addConstraint("", storm::expressions::sum(outVisitsSummands) == lpModel.getConstant(storm::utility::one())); +} + +template +auto problematicMecConstraintsOrder(storm::solver::LpSolver& lpModel, bool const& indicatorConstraints, bool const& redundantConstraints, + storm::storage::SparseMatrix const& matrix, uint64_t mecIndex, + storm::storage::MaximalEndComponent const& problematicMec, std::vector const& relevantObjectiveIndices, + std::vector const& choiceVariables, + std::vector> const& objectiveValueVariables, + UpperBoundsGetterType const& objectiveStateUpperBoundGetter) { + // Create bscc indicator and order variables with basic lower/upper bounds + storm::storage::BitVector mecChoices(matrix.getRowCount(), false); + std::map bsccIndicatorVariables; // b^C_{s} + std::map orderVariables; // r^C_{s} + for (auto const& stateChoices : problematicMec) { + auto const state = stateChoices.first; + auto bsccIndicatorVar = lpModel.addBinaryVariable("b_" + std::to_string(mecIndex) + "_" + std::to_string(state)); + bsccIndicatorVariables.emplace(state, bsccIndicatorVar.getExpression()); + auto orderVar = lpModel + .addBoundedContinuousVariable("r_" + std::to_string(mecIndex) + "_" + std::to_string(state), storm::utility::zero(), + storm::utility::one()) + .getExpression(); + lpModel.update(); + orderVariables.emplace(state, orderVar); + for (auto const& ecChoice : stateChoices.second) { + mecChoices.set(ecChoice, true); + } + for (auto const& objIndex : relevantObjectiveIndices) { + if (indicatorConstraints) { + lpModel.addIndicatorConstraint("", bsccIndicatorVar, true, + objectiveValueVariables[objIndex][state] <= lpModel.getConstant(storm::utility::zero())); + } else { + auto const upperBnd = lpModel.getConstant(objectiveStateUpperBoundGetter(objIndex, state)); + lpModel.addConstraint("", objectiveValueVariables[objIndex][state] <= upperBnd - upperBnd * bsccIndicatorVar.getExpression()); + } + } } - return ecStateVars; + // Create order constraints + auto const minDiff = lpModel.getConstant(-storm::utility::one() / storm::utility::convertNumber(problematicMec.size())); + for (auto const& stateChoices : problematicMec) { + auto const state = stateChoices.first; + auto const& choices = stateChoices.second; + auto const& bsccIndicatorVar = bsccIndicatorVariables.at(state); + auto const& orderVar = orderVariables.at(state); + for (auto choice : choices) { + auto const& choiceVariable = choiceVariables.at(choice); + std::vector choiceConstraint; + choiceConstraint.push_back(bsccIndicatorVar); + std::string const transSelectPrefix = "d_" + std::to_string(mecIndex) + "_" + std::to_string(choice) + "_"; + for (auto const& postElem : matrix.getRow(choice)) { + if (storm::utility::isZero(postElem.getValue())) { + continue; + } + auto succ = postElem.getColumn(); + if (redundantConstraints) { + lpModel.addConstraint( + "", bsccIndicatorVar + choiceVariable <= lpModel.getConstant(storm::utility::one()) + bsccIndicatorVariables.at(succ)); + } + auto transVar = lpModel.addBinaryVariable(transSelectPrefix + std::to_string(succ)).getExpression(); + lpModel.update(); + choiceConstraint.push_back(transVar); + lpModel.addConstraint("", transVar <= choiceVariable); + lpModel.addConstraint("", orderVar <= minDiff + orderVariables.at(succ) + lpModel.getConstant(storm::utility::one()) - transVar); + } + lpModel.addConstraint("", choiceVariable <= storm::expressions::sum(choiceConstraint)); + } + } } -template -bool DeterministicSchedsLpChecker::processEndComponents(std::vector>& ecVars) { - uint64_t ecCounter = 0; - auto backwardTransitions = model.getBackwardTransitions(); - - // Get the choices that do not induce a value (i.e. reward) for all objectives. - // Only MECS consisting of these choices are relevant - storm::storage::BitVector choicesWithValueZero(model.getNumberOfChoices(), true); - for (auto const& objHelper : objectiveHelper) { - for (auto const& value : objHelper.getChoiceValueOffsets()) { - STORM_LOG_ASSERT(!storm::utility::isZero(value.second), "Expected non-zero choice-value offset."); - choicesWithValueZero.set(value.first, false); +template +std::vector expVisitsConstraints(storm::solver::LpSolver& lpModel, bool const& indicatorConstraints, + storm::storage::SparseMatrix const& matrix, + storm::storage::SparseMatrix const& backwardTransitions, + storm::storage::SparseMatrix const& backwardChoices, uint64_t initialState, + std::vector const& objectiveHelper, + std::vector const& choiceVariables) { + auto objHelpIt = objectiveHelper.begin(); + storm::storage::BitVector anyMaybeStates = objHelpIt->getMaybeStates(); + for (++objHelpIt; objHelpIt != objectiveHelper.end(); ++objHelpIt) { + anyMaybeStates |= objHelpIt->getMaybeStates(); + } + storm::storage::BitVector allZeroRewardChoices(matrix.getRowCount(), true); + for (auto const& oh : objectiveHelper) { + for (auto const& rew : oh.getChoiceRewards()) { + assert(!storm::utility::isZero(rew.second)); + allZeroRewardChoices.set(rew.first, false); } } - storm::storage::MaximalEndComponentDecomposition mecs(model.getTransitionMatrix(), backwardTransitions, - storm::storage::BitVector(model.getNumberOfStates(), true), choicesWithValueZero); + storm::storage::MaximalEndComponentDecomposition mecs(matrix, backwardTransitions, anyMaybeStates, allZeroRewardChoices); + storm::storage::BitVector mecStates(matrix.getRowGroupCount(), false); for (auto const& mec : mecs) { - // For each objective we might need to split this mec into several subECs, if the objective yields a non-zero scheduler-independent state value for some - // states of this ec. However, note that this split happens objective-wise which is why we can not consider a subsystem in the mec-decomposition above - std::map, std::vector> excludedStatesToObjIndex; - bool mecContainsSchedulerDependentValue = false; - for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { - std::set excludedStates; - for (auto const& stateChoices : mec) { - auto schedIndValueIt = objectiveHelper[objIndex].getSchedulerIndependentStateValues().find(stateChoices.first); - if (schedIndValueIt == objectiveHelper[objIndex].getSchedulerIndependentStateValues().end()) { - mecContainsSchedulerDependentValue = true; - } else if (!storm::utility::isZero(schedIndValueIt->second)) { - excludedStates.insert(stateChoices.first); + for (auto const& sc : mec) { + mecStates.set(sc.first, true); + } + } + std::vector maxVisits; + if (!indicatorConstraints) { + maxVisits = VisitingTimesHelper::computeUpperBoundsOnExpectedVisitingTimes(anyMaybeStates, matrix, backwardTransitions); + } + + // Create variables and basic bounds + std::vector choiceVisitsVars(matrix.getRowCount()), botVisitsVars(matrix.getRowGroupCount()), + bsccVars(matrix.getRowGroupCount()); + for (auto state : anyMaybeStates) { + assert(indicatorConstraints || maxVisits[state] >= storm::utility::zero()); + for (auto choice : matrix.getRowGroupIndices(state)) { + choiceVisitsVars[choice] = + lpModel.addLowerBoundedContinuousVariable("y_" + std::to_string(choice), storm::utility::zero()).getExpression(); + lpModel.update(); + if (indicatorConstraints) { + if (auto const& a = choiceVariables[choice]; a.isVariable()) { + auto aVar = a.getBaseExpression().asVariableExpression().getVariable(); + lpModel.addIndicatorConstraint("", aVar, false, choiceVisitsVars[choice] <= lpModel.getConstant(storm::utility::zero())); } + } else { + lpModel.addConstraint("", choiceVisitsVars[choice] <= choiceVariables.at(choice) * lpModel.getConstant(maxVisits[state])); + } + } + if (mecStates.get(state)) { + bsccVars[state] = lpModel.addBinaryVariable("b_" + std::to_string(state)).getExpression(); + botVisitsVars[state] = + lpModel.addLowerBoundedContinuousVariable("y_" + std::to_string(state) + "bot", storm::utility::zero()).getExpression(); + lpModel.update(); + if (indicatorConstraints) { + lpModel.addIndicatorConstraint("", bsccVars[state].getBaseExpression().asVariableExpression().getVariable(), false, + botVisitsVars[state] <= lpModel.getConstant(storm::utility::zero())); + } else { + lpModel.addConstraint("", botVisitsVars[state] <= bsccVars[state] * lpModel.getConstant(maxVisits[state])); + } + } + } + + // Add expected visiting times constraints + auto notMaybe = ~anyMaybeStates; + std::vector outSummands; + for (auto state : anyMaybeStates) { + std::vector visitsSummands; + if (mecStates.get(state)) { + visitsSummands.push_back(-botVisitsVars[state]); + outSummands.push_back(botVisitsVars[state]); + } + for (auto choice : matrix.getRowGroupIndices(state)) { + visitsSummands.push_back(-choiceVisitsVars[choice]); + if (auto outProb = matrix.getConstrainedRowSum(choice, notMaybe); !storm::utility::isZero(outProb)) { + outSummands.push_back(lpModel.getConstant(outProb) * choiceVisitsVars[choice]); } - excludedStatesToObjIndex[excludedStates].push_back(objIndex); - } - - // Skip this mec if all state values are independent of the scheduler (e.g. no reward is reachable from here). - if (mecContainsSchedulerDependentValue) { - for (auto const& exclStates : excludedStatesToObjIndex) { - if (exclStates.first.empty()) { - auto varsForMec = processEc(mec, model.getTransitionMatrix(), "", choiceVariables, *lpModel); - ++ecCounter; - for (auto const& stateVar : varsForMec) { - for (auto const& objIndex : exclStates.second) { - ecVars[objIndex][stateVar.first] = stateVar.second; + } + if (state == initialState) { + visitsSummands.push_back(lpModel.getConstant(storm::utility::one())); + } + for (auto const& preEntry : backwardChoices.getRow(state)) { + assert(choiceVisitsVars[preEntry.getColumn()].isInitialized()); + visitsSummands.push_back(lpModel.getConstant(preEntry.getValue()) * choiceVisitsVars[preEntry.getColumn()]); + } + lpModel.addConstraint("", storm::expressions::sum(visitsSummands) == lpModel.getConstant(storm::utility::zero())); + } + lpModel.addConstraint("", storm::expressions::sum(outSummands) == lpModel.getConstant(storm::utility::one())); + + // Add bscc constraints + for (auto const& mec : mecs) { + for (auto const& stateChoices : mec) { + auto const& state = stateChoices.first; + for (auto choice : matrix.getRowGroupIndices(state)) { + if (stateChoices.second.count(choice) != 0) { + for (auto const& succ : matrix.getRow(choice)) { + if (storm::utility::isZero(succ.getValue())) { + continue; } + assert(mecStates.get(succ.getColumn())); + lpModel.addConstraint("", bsccVars[state] <= bsccVars[succ.getColumn()] + lpModel.getConstant(storm::utility::one()) - + choiceVariables[choice]); } } else { - // Compute sub-end components - storm::storage::BitVector subEcStates(model.getNumberOfStates(), false), subEcChoices(model.getNumberOfChoices(), false); - for (auto const& stateChoices : mec) { - if (exclStates.first.count(stateChoices.first) == 0) { - subEcStates.set(stateChoices.first, true); - for (auto const& choice : stateChoices.second) { - subEcChoices.set(choice, true); - } - } - } - storm::storage::MaximalEndComponentDecomposition subEcs(model.getTransitionMatrix(), backwardTransitions, subEcStates, - subEcChoices); - for (auto const& subEc : subEcs) { - auto varsForSubEc = - processEc(subEc, model.getTransitionMatrix(), "o" + std::to_string(*exclStates.second.begin()), choiceVariables, *lpModel); - ++ecCounter; - for (auto const& stateVar : varsForSubEc) { - for (auto const& objIndex : exclStates.second) { - ecVars[objIndex][stateVar.first] = stateVar.second; - } - } - } + lpModel.addConstraint("", bsccVars[state] <= lpModel.getConstant(storm::utility::one()) - choiceVariables[choice]); } } } } - bool hasEndComponents = - ecCounter > 0 || storm::utility::graph::checkIfECWithChoiceExists(model.getTransitionMatrix(), backwardTransitions, - storm::storage::BitVector(model.getNumberOfStates(), true), ~choicesWithValueZero); - STORM_LOG_WARN_COND(!hasEndComponents, "Processed " << ecCounter << " End components."); - return hasEndComponents; + + // Add objective values + std::vector objectiveValueVariables; + for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { + if (objectiveHelper[objIndex].getMaybeStates().get(initialState)) { + objectiveValueVariables.push_back(lpModel.addUnboundedContinuousVariable("x_" + std::to_string(objIndex)).getExpression()); + lpModel.update(); + std::vector summands; + for (auto const& objRew : objectiveHelper[objIndex].getChoiceRewards()) { + assert(choiceVisitsVars[objRew.first].isInitialized()); + summands.push_back(choiceVisitsVars[objRew.first] * lpModel.getConstant(objRew.second)); + } + lpModel.addConstraint("", objectiveValueVariables.back() == storm::expressions::sum(summands)); + } else { + objectiveValueVariables.push_back(lpModel.getConstant(objectiveHelper[objIndex].getConstantInitialStateValue())); + } + } + return objectiveValueVariables; +} + +template +bool useFlowEncoding(storm::Environment const& env, std::vector const& objectiveHelper) { + bool supportsFlowEncoding = std::all_of(objectiveHelper.begin(), objectiveHelper.end(), [](auto const& h) { return h.isTotalRewardObjective(); }); + switch (env.modelchecker().multi().getEncodingType()) { + case storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Auto: + return supportsFlowEncoding; + case storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow: + STORM_LOG_THROW(supportsFlowEncoding, storm::exceptions::InvalidOperationException, + "Flow encoding only applicable if all objectives are (transformable to) total reward objectives."); + return true; + default: + return false; + } } template void DeterministicSchedsLpChecker::initializeLpModel(Environment const& env) { STORM_LOG_INFO("Initializing LP model with " << model.getNumberOfStates() << " states."); - if (env.modelchecker().multi().getEncodingType() == storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Auto) { - flowEncoding = true; - for (auto const& helper : objectiveHelper) { - if (!helper.isTotalRewardObjective()) { - flowEncoding = false; - } - } - } else { - flowEncoding = env.modelchecker().multi().getEncodingType() == storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow; - } + flowEncoding = useFlowEncoding(env, objectiveHelper); STORM_LOG_INFO("Using " << (flowEncoding ? "flow" : "classical") << " encoding.\n"); - - uint64_t numStates = model.getNumberOfStates(); uint64_t initialState = *model.getInitialStates().begin(); + auto backwardTransitions = model.getBackwardTransitions(); + auto backwardChoices = model.getTransitionMatrix().transpose(); STORM_LOG_WARN_COND(!storm::settings::getModule().isLpSolverSetFromDefaultValue() || storm::settings::getModule().getLpSolver() == storm::solver::LpSolverType::Gurobi, "The selected MILP solver might not perform well. Consider installing / using Gurobi."); @@ -375,271 +589,38 @@ void DeterministicSchedsLpChecker::initializeLpMod lpModel->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); initialStateResults.clear(); - auto zero = lpModel->getConstant(storm::utility::zero()); - auto one = lpModel->getConstant(storm::utility::one()); - auto const& groups = model.getTransitionMatrix().getRowGroupIndices(); // Create choice variables. - choiceVariables.reserve(model.getNumberOfChoices()); - for (uint64_t state = 0; state < numStates; ++state) { - uint64_t numChoices = model.getNumberOfChoices(state); - if (numChoices == 1) { - choiceVariables.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()); - choiceVariables.push_back(localChoices.back()); - } - lpModel->addConstraint("", storm::expressions::sum(localChoices) == one); - } - } - // Create ec Variables for each state/objective - std::vector> ecVars(objectiveHelper.size(), - std::vector(model.getNumberOfStates())); - bool hasEndComponents = processEndComponents(ecVars); - + choiceVariables = createChoiceVariables(*lpModel, model.getTransitionMatrix()); if (flowEncoding) { - storm::storage::BitVector bottomStates(model.getNumberOfStates(), true); - for (auto const& helper : objectiveHelper) { - STORM_LOG_THROW(helper.isTotalRewardObjective(), storm::exceptions::InvalidOperationException, - "The given type of encoding is only supported if the objectives can be reduced to total reward objectives."); - storm::storage::BitVector objBottomStates(model.getNumberOfStates(), false); - for (auto const& stateVal : helper.getSchedulerIndependentStateValues()) { - STORM_LOG_THROW(storm::utility::isZero(stateVal.second), storm::exceptions::InvalidOperationException, - "Non-zero constant state-values not allowed for flow encoding."); - objBottomStates.set(stateVal.first, true); - } - bottomStates &= objBottomStates; - } - storm::storage::BitVector nonBottomStates = ~bottomStates; - STORM_LOG_TRACE("Found " << bottomStates.getNumberOfSetBits() << " bottom states."); - - // Compute upper bounds for each state - std::vector visitingTimesUpperBounds = DeterministicSchedsObjectiveHelper::computeUpperBoundOnExpectedVisitingTimes( - model.getTransitionMatrix(), bottomStates, nonBottomStates, hasEndComponents); - ValueType largestUpperBound = *std::max_element(visitingTimesUpperBounds.begin(), visitingTimesUpperBounds.end()); - STORM_LOG_WARN_COND(largestUpperBound < storm::utility::convertNumber(1e5), - "Found a large upper bound '" << storm::utility::convertNumber(largestUpperBound) - << "' in the LP encoding. This might trigger numerical instabilities."); - // create choiceValue variables and assert deterministic ones. - std::vector choiceValVars(model.getNumberOfChoices()); - for (auto state : nonBottomStates) { - for (uint64_t globalChoice = groups[state]; globalChoice < groups[state + 1]; ++globalChoice) { - choiceValVars[globalChoice] = - lpModel - ->addBoundedContinuousVariable("y" + std::to_string(globalChoice), storm::utility::zero(), visitingTimesUpperBounds[state]) - .getExpression(); - if (model.getNumberOfChoices(state) > 1) { - ; - lpModel->addConstraint( - "", choiceValVars[globalChoice] <= lpModel->getConstant(visitingTimesUpperBounds[state]) * choiceVariables[globalChoice]); - } - } - } - // create EC 'slack' variables for states that lie in an ec - std::vector ecValVars(model.getNumberOfStates()); - if (hasEndComponents) { - for (auto state : nonBottomStates) { - // For the in-out-encoding, all objectives have the same ECs (because there are no non-zero scheduler independend state values). - // Hence, we only care for the variables of the first objective. - if (ecVars.front()[state].isInitialized()) { - ecValVars[state] = - lpModel->addBoundedContinuousVariable("z" + std::to_string(state), storm::utility::zero(), visitingTimesUpperBounds[state]) - .getExpression(); - lpModel->addConstraint("", ecValVars[state] <= lpModel->getConstant(visitingTimesUpperBounds[state]) * ecVars.front()[state]); - } - } - } - // Get 'in' and 'out' expressions - std::vector bottomStatesIn; - std::vector> ins(numStates), outs(numStates); - ins[initialState].push_back(one); - for (auto state : nonBottomStates) { - for (uint64_t globalChoice = groups[state]; globalChoice < groups[state + 1]; ++globalChoice) { - for (auto const& transition : model.getTransitionMatrix().getRow(globalChoice)) { - uint64_t successor = transition.getColumn(); - storm::expressions::Expression exp = lpModel->getConstant(transition.getValue()) * choiceValVars[globalChoice]; - if (bottomStates.get(successor)) { - bottomStatesIn.push_back(exp); - } else { - ins[successor].push_back(exp); - } - } - outs[state].push_back(choiceValVars[globalChoice]); - } - if (ecValVars[state].isInitialized()) { - outs[state].push_back(ecValVars[state]); - bottomStatesIn.push_back(ecValVars[state]); - } - } - - // Assert 'in == out' at each state - for (auto state : nonBottomStates) { - lpModel->addConstraint("", storm::expressions::sum(ins[state]) == storm::expressions::sum(outs[state])); - } - // Assert the sum for the bottom states - lpModel->addConstraint("", storm::expressions::sum(bottomStatesIn) == one); - - // create initial state results for each objective + initialStateResults = expVisitsConstraints(*lpModel, env.modelchecker().multi().getUseIndicatorConstraints(), model.getTransitionMatrix(), + backwardTransitions, backwardChoices, initialState, objectiveHelper, choiceVariables); + } else { + std::vector> objectiveValueVariables; + initialStateResults.clear(); for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { - auto choiceValueOffsets = objectiveHelper[objIndex].getChoiceValueOffsets(); - std::vector objValue; - for (auto state : nonBottomStates) { - for (uint64_t globalChoice = groups[state]; globalChoice < groups[state + 1]; ++globalChoice) { - auto choiceValueIt = choiceValueOffsets.find(globalChoice); - if (choiceValueIt != choiceValueOffsets.end()) { - assert(!storm::utility::isZero(choiceValueIt->second)); - objValue.push_back(lpModel->getConstant(choiceValueIt->second) * choiceValVars[globalChoice]); - } - } - } - ValueType lowerBound = objectiveHelper[objIndex].getLowerValueBoundAtState(env, initialState); - ValueType upperBound = objectiveHelper[objIndex].getUpperValueBoundAtState(env, initialState); - storm::expressions::Expression objValueVariable; - if (lowerBound == upperBound) { - // GLPK does not like point-intervals...... - objValueVariable = lpModel->getConstant(lowerBound); - } else { - objValueVariable = lpModel->addBoundedContinuousVariable("x" + std::to_string(objIndex), lowerBound, upperBound).getExpression(); - } - if (objValue.empty()) { - lpModel->addConstraint("", objValueVariable == zero); + if (objectiveHelper[objIndex].getMaybeStates().get(initialState)) { + objectiveValueVariables.push_back(classicConstraints(*lpModel, env.modelchecker().multi().getUseIndicatorConstraints(), + model.getTransitionMatrix(), initialState, objIndex, objectiveHelper[objIndex], + choiceVariables)); + initialStateResults.push_back(objectiveValueVariables.back()[initialState]); } else { - lpModel->addConstraint("", objValueVariable == storm::expressions::sum(objValue)); + initialStateResults.push_back(lpModel->getConstant(objectiveHelper[objIndex].getConstantInitialStateValue())); } - initialStateResults.push_back(objValueVariable); } - } else { - // 'classic' backward encoding. - for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { - STORM_LOG_WARN_COND(objectiveHelper[objIndex].getLargestUpperBound(env) < storm::utility::convertNumber(1e5), - "Found a large upper value bound '" - << storm::utility::convertNumber(objectiveHelper[objIndex].getLargestUpperBound(env)) - << "' in the LP encoding. This might trigger numerical instabilities."); - auto const& schedulerIndependentStates = objectiveHelper[objIndex].getSchedulerIndependentStateValues(); - // Create state variables and store variables of ecs which contain a state with a scheduler independent value - std::vector stateVars; - stateVars.reserve(numStates); - for (uint64_t state = 0; state < numStates; ++state) { - auto valIt = schedulerIndependentStates.find(state); - if (valIt == schedulerIndependentStates.end()) { - ValueType lowerBound = objectiveHelper[objIndex].getLowerValueBoundAtState(env, state); - ValueType upperBound = objectiveHelper[objIndex].getUpperValueBoundAtState(env, state); - if (lowerBound == upperBound) { - // glpk does not like variables with point-interval bounds... - if (objectiveHelper[objIndex].minimizing()) { - stateVars.push_back(lpModel->getConstant(-lowerBound)); - } else { - stateVars.push_back(lpModel->getConstant(lowerBound)); - } - } else { - if (objectiveHelper[objIndex].minimizing()) { - stateVars.push_back( - lpModel->addBoundedContinuousVariable("x" + std::to_string(objIndex) + "_" + std::to_string(state), -upperBound, -lowerBound) - .getExpression()); - } else { - stateVars.push_back( - lpModel->addBoundedContinuousVariable("x" + std::to_string(objIndex) + "_" + std::to_string(state), lowerBound, upperBound) - .getExpression()); - } - } - } else { - ValueType value = valIt->second; - if (objectiveHelper[objIndex].minimizing()) { - value = -value; - } - stateVars.push_back(lpModel->getConstant(value)); - } - if (state == initialState) { - 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; - } - storm::expressions::Expression stateValue; - uint64_t numChoices = model.getNumberOfChoices(state); - uint64_t choiceOffset = groups[state]; - storm::expressions::Expression upperValueBoundAtState = lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state)); - for (uint64_t choice = 0; choice < numChoices; ++choice) { - storm::expressions::Expression choiceValue; - auto valIt = choiceValueOffsets.find(choiceOffset + choice); - if (valIt != choiceValueOffsets.end()) { - if (objectiveHelper[objIndex].minimizing()) { - choiceValue = lpModel->getConstant(-valIt->second); - } else { - choiceValue = lpModel->getConstant(valIt->second); - } - } - for (auto const& transition : model.getTransitionMatrix().getRow(state, choice)) { - if (!storm::utility::isZero(transition.getValue())) { - 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) { - stateValue = choiceValue; - } else { - uint64_t globalChoiceIndex = groups[state] + choice; - storm::expressions::Expression choiceValVar; - if (objectiveHelper[objIndex].minimizing()) { - choiceValVar = lpModel - ->addBoundedContinuousVariable( - "x" + std::to_string(objIndex) + "_" + std::to_string(state) + "_" + std::to_string(choice), - -objectiveHelper[objIndex].getUpperValueBoundAtState(env, state), storm::utility::zero()) - .getExpression(); - } else { - choiceValVar = lpModel - ->addBoundedContinuousVariable( - "x" + std::to_string(objIndex) + "_" + std::to_string(state) + "_" + std::to_string(choice), - storm::utility::zero(), objectiveHelper[objIndex].getUpperValueBoundAtState(env, state)) - .getExpression(); - } - lpModel->addConstraint("", choiceValVar <= choiceValue); - if (objectiveHelper[objIndex].minimizing()) { - lpModel->addConstraint("", choiceValVar <= -upperValueBoundAtState * (one - choiceVariables[globalChoiceIndex])); - } else { - lpModel->addConstraint("", choiceValVar <= upperValueBoundAtState * choiceVariables[globalChoiceIndex]); - } - if (choice == 0) { - stateValue = choiceValVar; - } else { - stateValue = stateValue + choiceValVar; - } - } - } - stateValue.simplify().reduceNesting(); - if (objectiveHelper[objIndex].minimizing()) { - lpModel->addConstraint( - "", stateVars[state] <= - stateValue + (lpModel->getConstant(storm::utility::convertNumber(numChoices - 1)) * upperValueBoundAtState)); - } else { - lpModel->addConstraint("", stateVars[state] <= stateValue); - } - if (numChoices > 1 && hasEndComponents) { - auto& ecVar = ecVars[objIndex][state]; - if (ecVar.isInitialized()) { - // if this state is part of an ec, make sure to assign a value of zero. - if (objectiveHelper[objIndex].minimizing()) { - // This part is optional - lpModel->addConstraint( - "", stateVars[state] >= (ecVar - one) * lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state))); - } else { - lpModel->addConstraint( - "", stateVars[state] <= (one - ecVar) * lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state))); - } - } - } + auto problematicMecs = computeProblematicMecs(model.getTransitionMatrix(), backwardTransitions, objectiveHelper); + uint64_t mecIndex = 0; + auto upperBoundsGetter = [&](uint64_t objIndex, uint64_t state) -> ValueType { return objectiveHelper[objIndex].getUpperValueBoundAtState(state); }; + for (auto const& mecObj : problematicMecs) { + if (env.modelchecker().multi().getUseBsccOrderEncoding()) { + problematicMecConstraintsOrder(*lpModel, env.modelchecker().multi().getUseIndicatorConstraints(), + env.modelchecker().multi().getUseRedundantBsccConstraints(), model.getTransitionMatrix(), mecIndex, mecObj.first, + mecObj.second, choiceVariables, objectiveValueVariables, upperBoundsGetter); + } else { + problematicMecConstraintsExpVisits(*lpModel, env.modelchecker().multi().getUseIndicatorConstraints(), + env.modelchecker().multi().getUseRedundantBsccConstraints(), model.getTransitionMatrix(), backwardChoices, + mecIndex, mecObj.first, mecObj.second, objectiveValueVariables, choiceVariables, upperBoundsGetter); } + ++mecIndex; } } lpModel->update(); @@ -785,36 +766,27 @@ void DeterministicSchedsLpChecker::checkRecursive( template typename DeterministicSchedsLpChecker::Point DeterministicSchedsLpChecker::validateCurrentModel( Environment const& env) const { - storm::storage::Scheduler scheduler(model.getNumberOfStates()); + storm::storage::BitVector selectedChoices(model.getNumberOfChoices(), false); for (uint64_t state = 0; state < model.getNumberOfStates(); ++state) { - uint64_t numChoices = model.getNumberOfChoices(state); - if (numChoices == 1) { - scheduler.setChoice(0, state); + auto choices = model.getTransitionMatrix().getRowGroupIndices(state); + if (choices.size() == 1) { + selectedChoices.set(*choices.begin()); } else { - uint64_t globalChoiceOffset = model.getTransitionMatrix().getRowGroupIndices()[state]; bool choiceFound = false; - for (uint64_t localChoice = 0; localChoice < numChoices; ++localChoice) { - bool localChoiceEnabled = false; - localChoiceEnabled = - (lpModel->getIntegerValue(choiceVariables[globalChoiceOffset + localChoice].getBaseExpression().asVariableExpression().getVariable()) == 1); - if (localChoiceEnabled) { + for (auto choice : choices) { + assert(choiceVariables[choice].isVariable()); + if (lpModel->getBinaryValue(choiceVariables[choice].getBaseExpression().asVariableExpression().getVariable())) { STORM_LOG_THROW(!choiceFound, storm::exceptions::UnexpectedException, "Multiple choices selected at state " << state); - scheduler.setChoice(localChoice, state); + selectedChoices.set(choice, true); choiceFound = true; } } - STORM_LOG_THROW(choiceFound, storm::exceptions::UnexpectedException, "No choice selected at state " << state); } } - bool dropUnreachableStates = true; - bool preserveModelType = true; - auto inducedModel = model.applyScheduler(scheduler, dropUnreachableStates, preserveModelType)->template as(); + Point inducedPoint; for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { - ValueType inducedValue = objectiveHelper[objIndex].evaluateOnModel(env, *inducedModel); - if (objectiveHelper[objIndex].minimizing()) { - inducedValue = -inducedValue; - } + ValueType inducedValue = objectiveHelper[objIndex].evaluateScheduler(env, selectedChoices); inducedPoint.push_back(storm::utility::convertNumber(inducedValue)); // If this objective has weight zero, the lp solution is not necessarily correct if (!storm::utility::isZero(currentWeightVector[objIndex])) { @@ -832,6 +804,4 @@ template class DeterministicSchedsLpChecker, template class DeterministicSchedsLpChecker, storm::RationalNumber>; template class DeterministicSchedsLpChecker, storm::RationalNumber>; template class DeterministicSchedsLpChecker, storm::RationalNumber>; -} // namespace multiobjective -} // namespace modelchecker -} // namespace storm +} // namespace storm::modelchecker::multiobjective \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h index 820350da68..b01ff68d4c 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include "storm/modelchecker/multiobjective/Objective.h" @@ -16,6 +17,10 @@ class Environment; namespace modelchecker { namespace multiobjective { +/*! + * Represents the LP Encoding for achievability under simple strategies + * @see http://doi.org/10.18154/RWTH-2023-09669 Chapter 8.3 and 8.4 + */ template class DeterministicSchedsLpChecker { public: @@ -32,9 +37,17 @@ class DeterministicSchedsLpChecker { /*! * Optimizes in the currently given direction - * @return some optimal point found in that direction. + * + * If there is a point in the given area, an achievable point p and a value v are returned such that + * p*w≈q*w≈v, where w is the weight vector and q is an (unknown) optimal point in the overapproximation. + * - if eps is present, the requirement is p*w ≤ q*w ≤ v and |p*w-v| ≤ |eps*w|, and + * - if eps is not present (i.e. empty), the default precision of the underlying LP solver is used. + * If there is no achievable point in the given area, the returned object remains uninitialized. + * + * @param eps if not empty, defines the precision requirement + * */ - boost::optional check(storm::Environment const& env, Polytope overapproximation); + std::optional> check(storm::Environment const& env, Polytope overapproximation, Point const& eps = {}); /*! * Optimizes in the currently given direction, recursively checks for points in the given area. @@ -53,7 +66,6 @@ class DeterministicSchedsLpChecker { private: void initialize(Environment const& env); - bool processEndComponents(std::vector>& ecVars); void initializeLpModel(Environment const& env); // Builds the induced markov chain of the current model and checks whether the resulting value coincide with the result of the lp solver. diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp index 6593c2a3f2..535f534b4a 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp @@ -1,39 +1,36 @@ #include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" -#include "storm/logic/CloneVisitor.h" #include "storm/logic/Formulas.h" -#include "storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h" -#include "storm/modelchecker/prctl/SparseMdpPrctlModelChecker.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h" #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" -#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/models/sparse/MarkovAutomaton.h" #include "storm/models/sparse/Mdp.h" #include "storm/models/sparse/StandardRewardModel.h" -#include "storm/solver/GurobiLpSolver.h" -#include "storm/solver/LpSolver.h" -#include "storm/solver/SmtSolver.h" +#include "storm/solver/LinearEquationSolver.h" +#include "storm/solver/MinMaxLinearEquationSolver.h" #include "storm/storage/BitVector.h" #include "storm/storage/MaximalEndComponentDecomposition.h" -#include "storm/storage/expressions/Expressions.h" +#include "storm/storage/StronglyConnectedComponentDecomposition.h" #include "storm/transformer/EndComponentEliminator.h" +#include "storm/utility/Extremum.h" #include "storm/utility/FilteredRewardModel.h" #include "storm/utility/graph.h" -#include "storm/utility/solver.h" #include "storm/utility/vector.h" +#include "storm/exceptions/UncheckedRequirementException.h" #include "storm/exceptions/UnexpectedException.h" -namespace storm { -namespace modelchecker { -namespace multiobjective { + +namespace storm::modelchecker::multiobjective { template DeterministicSchedsObjectiveHelper::DeterministicSchedsObjectiveHelper(ModelType const& model, storm::modelchecker::multiobjective::Objective const& objective) : model(model), objective(objective) { - // Intentionally left empty + initialize(); } template @@ -47,215 +44,227 @@ storm::storage::BitVector evaluatePropositionalFormula(ModelType const& model, s template std::vector getTotalRewardVector(storm::models::sparse::MarkovAutomaton const& model, storm::logic::Formula const& formula) { - boost::optional rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName(); - typename storm::models::sparse::MarkovAutomaton::RewardModelType const& rewardModel = - rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel(); - - // Get a reward model where the state rewards are scaled accordingly - std::vector stateRewardWeights(model.getNumberOfStates(), storm::utility::zero()); - for (auto const markovianState : model.getMarkovianStates()) { - stateRewardWeights[markovianState] = storm::utility::one() / model.getExitRate(markovianState); + if (formula.isRewardOperatorFormula()) { + boost::optional rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName(); + typename storm::models::sparse::MarkovAutomaton::RewardModelType const& rewardModel = + rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel(); + + // Get a reward model where the state rewards are scaled accordingly + std::vector stateRewardWeights(model.getNumberOfStates(), storm::utility::zero()); + for (auto const markovianState : model.getMarkovianStates()) { + stateRewardWeights[markovianState] = storm::utility::one() / model.getExitRate(markovianState); + } + return rewardModel.getTotalActionRewardVector(model.getTransitionMatrix(), stateRewardWeights); + } else { + assert(formula.isTimeOperatorFormula()); + std::vector result(model.getNumberOfChoices(), storm::utility::zero()); + for (auto const markovianState : model.getMarkovianStates()) { + for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(markovianState)) { + result[choice] = storm::utility::one() / model.getExitRate(markovianState); + } + } + return result; } - return rewardModel.getTotalActionRewardVector(model.getTransitionMatrix(), stateRewardWeights); } template std::vector getTotalRewardVector(storm::models::sparse::Mdp const& model, storm::logic::Formula const& formula) { - boost::optional rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName(); - typename storm::models::sparse::Mdp::RewardModelType const& rewardModel = - rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel(); - return rewardModel.getTotalRewardVector(model.getTransitionMatrix()); + if (formula.isRewardOperatorFormula()) { + boost::optional rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName(); + typename storm::models::sparse::Mdp::RewardModelType const& rewardModel = + rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel(); + return rewardModel.getTotalRewardVector(model.getTransitionMatrix()); + } else { + assert(formula.isTimeOperatorFormula()); + return std::vector(model.getNumberOfChoices(), storm::utility::one()); + } } 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 prob1State : prob1States) { - result[prob1State] = storm::utility::one(); - } - } - { - storm::storage::BitVector prob0States = storm::utility::graph::performProb0A(backwardTransitions, phiStates, psiStates); - for (auto prob0State : prob0States) { - result[prob0State] = storm::utility::zero(); +void DeterministicSchedsObjectiveHelper::initialize() { + STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "Expected a single initial state."); + uint64_t initialState = *model.getInitialStates().begin(); + relevantZeroRewardChoices = storm::storage::BitVector(model.getTransitionMatrix().getRowCount(), false); + auto const& formula = *objective.formula; + 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(); + auto prob1States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), backwardTransitions, + phiStates, psiStates); + auto prob0States = storm::utility::graph::performProb0A(backwardTransitions, phiStates, psiStates); + if (prob0States.get(initialState)) { + constantInitialStateValue = storm::utility::zero(); + } else if (prob1States.get(initialState)) { + constantInitialStateValue = storm::utility::one(); + } + maybeStates = ~(prob0States | prob1States); + // Cut away those states that are not reachable, or only reachable via non-maybe states. + maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates); + for (auto const& state : maybeStates) { + for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) { + auto rowSum = model.getTransitionMatrix().getConstrainedRowSum(choice, prob1States); + if (storm::utility::isZero(rowSum)) { + relevantZeroRewardChoices.set(choice); + } else { + choiceRewards.emplace(choice, rowSum); } } - } 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 rew0State : rew0States) { - result[rew0State] = storm::utility::zero(); - } - } else if (formula.isRewardOperatorFormula() && formula.getSubformula().isTotalRewardFormula()) { + } + } 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().asTotalRewardFormula()); + storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), formula.getSubformula().asEventuallyFormula()); 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 rew0State : rew0States) { - result[rew0State] = storm::utility::zero(); - } - } else { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported."); + rew0States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), + model.getBackwardTransitions(), statesWithoutReward, rew0States); } - 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 = getTotalRewardVector(model, *objective.formula); - - // 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]; - } - } + if (rew0States.get(initialState)) { + constantInitialStateValue = storm::utility::zero(); + } + maybeStates = ~rew0States; + // Cut away those states that are not reachable, or only reachable via non-maybe states. + maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates); + std::vector choiceBasedRewards = getTotalRewardVector(model, *objective.formula); + for (auto const& state : maybeStates) { + for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) { + auto const& value = choiceBasedRewards[choice]; + if (storm::utility::isZero(value)) { + relevantZeroRewardChoices.set(choice); + } else { + choiceRewards.emplace(choice, value); } } - } 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(); - } - // 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 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(); + if (rew0States.get(initialState)) { + constantInitialStateValue = storm::utility::zero(); + } + maybeStates = ~rew0States; + // Cut away those states that are not reachable, or only reachable via non-maybe states. + maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates); + std::vector choiceBasedRewards = getTotalRewardVector(model, *objective.formula); + for (auto const& state : maybeStates) { + for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) { + auto const& value = choiceBasedRewards[choice]; + if (storm::utility::isZero(value)) { + relevantZeroRewardChoices.set(choice); + } else { + choiceRewards.emplace(choice, value); } } + } + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported."); + } + + // negate choice rewards for minimizing objectives + if (storm::solver::minimize(formula.getOptimalityType())) { + for (auto& entry : choiceRewards) { + entry.second *= -storm::utility::one(); + } + } + + // Find out if we have to deal with infinite positive or negative rewards + storm::storage::BitVector negativeRewardChoices(model.getNumberOfChoices(), false); + storm::storage::BitVector positiveRewardChoices(model.getNumberOfChoices(), false); + for (auto const& rew : getChoiceRewards()) { + if (rew.second > storm::utility::zero()) { + positiveRewardChoices.set(rew.first, true); } else { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported."); + assert(rew.second < storm::utility::zero()); + negativeRewardChoices.set(rew.first, true); } - choiceValueOffsets = std::move(result); } - return choiceValueOffsets.get(); + auto backwardTransitions = model.getBackwardTransitions(); + bool hasNegativeEC = + storm::utility::graph::checkIfECWithChoiceExists(model.getTransitionMatrix(), backwardTransitions, getMaybeStates(), negativeRewardChoices); + bool hasPositiveEc = + storm::utility::graph::checkIfECWithChoiceExists(model.getTransitionMatrix(), backwardTransitions, getMaybeStates(), positiveRewardChoices); + STORM_LOG_THROW(!(hasNegativeEC && hasPositiveEc), storm::exceptions::NotSupportedException, + "Objective is not convergent: Infinite positive and infinite negative reward is possible."); + infinityCase = hasNegativeEC ? InfinityCase::HasNegativeInfinite : (hasPositiveEc ? InfinityCase::HasPositiveInfinite : InfinityCase::AlwaysFinite); + + if (infinityCase == InfinityCase::HasNegativeInfinite) { + storm::storage::BitVector negativeEcStates(maybeStates.size(), false); + storm::storage::MaximalEndComponentDecomposition mecs(model.getTransitionMatrix(), backwardTransitions, getMaybeStates()); + for (auto const& mec : mecs) { + if (std::any_of(mec.begin(), mec.end(), [&negativeRewardChoices](auto const& sc) { + return std::any_of(sc.second.begin(), sc.second.end(), [&negativeRewardChoices](auto const& c) { return negativeRewardChoices.get(c); }); + })) { + for (auto const& sc : mec) { + negativeEcStates.set(sc.first); + } + } + } + STORM_LOG_ASSERT(!negativeEcStates.empty(), "Expected some negative ec"); + rewMinusInfEStates = storm::utility::graph::performProbGreater0E(backwardTransitions, getMaybeStates(), negativeEcStates); + STORM_LOG_ASSERT(model.getInitialStates().isSubsetOf(rewMinusInfEStates), "Initial state does not reach all maybestates"); + } } -template -std::vector evaluateOperatorFormula(Environment const& env, storm::models::sparse::Mdp const& model, - storm::logic::Formula const& formula) { - storm::modelchecker::SparseMdpPrctlModelChecker> mc(model); - storm::modelchecker::CheckTask task(formula, false); - auto checkResult = mc.check(env, task); - STORM_LOG_THROW(checkResult && checkResult->isExplicitQuantitativeCheckResult(), storm::exceptions::UnexpectedException, - "Unexpected type of check result for subformula " << formula << "."); - return checkResult->template asExplicitQuantitativeCheckResult().getValueVector(); +template +storm::storage::BitVector const& DeterministicSchedsObjectiveHelper::getMaybeStates() const { + return maybeStates; } -template -std::vector evaluateOperatorFormula(Environment const& env, storm::models::sparse::MarkovAutomaton const& model, - storm::logic::Formula const& formula) { - storm::modelchecker::SparseMarkovAutomatonCslModelChecker> mc(model); - storm::modelchecker::CheckTask task(formula, false); - auto checkResult = mc.check(env, task); - STORM_LOG_THROW(checkResult && checkResult->isExplicitQuantitativeCheckResult(), storm::exceptions::UnexpectedException, - "Unexpected type of check result for subformula " << formula << "."); - return checkResult->template asExplicitQuantitativeCheckResult().getValueVector(); +template +storm::storage::BitVector const& DeterministicSchedsObjectiveHelper::getRewMinusInfEStates() const { + STORM_LOG_ASSERT(getInfinityCase() == InfinityCase::HasNegativeInfinite, "Tried to get -inf states, but there are none"); + return rewMinusInfEStates; } template -std::vector computeValueBounds(Environment const& env, bool lowerValueBounds, ModelType const& model, - storm::logic::Formula const& formula) { - // Change the optimization direction in the formula. - auto newFormula = storm::logic::CloneVisitor().clone(formula); - newFormula->asOperatorFormula().setOptimalityType(lowerValueBounds ? storm::solver::OptimizationDirection::Minimize - : storm::solver::OptimizationDirection::Maximize); - - if (std::is_same::value) { - // don't have to worry about precision in exact mode. - return evaluateOperatorFormula(env, model, *newFormula); - } else { - // Create an environment where sound results are enforced - storm::Environment soundEnv(env); - soundEnv.solver().setForceSoundness(true); - auto result = evaluateOperatorFormula(soundEnv, model, *newFormula); +bool DeterministicSchedsObjectiveHelper::hasConstantInitialStateValue() const { + return constantInitialStateValue.has_value(); +} - auto eps = storm::utility::convertNumber(soundEnv.solver().minMax().getPrecision()); - // Add/substract eps to all entries to make up for precision errors - if (lowerValueBounds) { - eps = -eps; - } - for (auto& v : result) { - v += eps; - } - return result; - } +template +typename DeterministicSchedsObjectiveHelper::ValueType DeterministicSchedsObjectiveHelper::getConstantInitialStateValue() const { + assert(hasConstantInitialStateValue()); + return *constantInitialStateValue; } template -typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getUpperValueBoundAtState(Environment const& env, uint64_t state) const { - computeUpperBounds(env); - return upperResultBounds.get()[state]; +std::map const& DeterministicSchedsObjectiveHelper::getChoiceRewards() const { + return choiceRewards; } template -typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getLowerValueBoundAtState(Environment const& env, uint64_t state) const { - computeLowerBounds(env); - return lowerResultBounds.get()[state]; +storm::storage::BitVector const& DeterministicSchedsObjectiveHelper::getRelevantZeroRewardChoices() const { + return relevantZeroRewardChoices; } template -bool DeterministicSchedsObjectiveHelper::minimizing() const { - return storm::solver::minimize(objective.formula->getOptimalityType()); +typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getUpperValueBoundAtState(uint64_t state) const { + STORM_LOG_ASSERT(maybeStates.get(state), "Expected a maybestate."); + STORM_LOG_ASSERT(upperResultBounds.has_value(), "requested upper value bounds but they were not computed."); + return upperResultBounds->at(state); +} + +template +typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getLowerValueBoundAtState(uint64_t state) const { + STORM_LOG_ASSERT(maybeStates.get(state), "Expected a maybestate."); + STORM_LOG_ASSERT(lowerResultBounds.has_value(), "requested lower value bounds but they were not computed."); + return lowerResultBounds->at(state); +} + +template +typename DeterministicSchedsObjectiveHelper::InfinityCase const& DeterministicSchedsObjectiveHelper::getInfinityCase() const { + return infinityCase; } template @@ -263,308 +272,431 @@ bool DeterministicSchedsObjectiveHelper::isTotalRewardObjective() con return objective.formula->isRewardOperatorFormula() && objective.formula->getSubformula().isTotalRewardFormula(); } -/* -template -ValueType getLpathDfs(uint64_t currentState, ValueType currentValue, storm::storage::BitVector& visited, storm::storage::BitVector const& mecStates, -storm::storage::SparseMatrix const& transitions, uint64_t& counter) { - // exhausive dfs - if (visited.get(currentState)) { - if (++counter % 1000000 == 0) { - std::cout << "\t\t checked " << counter << " paths\n"; - } - return currentValue; +template +bool DeterministicSchedsObjectiveHelper::hasThreshold() const { + return objective.formula->hasBound(); +} + +template +typename DeterministicSchedsObjectiveHelper::ValueType DeterministicSchedsObjectiveHelper::getThreshold() const { + STORM_LOG_ASSERT(hasThreshold(), "Trying to get a threshold but there is none"); + STORM_LOG_THROW(!storm::logic::isStrict(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException, + "Objective " << *objective.originalFormula << ": Strict objective thresholds are not supported."); + ValueType threshold = objective.formula->template getThresholdAs(); + if (storm::solver::minimize(objective.formula->getOptimalityType())) { + threshold *= -storm::utility::one(); // negate minimizing thresholds + STORM_LOG_THROW(!storm::logic::isLowerBound(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException, + "Objective " << *objective.originalFormula << ": Minimizing objective should specify an upper bound."); } else { - visited.set(currentState); - ValueType result = storm::utility::one(); - for (auto const& transition : transitions.getRowGroup(currentState)) { - result = std::min(result, getLpathDfs(transition.getColumn(), currentValue * transition.getValue(), visited, mecStates, transitions, -counter)); - } - visited.set(currentState, false); - return result; + STORM_LOG_THROW(storm::logic::isLowerBound(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException, + "Objective " << *objective.originalFormula << ": Maximizing objective should specify a lower bound."); } + return threshold; +} + +template +bool DeterministicSchedsObjectiveHelper::minimizing() const { + return storm::solver::minimize(objective.formula->getOptimalityType()); } -void oneOfManyEncoding(std::vector& operands, storm::solver::SmtSolver& solver) { - assert(!operands.empty()); - uint64_t currentOperand = 0; - while (currentOperand + 2 < operands.size()) { - solver.add(!(operands[currentOperand] && operands[currentOperand + 1])); - auto auxVar = solver.getManager().declareFreshBooleanVariable(true).getExpression(); - solver.add(storm::expressions::iff(auxVar, (operands[currentOperand] || operands[currentOperand + 1]))); - operands.push_back(auxVar); - currentOperand += 2; +template +void setLowerUpperTotalRewardBoundsToSolver(storm::solver::AbstractEquationSolver& solver, storm::storage::SparseMatrix const& matrix, + std::vector const& rewards, std::vector const& exitProbabilities, + std::optional dir, bool reqLower, bool reqUpper) { + if (!reqLower && !reqUpper) { + return; // nothing to be done! } - if (currentOperand + 1 == operands.size()) { - solver.add(operands[currentOperand]); - } else if (currentOperand + 2 == operands.size()) { - solver.add(storm::expressions::xclusiveor(operands[currentOperand], operands[currentOperand + 1])); + STORM_LOG_ASSERT(!rewards.empty(), "empty reward vector,"); + + auto [minIt, maxIt] = std::minmax_element(rewards.begin(), rewards.end()); + bool const hasNegativeValues = *minIt < storm::utility::zero(); + bool const hasPositiveValues = *maxIt > storm::utility::zero(); + + std::optional lowerBound, upperBound; + + // Get 0 as a trivial lower/upper bound if possible + if (!hasNegativeValues) { + lowerBound = storm::utility::zero(); } -} -*/ - -/* -template -ValueType getExpVisitsUpperBoundForMec(storm::storage::BitVector const& mecStates, storm::storage::SparseMatrix const& transitions) { - auto generalSolver = storm::utility::solver::getLpSolver("mecBounds", storm::solver::LpSolverTypeSelection::Gurobi); - auto solver = dynamic_cast*>(generalSolver.get()); - solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); - - auto one = solver->getConstant(storm::utility::one()); - auto zero = solver->getConstant(storm::utility::zero()); - - std::vector> ins(mecStates.size()), outs(mecStates.size()); - std::vector choiceVars(transitions.getRowCount()), ecVars(mecStates.size()); - std::vector choicesDisjunction; - std::vector leavingSum; - for (auto const& mecState : mecStates) { - std::vector choiceIndicatorVars; - ecVars[mecState] = solver->addLowerBoundedContinuousVariable("e" + std::to_string(mecState), storm::utility::zero()).getExpression(); - ins[mecState].push_back(one); - outs[mecState].push_back(ecVars[mecState]); - leavingSum.push_back(ecVars[mecState]); - choiceIndicatorVars.push_back(solver->addBoundedIntegerVariable("c_mec" + std::to_string(mecState), storm::utility::zero(), -storm::utility::one())); solver->addIndicatorConstraint("", choiceIndicatorVars.back().getBaseExpression().asVariableExpression().getVariable(), 0, -ecVars[mecState] <= zero); for (uint64_t choice = transitions.getRowGroupIndices()[mecState]; choice < transitions.getRowGroupIndices()[mecState + 1]; ++choice) -{ choiceVars[choice] = solver->addLowerBoundedContinuousVariable("y" + std::to_string(choice), storm::utility::zero()).getExpression(); - choiceIndicatorVars.push_back(solver->addBoundedIntegerVariable("c" + std::to_string(choice), storm::utility::zero(), -storm::utility::one()).getExpression()); solver->addIndicatorConstraint("", -choiceIndicatorVars.back().getBaseExpression().asVariableExpression().getVariable(), 0, choiceVars[choice] <= zero); - choicesDisjunction.push_back(choiceVars[choice]); - outs[mecState].push_back(choiceVars[choice]); - for (auto const& entry : transitions.getRow(choice)) { - if (!storm::utility::isZero(entry.getValue())) { - if (mecStates.get(entry.getColumn())) { - ins[entry.getColumn()].push_back(choiceVars[choice] * solver->getConstant(entry.getValue())); - } else { - leavingSum.push_back(choiceVars[choice] * solver->getConstant(entry.getValue())); - } - } - } + if (!hasPositiveValues) { + upperBound = storm::utility::zero(); + } + + // Invoke the respective reward bound computers if needed + std::vector tmpRewards; + if (reqUpper && !upperBound.has_value()) { + if (hasNegativeValues) { + tmpRewards.resize(rewards.size()); + storm::utility::vector::applyPointwise(rewards, tmpRewards, [](ValueType const& v) { return std::max(storm::utility::zero(), v); }); + } + if (dir.has_value() && maximize(*dir)) { + solver.setUpperBound( + storm::modelchecker::helper::BaierUpperRewardBoundsComputer(matrix, hasNegativeValues ? tmpRewards : rewards, exitProbabilities) + .computeUpperBound()); + } else { + solver.setUpperBounds( + storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer(matrix, hasNegativeValues ? tmpRewards : rewards, exitProbabilities) + .computeUpperBounds()); } - //oneOfManyEncoding(choicesSum, *solver); - solver->addConstraint("", storm::expressions::sum(choiceIndicatorVars) == one); } - for (auto const& mecState : mecStates) { - STORM_LOG_ASSERT(!ins[mecState].empty(), "empty in set at a state"); - STORM_LOG_ASSERT(!outs[mecState].empty(), "empty out set at a state"); - solver->addConstraint("", storm::expressions::sum(ins[mecState]) == storm::expressions::sum(outs[mecState])); + if (reqLower && !upperBound.has_value()) { + // For lower bounds we actually compute upper bounds for the negated rewards. + // We therefore need tmpRewards in any way. + tmpRewards.resize(rewards.size()); + storm::utility::vector::applyPointwise(rewards, tmpRewards, + [](ValueType const& v) { return std::max(storm::utility::zero(), -v); }); + if (dir.has_value() && minimize(*dir)) { + solver.setLowerBound( + -storm::modelchecker::helper::BaierUpperRewardBoundsComputer(matrix, tmpRewards, exitProbabilities).computeUpperBound()); + } else { + auto lowerBounds = + storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer(matrix, tmpRewards, exitProbabilities).computeUpperBounds(); + storm::utility::vector::applyPointwise(lowerBounds, lowerBounds, [](ValueType const& v) { return -v; }); + solver.setLowerBounds(std::move(lowerBounds)); + } } - STORM_LOG_ASSERT(!leavingSum.empty(), "empty leaving sum at a mec"); - solver->addConstraint("", storm::expressions::sum(leavingSum) == solver->getConstant(storm::utility::convertNumber(mecStates.getNumberOfSetBits()))); choicesDisjunction.push_back(one); auto boundVar = solver->addUnboundedContinuousVariable("bound", -storm::utility::one()); solver->addGeneralConstraint("", boundVar, storm::solver::GurobiLpSolver::GeneralConstraintOperator::Max, -choicesDisjunction); solver->optimize(); STORM_LOG_THROW(!solver->isInfeasible(), storm::exceptions::UnexpectedException, "MEC LP has infeasable solution"); - STORM_LOG_THROW(!solver->isUnbounded(), storm::exceptions::UnexpectedException, "MEC LP has unbounded solution"); - return solver->getObjectiveValue(); } -*/ -template -std::vector DeterministicSchedsObjectiveHelper::computeUpperBoundOnExpectedVisitingTimes( - storm::storage::SparseMatrix const& modelTransitions, storm::storage::BitVector const& bottomStates, - storm::storage::BitVector const& nonBottomStates, bool hasEndComponents) { - storm::storage::SparseMatrix transitions; - std::vector probabilitiesToBottomStates; - boost::optional> modelToSubsystemStateMapping; - if (hasEndComponents) { - // We assume that end components will always be left (or form a sink state). - // The approach is to give a lower bound lpath on a path that leaves the end component. - // Then we use end component elimination and add a self loop on the 'ec' states with probability 1-lpath - storm::storage::MaximalEndComponentDecomposition mecs(modelTransitions, modelTransitions.transpose(true), nonBottomStates); - auto mecElimination = storm::transformer::EndComponentEliminator::transform(modelTransitions, mecs, nonBottomStates, nonBottomStates); - - probabilitiesToBottomStates.reserve(mecElimination.matrix.getRowCount()); - for (uint64_t row = 0; row < mecElimination.matrix.getRowCount(); ++row) { - if (mecElimination.matrix.getRow(row).getNumberOfEntries() == 0) { - probabilitiesToBottomStates.push_back(storm::utility::one()); - } else { - probabilitiesToBottomStates.push_back(modelTransitions.getConstrainedRowSum(mecElimination.newToOldRowMapping[row], bottomStates)); +template +std::vector computeValuesOfReducedSystem(Environment const& env, storm::storage::SparseMatrix const& submatrix, + std::vector const& exitProbs, std::vector const& rewards, + storm::OptimizationDirection const& dir) { + auto minMaxSolver = storm::solver::GeneralMinMaxLinearEquationSolverFactory().create(env, submatrix); + minMaxSolver->setHasUniqueSolution(true); + minMaxSolver->setOptimizationDirection(dir); + auto req = minMaxSolver->getRequirements(env, dir); + setLowerUpperTotalRewardBoundsToSolver(*minMaxSolver, submatrix, rewards, exitProbs, dir, req.lowerBounds(), req.upperBounds()); + req.clearBounds(); + if (req.validInitialScheduler()) { + std::vector initSched(submatrix.getRowGroupCount()); + auto backwardsTransitions = submatrix.transpose(true); + storm::storage::BitVector statesWithChoice(initSched.size(), false); + storm::storage::BitVector foundStates(initSched.size(), false); + std::vector stack; + for (uint64_t state = 0; state < initSched.size(); ++state) { + if (foundStates.get(state)) { + continue; } - } - - modelToSubsystemStateMapping = std::move(mecElimination.oldToNewStateMapping); - transitions = storm::storage::SparseMatrix(mecElimination.matrix, true); // insert diagonal entries - - // slow down mec states by adding self loop probability 1-lpath - for (auto const& mec : mecs) { - ValueType lpath = storm::utility::one(); - /* //smt/lp - storm::storage::BitVector mecStates(modelTransitions.getRowGroupCount(), false); - uint64_t numTransitions = 0; - uint64_t numChoices = 0; - for (auto const& stateChoices : mec) { - mecStates.set(stateChoices.first); - numTransitions += modelTransitions.getRowGroup(stateChoices.first).getNumberOfEntries(); - numChoices += modelTransitions.getRowGroupSize(stateChoices.first); - } - std::cout << "Checking a mec with " << mecStates.getNumberOfSetBits() << " states " << numChoices << " choices and " << numTransitions << " - transitions.\n"; lpath = storm::utility::one() / getExpVisitsUpperBoundForMec(mecStates, modelTransitions); - }*/ - // We multiply the smallest transition probabilities occurring at each state and MEC-Choice - // as well as the smallest 'exit' probability - ValueType minExitProbability = storm::utility::one(); - for (auto const& stateChoices : mec) { - auto state = stateChoices.first; - ValueType minProb = storm::utility::one(); - for (uint64_t choice = modelTransitions.getRowGroupIndices()[state]; choice < modelTransitions.getRowGroupIndices()[state + 1]; ++choice) { - if (stateChoices.second.count(choice) == 0) { - // The choice leaves the EC, so we take the sum over the exiting probabilities - ValueType exitProbabilitySum = storm::utility::zero(); - for (auto const& transition : modelTransitions.getRow(choice)) { - if (!mec.containsState(transition.getColumn())) { - exitProbabilitySum += transition.getValue(); - } - } - minExitProbability = std::min(minExitProbability, exitProbabilitySum); - } else { - // Get the minimum over all transition probabilities - for (auto const& transition : modelTransitions.getRow(choice)) { - if (!storm::utility::isZero(transition.getValue())) { - minProb = std::min(minProb, transition.getValue()); - } + for (auto choice : submatrix.getRowGroupIndices(state)) { + if (!storm::utility::isZero(exitProbs[choice])) { + initSched[state] = choice - submatrix.getRowGroupIndices()[state]; + statesWithChoice.set(state); + foundStates.set(state); + for (auto const& predecessor : backwardsTransitions.getRow(state)) { + if (!foundStates.get(predecessor.getColumn())) { + stack.push_back(predecessor.getColumn()); + foundStates.set(predecessor.getColumn()); } } + break; } - lpath *= minProb; } - lpath *= minExitProbability; - /* other ideas... - std::vector leastPathProbs(modelTransitions.getRowGroupCount(), storm::utility::one()); - std::vector prevLeastPathProbs(modelTransitions.getRowGroupCount(), storm::utility::one()); - uint64_t mecSize = std::distance(mec.begin(), mec.end()); - for (uint64_t i = 0; i < mecSize; ++i) { - for (auto const& stateChoices : mec) { - auto state = stateChoices.first; - auto currVal = prevLeastPathProbs[state]; - for (auto const& transition : modelTransitions.getRowGroup(state)) { - if (!storm::utility::isZero(transition.getValue()) && transition.getColumn() != state) { - currVal = std::min(currVal, transition.getValue() * prevLeastPathProbs[transition.getColumn()]); - } - } - leastPathProbs[state] = currVal; - } - leastPathProbs.swap(prevLeastPathProbs); + } + while (!stack.empty()) { + auto state = stack.back(); + stack.pop_back(); + for (auto choice : submatrix.getRowGroupIndices(state)) { + auto row = submatrix.getRow(choice); + if (std::any_of(row.begin(), row.end(), [&statesWithChoice](auto const& entry) { + return !storm::utility::isZero(entry.getValue()) && statesWithChoice.get(entry.getColumn()); + })) { + initSched[state] = choice - submatrix.getRowGroupIndices()[state]; + statesWithChoice.set(state); + break; } - lpath = *std::min_element(leastPathProbs.begin(), leastPathProbs.end()); - //lpath = std::max(lpath, storm::utility::convertNumber(1e-2)); // TODO } - if (false) { - storm::storage::BitVector mecStates(modelTransitions.getRowGroupCount(), false); - uint64_t numTransitions = 0; - uint64_t numChoices = 0; - for (auto const& stateChoices : mec) { - mecStates.set(stateChoices.first); - numTransitions += modelTransitions.getRowGroup(stateChoices.first).getNumberOfEntries(); - numChoices += modelTransitions.getRowGroupSize(stateChoices.first); - } - std::cout << "Checking a mec with " << mecStates.getNumberOfSetBits() << " states " << numChoices << " choices and " << numTransitions << " - transitions.\n"; lpath = storm::utility::one(); for (auto const& stateChoices : mec) { storm::storage::BitVector visited = - ~mecStates; uint64_t counter = 0; ValueType lpathOld = lpath; lpath = std::min(lpath, getLpathDfs(stateChoices.first, - storm::utility::one(), visited, mecStates, modelTransitions, counter)); if (lpathOld > lpath) { std::cout << "\tnew lpath is " << - storm::utility::convertNumber(lpath) << ". checked " << counter << " paths.\n"; - } - } - }*/ - STORM_LOG_ASSERT(!storm::utility::isZero(lpath), "unexpected value of lpath"); - STORM_LOG_WARN_COND( - lpath >= storm::utility::convertNumber(0.001), - "Small lower bound for the probability to leave a mec: " << storm::utility::convertNumber(lpath) << ". Numerical issues might occur."); - uint64_t mecState = modelToSubsystemStateMapping.get()[mec.begin()->first]; - // scale all the probabilities at this state with lpath - for (uint64_t mecChoice = transitions.getRowGroupIndices()[mecState]; mecChoice < transitions.getRowGroupIndices()[mecState + 1]; ++mecChoice) { - for (auto& entry : transitions.getRow(mecChoice)) { - if (entry.getColumn() == mecState) { - entry.setValue(entry.getValue() * lpath + storm::utility::one() - lpath); - } else { - entry.setValue(entry.getValue() * lpath); - } + assert(statesWithChoice.get(state)); + for (auto const& predecessor : backwardsTransitions.getRow(state)) { + if (!foundStates.get(predecessor.getColumn())) { + stack.push_back(predecessor.getColumn()); + foundStates.set(predecessor.getColumn()); } - probabilitiesToBottomStates[mecChoice] *= lpath; } } - } else { - transitions = modelTransitions.getSubmatrix(true, nonBottomStates, nonBottomStates); - probabilitiesToBottomStates = modelTransitions.getConstrainedRowGroupSumVector(nonBottomStates, bottomStates); + assert(statesWithChoice.full()); + minMaxSolver->setInitialScheduler(std::move(initSched)); + req.clearValidInitialScheduler(); } + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, + "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); + minMaxSolver->setRequirementsChecked(true); - auto subsystemBounds = storm::modelchecker::helper::BaierUpperRewardBoundsComputer::computeUpperBoundOnExpectedVisitingTimes( - transitions, probabilitiesToBottomStates); - uint64_t subsystemState = 0; + std::vector result(submatrix.getRowGroupCount()); + minMaxSolver->solveEquations(env, result, rewards); - std::vector visitingTimesUpperBounds; - visitingTimesUpperBounds.reserve(bottomStates.size()); - for (uint64_t state = 0; state < bottomStates.size(); ++state) { - if (bottomStates.get(state)) { - visitingTimesUpperBounds.push_back(storm::utility::zero()); + return result; +} + +template +void plusMinMaxSolverPrecision(Environment const& env, ValueType& value) { + if (!storm::NumberTraits::IsExact) { + auto eps = storm::utility::convertNumber(env.solver().minMax().getPrecision()); + if (env.solver().minMax().getRelativeTerminationCriterion()) { + value += value * eps; } else { - if (modelToSubsystemStateMapping) { - visitingTimesUpperBounds.push_back(subsystemBounds[modelToSubsystemStateMapping.get()[state]]); - } else { - visitingTimesUpperBounds.push_back(subsystemBounds[subsystemState]); - } - ++subsystemState; + value += eps; } } - return visitingTimesUpperBounds; } -template -typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getLargestUpperBound(Environment const& env) const { - computeUpperBounds(env); - return *std::max_element(upperResultBounds->begin(), upperResultBounds->end()); +template +void minusMinMaxSolverPrecision(Environment const& env, ValueType& value) { + if (!storm::NumberTraits::IsExact) { + auto eps = storm::utility::convertNumber(env.solver().minMax().getPrecision()); + if (env.solver().minMax().getRelativeTerminationCriterion()) { + value -= value * eps; + } else { + value -= eps; + } + } } -template -void DeterministicSchedsObjectiveHelper::computeUpperBounds(Environment const& env) const { - if (!upperResultBounds) { - upperResultBounds = computeValueBounds(env, false, model, *objective.formula); - auto upperResultBound = objective.upperResultBound; - if (!upperResultBound.is_initialized() && storm::utility::vector::hasInfinityEntry(upperResultBounds.get())) { - STORM_LOG_THROW(objective.formula->isRewardOperatorFormula(), storm::exceptions::NotSupportedException, - "The upper bound for objective " << *objective.originalFormula - << " is infinity at some state. This is only supported for reachability rewards / total rewards."); - STORM_LOG_THROW(objective.formula->getSubformula().isTotalRewardFormula() || objective.formula->getSubformula().isEventuallyFormula(), - storm::exceptions::NotSupportedException, - "The upper bound for objective " << *objective.originalFormula - << " is infinity at some state. This is only supported for reachability rewards / total rewards."); - auto rewards = getTotalRewardVector(model, *objective.formula); - auto zeroValuedStates = storm::utility::vector::filterZero(upperResultBounds.get()); - auto expVisits = computeUpperBoundOnExpectedVisitingTimes(model.getTransitionMatrix(), zeroValuedStates, ~zeroValuedStates, true); - ValueType upperBound = storm::utility::zero(); - for (uint64_t state = 0; state < expVisits.size(); ++state) { - ValueType maxReward = storm::utility::zero(); - for (auto row = model.getTransitionMatrix().getRowGroupIndices()[state], endRow = model.getTransitionMatrix().getRowGroupIndices()[state + 1]; - row < endRow; ++row) { - maxReward = std::max(maxReward, rewards[row]); - } - upperBound += expVisits[state] * maxReward; +template +ValueType sumOfMecRewards(storm::storage::MaximalEndComponent const& mec, std::vector const& rewards) { + auto sum = storm::utility::zero(); + for (auto const& stateChoices : mec) { + storm::utility::Extremum optimalStateValue; + for (auto const& choice : stateChoices.second) { + optimalStateValue &= rewards.at(choice); + } + sum += *optimalStateValue; + } + return sum; +} + +template +ValueType getLowerBoundForNonZeroReachProb(storm::storage::SparseMatrix const& transitions, uint64_t init, uint64_t target) { + storm::storage::BitVector allStates(transitions.getRowGroupCount(), true); + storm::storage::BitVector initialAsBitVector(transitions.getRowGroupCount(), false); + initialAsBitVector.set(init, true); + storm::storage::BitVector targetAsBitVector(transitions.getRowGroupCount(), false); + targetAsBitVector.set(target, true); + auto relevantStates = storm::utility::graph::getReachableStates(transitions, initialAsBitVector, allStates, targetAsBitVector); + auto product = storm::utility::one(); + for (auto state : relevantStates) { + if (state == target) { + continue; + } + storm::utility::Minimum minProb; + for (auto const& entry : transitions.getRowGroup(state)) { + if (relevantStates.get(entry.getColumn())) { + minProb &= entry.getValue(); } - upperResultBound = upperBound; } - storm::utility::vector::clip(upperResultBounds.get(), objective.lowerResultBound, upperResultBound); + product *= *minProb; } + return product; } template -void DeterministicSchedsObjectiveHelper::computeLowerBounds(Environment const& env) const { - if (!lowerResultBounds) { - lowerResultBounds = computeValueBounds(env, true, model, *objective.formula); - storm::utility::vector::clip(lowerResultBounds.get(), objective.lowerResultBound, objective.upperResultBound); - STORM_LOG_THROW(!storm::utility::vector::hasInfinityEntry(lowerResultBounds.get()), storm::exceptions::NotSupportedException, - "The lower bound for objective " << *objective.originalFormula << " is infinity at some state. This is not supported."); +void DeterministicSchedsObjectiveHelper::computeLowerUpperBounds(Environment const& env) const { + assert(!upperResultBounds.has_value() && !lowerResultBounds.has_value()); + auto backwardTransitions = model.getBackwardTransitions(); + auto nonMaybeStates = ~maybeStates; + // Eliminate problematic mecs + storm::storage::MaximalEndComponentDecomposition problMecs(model.getTransitionMatrix(), backwardTransitions, maybeStates, + getRelevantZeroRewardChoices()); + auto quotient1 = storm::transformer::EndComponentEliminator::transform(model.getTransitionMatrix(), problMecs, maybeStates, maybeStates); + auto backwardTransitions1 = quotient1.matrix.transpose(true); + std::vector rewards1(quotient1.matrix.getRowCount(), storm::utility::zero()); + std::vector exitProbs1(quotient1.matrix.getRowCount(), storm::utility::zero()); + storm::storage::BitVector subsystemChoices1(quotient1.matrix.getRowCount(), true); + storm::storage::BitVector allQuotient1States(quotient1.matrix.getRowGroupCount(), true); + for (uint64_t choice = 0; choice < quotient1.matrix.getRowCount(); ++choice) { + auto const& oldChoice = quotient1.newToOldRowMapping[choice]; + if (auto findRes = choiceRewards.find(oldChoice); findRes != choiceRewards.end()) { + rewards1[choice] = findRes->second; + } + if (quotient1.sinkRows.get(choice)) { + subsystemChoices1.set(choice, false); + exitProbs1[choice] = storm::utility::one(); + } else if (auto prob = model.getTransitionMatrix().getConstrainedRowSum(oldChoice, nonMaybeStates); !storm::utility::isZero(prob)) { + exitProbs1[choice] = prob; + subsystemChoices1.set(choice, false); + } + } + // Assert that every state can eventually exit the subsystem. If that is not the case, we get infinite reward as problematic (aka zero mecs) are removed. + auto exitStates1 = quotient1.matrix.getRowGroupFilter(~subsystemChoices1, false); + STORM_LOG_THROW(storm::utility::graph::performProbGreater0E(backwardTransitions1, allQuotient1States, exitStates1).full(), + storm::exceptions::InvalidOperationException, "Objective " << *objective.originalFormula << " does not induce finite value."); + + // Compute bounds for finite cases + if (getInfinityCase() != InfinityCase::HasNegativeInfinite) { + auto result1 = computeValuesOfReducedSystem(env, quotient1.matrix, exitProbs1, rewards1, storm::OptimizationDirection::Minimize); + lowerResultBounds = std::vector(model.getNumberOfStates(), storm::utility::zero()); + for (auto const& state : maybeStates) { + ValueType val = result1.at(quotient1.oldToNewStateMapping.at(state)); + minusMinMaxSolverPrecision(env, val); + (*lowerResultBounds)[state] = val; + } + } + if (getInfinityCase() != InfinityCase::HasPositiveInfinite) { + auto result1 = computeValuesOfReducedSystem(env, quotient1.matrix, exitProbs1, rewards1, storm::OptimizationDirection::Maximize); + upperResultBounds = std::vector(model.getNumberOfStates(), storm::utility::zero()); + for (auto const& state : maybeStates) { + ValueType val = result1.at(quotient1.oldToNewStateMapping.at(state)); + plusMinMaxSolverPrecision(env, val); + if (infinityCase == InfinityCase::HasNegativeInfinite) { // Upper bound has to be at least 0 to allow for "trivial solution" in encoding + val = std::max(val, storm::utility::zero()); + } + (*upperResultBounds)[state] = val; + } + } + if (getInfinityCase() != InfinityCase::AlwaysFinite) { + assert(getInfinityCase() == InfinityCase::HasPositiveInfinite || getInfinityCase() == InfinityCase::HasNegativeInfinite); + STORM_LOG_THROW( + hasThreshold() || getInfinityCase() == InfinityCase::HasNegativeInfinite, storm::exceptions::NotSupportedException, + "The upper bound for objective " << *objective.originalFormula << " is infinity at some state. This is only supported for thresholded objectives"); + // Eliminate remaining mecs + storm::storage::MaximalEndComponentDecomposition remainingMecs(quotient1.matrix, backwardTransitions1, allQuotient1States, + subsystemChoices1); + STORM_LOG_ASSERT(!remainingMecs.empty(), "Incorrect infinityCase: There is no MEC with rewards."); + auto quotient2 = + storm::transformer::EndComponentEliminator::transform(quotient1.matrix, remainingMecs, allQuotient1States, allQuotient1States); + std::vector rewards2(quotient2.matrix.getRowCount(), storm::utility::zero()); + std::vector exitProbs2(quotient2.matrix.getRowCount(), storm::utility::zero()); + for (uint64_t choice = 0; choice < quotient2.matrix.getRowCount(); ++choice) { + auto const& oldChoice = quotient2.newToOldRowMapping[choice]; + if (quotient2.sinkRows.get(choice)) { + exitProbs2[choice] = storm::utility::one(); + } else { + rewards2[choice] = rewards1[oldChoice]; + exitProbs2[choice] = exitProbs1[oldChoice]; + } + } + + // Add rewards within ECs + auto initialState2 = quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(*model.getInitialStates().begin())); + ValueType rewardValueForPosInfCase; + if (getInfinityCase() == InfinityCase::HasPositiveInfinite) { + rewardValueForPosInfCase = getThreshold(); + ; + // We need to substract a lower bound for the value at the initial state + std::vector rewards2Negative; + rewards2Negative.reserve(rewards2.size()); + for (auto const& rew : rewards2) { + rewards2Negative.push_back(std::min(rew, storm::utility::zero())); + } + auto lower2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2Negative, storm::OptimizationDirection::Minimize); + rewardValueForPosInfCase -= lower2.at(initialState2); + } + for (auto const& mec : remainingMecs) { + auto mecState2 = quotient2.oldToNewStateMapping[mec.begin()->first]; + ValueType mecReward = getInfinityCase() == InfinityCase::HasPositiveInfinite + ? sumOfMecRewards(mec, rewards1) + : sumOfMecRewards(mec, rewards1); + mecReward /= VisitingTimesHelper::computeMecTraversalLowerBound(mec, quotient1.matrix); + auto groupIndices = quotient2.matrix.getRowGroupIndices(mecState2); + for (auto const choice2 : groupIndices) { + rewards2[choice2] += mecReward; + } + if (getInfinityCase() == InfinityCase::HasPositiveInfinite) { + auto sinkChoice = quotient2.sinkRows.getNextSetIndex(*groupIndices.begin()); + STORM_LOG_ASSERT(sinkChoice < quotient2.matrix.getRowGroupIndices()[mecState2 + 1], "EC state in quotient has no sink row."); + rewards2[sinkChoice] += rewardValueForPosInfCase / getLowerBoundForNonZeroReachProb(quotient2.matrix, initialState2, mecState2); + } + } + + // Compute and insert missing bounds + if (getInfinityCase() == InfinityCase::HasNegativeInfinite) { + auto result2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2, storm::OptimizationDirection::Minimize); + lowerResultBounds = std::vector(model.getNumberOfStates(), storm::utility::zero()); + for (auto const& state : maybeStates) { + ValueType val = result2.at(quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(state))); + minusMinMaxSolverPrecision(env, val); + (*lowerResultBounds)[state] = val; + } + } else { + assert(getInfinityCase() == InfinityCase::HasPositiveInfinite); + auto result2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2, storm::OptimizationDirection::Maximize); + upperResultBounds = std::vector(model.getNumberOfStates(), storm::utility::zero()); + for (auto const& state : maybeStates) { + ValueType val = result2.at(quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(state))); + plusMinMaxSolverPrecision(env, val); + (*upperResultBounds)[state] = val; + } + } } + STORM_LOG_THROW( + std::all_of(maybeStates.begin(), maybeStates.end(), [this](auto const& s) { return this->lowerResultBounds->at(s) <= this->upperResultBounds->at(s); }), + storm::exceptions::UnexpectedException, "Pre-computed lower bound exceeds upper bound."); } template -typename ModelType::ValueType DeterministicSchedsObjectiveHelper::evaluateOnModel(Environment const& env, ModelType const& evaluatedModel) const { - return evaluateOperatorFormula(env, evaluatedModel, *objective.formula)[*evaluatedModel.getInitialStates().begin()]; +typename DeterministicSchedsObjectiveHelper::ValueType DeterministicSchedsObjectiveHelper::evaluateScheduler( + Environment const& env, storm::storage::BitVector const& selectedChoices) const { + STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1u, "Expected a single initial state."); + STORM_LOG_ASSERT(model.getInitialStates().isSubsetOf(getMaybeStates()), "Expected initial state to be maybestate."); + STORM_LOG_ASSERT(selectedChoices.getNumberOfSetBits() == model.getNumberOfStates(), "invalid choice selection."); + storm::storage::BitVector allStates(model.getNumberOfStates(), true); + auto selectedMatrix = model.getTransitionMatrix().getSubmatrix(false, selectedChoices, allStates); + assert(selectedMatrix.getRowCount() == selectedMatrix.getRowGroupCount()); + selectedMatrix.makeRowGroupingTrivial(); + auto subMaybeStates = + getMaybeStates() & storm::utility::graph::getReachableStates(selectedMatrix, model.getInitialStates(), getMaybeStates(), ~getMaybeStates()); + auto bsccCandidates = storm::utility::graph::performProbGreater0(selectedMatrix.transpose(true), subMaybeStates, ~subMaybeStates); + bsccCandidates.complement(); // i.e. states that can not reach non-maybe states + if (!bsccCandidates.empty()) { + storm::storage::StronglyConnectedComponentDecompositionOptions bsccOptions; + bsccOptions.onlyBottomSccs(true).subsystem(bsccCandidates); + storm::storage::StronglyConnectedComponentDecomposition bsccs(selectedMatrix, bsccOptions); + for (auto const& bscc : bsccs) { + for (auto const& s : bscc) { + subMaybeStates.set(s, false); + } + } + } + + if (subMaybeStates.get(*model.getInitialStates().begin())) { + storm::solver::GeneralLinearEquationSolverFactory factory; + bool const useEqSysFormat = factory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + auto eqSysMatrix = selectedMatrix.getSubmatrix(true, subMaybeStates, subMaybeStates, useEqSysFormat); + assert(eqSysMatrix.getRowCount() == eqSysMatrix.getRowGroupCount()); + auto exitProbs = selectedMatrix.getConstrainedRowSumVector(subMaybeStates, ~subMaybeStates); + std::vector rewards; + rewards.reserve(exitProbs.size()); + auto const& allRewards = getChoiceRewards(); + for (auto const& state : getMaybeStates()) { + auto choice = selectedChoices.getNextSetIndex(model.getTransitionMatrix().getRowGroupIndices()[state]); + STORM_LOG_ASSERT(choice < model.getTransitionMatrix().getRowGroupIndices()[state + 1], " no choice selected at state" << state); + if (subMaybeStates.get(state)) { + if (auto findRes = allRewards.find(choice); findRes != allRewards.end()) { + rewards.push_back(findRes->second); + } else { + rewards.push_back(storm::utility::zero()); + } + } else { + STORM_LOG_ASSERT(!bsccCandidates.get(state) || allRewards.count(choice) == 0, "Strategy selected a bscc with rewards"); + } + } + if (useEqSysFormat) { + eqSysMatrix.convertToEquationSystem(); + } + auto solver = factory.create(env, eqSysMatrix); + storm::storage::BitVector init = model.getInitialStates() % subMaybeStates; + assert(init.getNumberOfSetBits() == 1); + auto const initState = *init.begin(); + solver->setRelevantValues(std::move(init)); + auto req = solver->getRequirements(env); + if (!useEqSysFormat) { + setLowerUpperTotalRewardBoundsToSolver(*solver, eqSysMatrix, rewards, exitProbs, std::nullopt, req.lowerBounds(), req.upperBounds()); + req.clearUpperBounds(); + req.clearLowerBounds(); + } + STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, + "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked."); + std::vector x(rewards.size()); + solver->solveEquations(env, x, rewards); + + return x[initState]; + } else { + // initial state is on a bscc + return storm::utility::zero(); + } } template class DeterministicSchedsObjectiveHelper>; template class DeterministicSchedsObjectiveHelper>; template class DeterministicSchedsObjectiveHelper>; template class DeterministicSchedsObjectiveHelper>; -} // namespace multiobjective -} // namespace modelchecker -} // namespace storm +} // namespace storm::modelchecker::multiobjective \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h index fa2524f3af..a232693e92 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h @@ -1,18 +1,21 @@ #pragma once -#include #include +#include #include "storm/modelchecker/multiobjective/Objective.h" #include "storm/storage/BitVector.h" -#include "storm/storage/SparseMatrix.h" namespace storm { +namespace storage { +template +class MaximalEndComponentDecomposition; +} + class Environment; -namespace modelchecker { -namespace multiobjective { +namespace modelchecker::multiobjective { template class DeterministicSchedsObjectiveHelper { @@ -21,46 +24,89 @@ class DeterministicSchedsObjectiveHelper { DeterministicSchedsObjectiveHelper(ModelType const& model, Objective const& objective); /*! - * Returns states and values for states that are independent of the scheduler. + * Returns true iff the value at the (unique) initial state is a constant (that is independent of the scheduler) + */ + bool hasConstantInitialStateValue() const; + + /*! + * If hasConstantInitialStateValue is true, this returns that constant value. + */ + ValueType getConstantInitialStateValue() const; + + /*! + * Returns those states for which the scheduler potentially influences the value. + * + * @return + */ + storm::storage::BitVector const& getMaybeStates() const; + + /*! + * Returns the set of states for which value -inf is possible + * @pre getInfinityCase() == HasNegativeInfinite has to hold + * @return */ - std::map const& getSchedulerIndependentStateValues() const; + storm::storage::BitVector const& getRewMinusInfEStates() 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 + * Returns the choice rewards if non-zero. + * This does not include choices of non-maybe states */ - std::map const& getChoiceValueOffsets() const; + std::map const& getChoiceRewards() const; - ValueType const& getUpperValueBoundAtState(Environment const& env, uint64_t state) const; - ValueType const& getLowerValueBoundAtState(Environment const& env, uint64_t state) const; + /*! + * Returns the choices that (i) originate from a maybestate and (ii) have zero value + * + */ + storm::storage::BitVector const& getRelevantZeroRewardChoices() const; - ValueType const& getLargestUpperBound(Environment const& env) const; + ValueType const& getUpperValueBoundAtState(uint64_t state) const; + ValueType const& getLowerValueBoundAtState(uint64_t state) const; + + enum class InfinityCase { AlwaysFinite, HasPositiveInfinite, HasNegativeInfinite }; + + InfinityCase const& getInfinityCase() const; bool minimizing() const; + void computeLowerUpperBounds(Environment const& env) const; + /*! * Returns true if this is a total reward objective */ bool isTotalRewardObjective() const; - ValueType evaluateOnModel(Environment const& env, ModelType const& evaluatedModel) const; + /*! + * Returns true, if this objective specifies a threshold + */ + bool hasThreshold() const; - static std::vector computeUpperBoundOnExpectedVisitingTimes(storm::storage::SparseMatrix const& modelTransitions, - storm::storage::BitVector const& bottomStates, - storm::storage::BitVector const& nonBottomStates, bool hasEndComponents); + /*! + * Returns the threshold (if specified) + */ + ValueType getThreshold() const; + + /*! + * Computes the value at the initial state under the given scheduler. + * @return + */ + ValueType evaluateScheduler(Environment const& env, storm::storage::BitVector const& selectedChoices) const; private: - void computeUpperBounds(Environment const& env) const; - void computeLowerBounds(Environment const& env) const; + void initialize(); + + storm::storage::BitVector maybeStates; // S_?^j + storm::storage::BitVector rewMinusInfEStates; // S_{-infty}^j + std::optional constantInitialStateValue; + storm::storage::BitVector relevantZeroRewardChoices; + std::map choiceRewards; + InfinityCase infinityCase; - mutable boost::optional> schedulerIndependentStateValues; - mutable boost::optional> choiceValueOffsets; - mutable boost::optional> lowerResultBounds; - mutable boost::optional> upperResultBounds; + // Computed on demand: + mutable std::optional> lowerResultBounds; + mutable std::optional> upperResultBounds; ModelType const& model; Objective const& objective; }; -} // namespace multiobjective -} // namespace modelchecker +} // namespace modelchecker::multiobjective } // namespace storm \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp index ee006f664e..e475736ce6 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.cpp @@ -302,6 +302,7 @@ DeterministicSchedsParetoExplorer::Determini objectiveHelper.reserve(objectives.size()); for (auto const& obj : objectives) { objectiveHelper.emplace_back(*model, obj); + STORM_LOG_ASSERT(!objectiveHelper.back().hasThreshold(), "Unexpected input: got a Pareto query with a thresholded objective."); } lpChecker = std::make_shared>(*model, objectiveHelper); if (preprocessorResult.containsOnlyTotalRewardFormulas()) { @@ -309,54 +310,14 @@ DeterministicSchedsParetoExplorer::Determini } else { wvChecker = nullptr; } - if (storm::settings::getModule().isShowStatisticsSet()) { - storm::utility::Stopwatch sw(true); - std::string modelname = "original-model"; - std::vector models; - models.push_back(&preprocessorResult.originalModel); - models.push_back(model.get()); - for (SparseModelType const* m : models) { - STORM_PRINT_AND_LOG("#STATS " << m->getNumberOfStates() << " states in " << modelname << '\n'); - STORM_PRINT_AND_LOG("#STATS " << m->getNumberOfChoices() << " choices in " << modelname << '\n'); - STORM_PRINT_AND_LOG("#STATS " << m->getNumberOfTransitions() << " transitions in " << modelname << '\n'); - storm::RationalNumber numScheds = storm::utility::one(); - for (uint64_t state = 0; state < m->getNumberOfStates(); ++state) { - storm::RationalNumber numChoices = storm::utility::convertNumber(m->getNumberOfChoices(state)); - numScheds *= storm::utility::max(storm::utility::one(), numChoices); - } - auto numSchedsStr = storm::utility::to_string(numScheds); - STORM_PRINT_AND_LOG("#STATS " << numSchedsStr.front() << "e" << (numSchedsStr.size() - 1) << " memoryless deterministic schedulers in " << modelname - << '\n'); - storm::storage::MaximalEndComponentDecomposition mecs(*m); - uint64_t nonConstMecCounter = 0; - uint64_t nonConstMecStateCounter = 0; - for (auto const& mec : mecs) { - bool mecHasNonConstValue = false; - for (auto const& stateChoices : mec) { - for (auto const& helper : objectiveHelper) { - if (helper.getSchedulerIndependentStateValues().count(stateChoices.first) == 0) { - ++nonConstMecStateCounter; - mecHasNonConstValue = true; - break; - } - } - } - if (mecHasNonConstValue) - ++nonConstMecCounter; - } - STORM_PRINT_AND_LOG("#STATS " << nonConstMecCounter << " non-constant MECS in " << modelname << '\n'); - STORM_PRINT_AND_LOG("#STATS " << nonConstMecStateCounter << " non-constant MEC States in " << modelname << '\n'); - // Print the same statistics for the unfolded model as well. - modelname = "unfolded-model"; - } - sw.stop(); - STORM_PRINT_AND_LOG("#STATS " << sw << " seconds for computing these statistics.\n"); - } } template std::unique_ptr DeterministicSchedsParetoExplorer::check(Environment const& env) { clean(); + for (auto& obj : objectiveHelper) { + obj.computeLowerUpperBounds(env); + } initializeFacets(env); // Compute the relative precision in each dimension. @@ -385,7 +346,8 @@ std::unique_ptr DeterministicSchedsParetoExplorer(1e-8); } } - STORM_PRINT_AND_LOG("Relative precision is " << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(eps)) << '\n'); + STORM_LOG_INFO("Relative precision for deterministic scheduler Pareto explorer is " + << storm::utility::vector::toString(storm::utility::vector::convertNumericVector(eps)) << '\n'); } else { STORM_LOG_THROW(env.modelchecker().multi().getPrecisionType() == MultiObjectiveModelCheckerEnvironment::PrecisionType::Absolute, storm::exceptions::IllegalArgumentException, "Unknown multiobjective precision type."); @@ -407,12 +369,6 @@ std::unique_ptr DeterministicSchedsParetoExplorer(transformObjectiveValuesToOriginal(objectives, p.second.get()))); } } - if (storm::settings::getModule().isShowStatisticsSet()) { - STORM_PRINT_AND_LOG("#STATS " << paretoPoints.size() << " Pareto points\n"); - STORM_PRINT_AND_LOG("#STATS " << unachievableAreas.size() << " unachievable areas\n"); - STORM_PRINT_AND_LOG("#STATS " << overApproximation->getHalfspaces().size() << " unachievable halfspaces\n"); - STORM_PRINT_AND_LOG(lpChecker->getStatistics("#STATS ")); - } return std::make_unique>(originalModelInitialState, std::move(paretoPoints), nullptr, nullptr); } @@ -509,21 +465,25 @@ void DeterministicSchedsParetoExplorer::init for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { std::vector weightVector(objectives.size(), storm::utility::zero()); weightVector[objIndex] = storm::utility::one(); - std::vector point; + std::vector pointCoord; + GeometryValueType offset; if (wvChecker) { wvChecker->check(env, storm::utility::vector::convertNumericVector(weightVector)); - point = storm::utility::vector::convertNumericVector(wvChecker->getUnderApproximationOfInitialStateResults()); - negateMinObjectives(point); + pointCoord = storm::utility::vector::convertNumericVector(wvChecker->getUnderApproximationOfInitialStateResults()); + negateMinObjectives(pointCoord); + auto upperBoundPoint = storm::utility::vector::convertNumericVector(wvChecker->getOverApproximationOfInitialStateResults()); + negateMinObjectives(upperBoundPoint); + offset = storm::utility::vector::dotProduct(weightVector, upperBoundPoint); } else { lpChecker->setCurrentWeightVector(env, weightVector); auto optionalPoint = lpChecker->check(env, overApproximation); - STORM_LOG_THROW(optionalPoint.is_initialized(), storm::exceptions::UnexpectedException, "Unable to find a point in the current overapproximation."); - point = std::move(optionalPoint.get()); + STORM_LOG_THROW(optionalPoint.has_value(), storm::exceptions::UnexpectedException, "Unable to find a point in the current overapproximation."); + pointCoord = std::move(optionalPoint->first); + offset = std::move(optionalPoint->second); } - Point p(point); + Point p(pointCoord); p.setOnFacet(); // Adapt the overapproximation - GeometryValueType offset = storm::utility::vector::dotProduct(weightVector, p.get()); addHalfspaceToOverApproximation(env, weightVector, offset); pointset.addPoint(env, std::move(p)); } @@ -549,13 +509,8 @@ template std::vector DeterministicSchedsParetoExplorer::getReferenceCoordinates(Environment const& env) const { std::vector result; for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { - ModelValueType value; - if (objectiveHelper[objIndex].minimizing()) { - value = -objectiveHelper[objIndex].getUpperValueBoundAtState(env, *model->getInitialStates().begin()); - } else { - value = objectiveHelper[objIndex].getLowerValueBoundAtState(env, *model->getInitialStates().begin()); - } - result.push_back(storm::utility::convertNumber(value)); + result.push_back( + storm::utility::convertNumber(objectiveHelper[objIndex].getLowerValueBoundAtState(*model->getInitialStates().begin()))); } return result; } @@ -605,24 +560,29 @@ template bool DeterministicSchedsParetoExplorer::optimizeAndSplitFacet(Environment const& env, Facet& f) { // Invoke optimization and insert the explored points boost::optional optPointId; - std::vector point; + std::vector pointCoord; + GeometryValueType offset; if (wvChecker) { wvChecker->check(env, storm::utility::vector::convertNumericVector(f.getHalfspace().normalVector())); - point = storm::utility::vector::convertNumericVector(wvChecker->getUnderApproximationOfInitialStateResults()); - negateMinObjectives(point); + pointCoord = storm::utility::vector::convertNumericVector(wvChecker->getUnderApproximationOfInitialStateResults()); + negateMinObjectives(pointCoord); + auto upperBoundPoint = storm::utility::vector::convertNumericVector(wvChecker->getOverApproximationOfInitialStateResults()); + negateMinObjectives(upperBoundPoint); + offset = storm::utility::vector::dotProduct(f.getHalfspace().normalVector(), upperBoundPoint); } else { auto currentArea = overApproximation->intersection(f.getHalfspace().invert()); - auto optionalPoint = lpChecker->check(env, currentArea); - if (optionalPoint.is_initialized()) { - point = std::move(optionalPoint.get()); + auto optionalPoint = lpChecker->check(env, overApproximation, eps); + if (optionalPoint.has_value()) { + pointCoord = std::move(optionalPoint->first); } else { // As we did not find any feasable solution in the given area, we take a point that lies on the facet - point = pointset.getPoint(f.getPoints().front()).get(); + pointCoord = pointset.getPoint(f.getPoints().front()).get(); } + offset = std::move(optionalPoint->second); } - Point p(point); + + Point p(pointCoord); p.setOnFacet(); - GeometryValueType offset = storm::utility::vector::dotProduct(f.getHalfspace().normalVector(), p.get()); addHalfspaceToOverApproximation(env, f.getHalfspace().normalVector(), offset); optPointId = pointset.addPoint(env, std::move(p)); diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h index f878b5dd00..69cc982fda 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h @@ -21,6 +21,10 @@ class Environment; namespace modelchecker { namespace multiobjective { +/*! + * Implements the exploration of the Pareto front. + * @see http://doi.org/10.18154/RWTH-2023-09669 Chapter 8.5 + */ template class DeterministicSchedsParetoExplorer { public: diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.cpp new file mode 100644 index 0000000000..eb70e18266 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.cpp @@ -0,0 +1,151 @@ +#include "storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h" + +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/transformer/EndComponentEliminator.h" +#include "storm/utility/Extremum.h" +#include "storm/utility/constants.h" + +namespace storm::modelchecker::multiobjective { + +template +ValueType VisitingTimesHelper::computeMecTraversalLowerBound(storm::storage::MaximalEndComponent const& mec, + storm::storage::SparseMatrix const& transitions, + bool assumeOptimalTransitionProbabilities) { + STORM_LOG_ASSERT(mec.size() > 0, "empty mec not expected"); + auto res = storm::utility::one(); + if (mec.size() == 1) { + return res; + } + // Whenever s' is reachable from s, there is at least one finite path that visits each state of the mec at most once. + // We compute a naive lower bound for the probability of such a single path using the product of the smallest probabilities occurring at each state + // (Excluding self-loop probabilities) + for (auto const& stateChoices : mec) { + storm::utility::Minimum v; + auto const& s = stateChoices.first; + for (auto const& c : stateChoices.second) { + auto row = transitions.getRow(c); + if (assumeOptimalTransitionProbabilities) { + // Count the number of non-selfloop entries and assume a uniform distribution (which gives us the largest minimal probability) + auto numEntries = row.getNumberOfEntries(); + for (auto const& entry : row) { + if (entry.getColumn() > s) { + if (entry.getColumn() == s) { + --numEntries; + } + break; + } + } + if (numEntries > 0) { + v &= storm::utility::one() / storm::utility::convertNumber(numEntries); + } + } else { + // actually determine the minimal probability + for (auto const& entry : row) { + if (entry.getColumn() != s && !storm::utility::isZero(entry.getValue())) { + v &= entry.getValue(); + } + } + } + } + STORM_LOG_ASSERT(!v.empty(), "self-loop state in non-singleton mec?"); + res *= *v; + } + return res; +} + +template +ValueType VisitingTimesHelper::computeMecVisitsUpperBound(storm::storage::MaximalEndComponent const& mec, + storm::storage::SparseMatrix const& transitions, + bool assumeOptimalTransitionProbabilities) { + auto const one = storm::utility::one(); + auto const traversalLowerBound = computeMecTraversalLowerBound(mec, transitions, assumeOptimalTransitionProbabilities); + if (assumeOptimalTransitionProbabilities) { + // We assume that the probability to go back to the MEC using an exiting choice is zero. + return one / traversalLowerBound; + } else { + // compute the largest probability to go back to the MEC when using an exiting choice + storm::utility::Maximum q(storm::utility::zero()); + for (auto const& stateChoices : mec) { + for (auto c : transitions.getRowGroupIndices(stateChoices.first)) { + if (stateChoices.second.contains(c)) { + continue; // not an exit choice! + } + auto choiceValue = storm::utility::zero(); + for (auto const& entry : transitions.getRow(c)) { + if (mec.containsState(entry.getColumn())) { + choiceValue += entry.getValue(); + } + } + q &= choiceValue; + } + } + return one / (traversalLowerBound * (one - *q)); + } +} + +template +std::vector VisitingTimesHelper::computeUpperBoundsOnExpectedVisitingTimes( + storm::storage::BitVector const& subsystem, storm::storage::SparseMatrix const& transitions, + storm::storage::SparseMatrix const& backwardTransitions) { + storm::storage::MaximalEndComponentDecomposition mecs(transitions, backwardTransitions, subsystem); + storm::storage::BitVector allStates(subsystem.size(), true); + auto quotientData = storm::transformer::EndComponentEliminator::transform(transitions, mecs, subsystem, subsystem); + auto notSubsystem = ~subsystem; + + // map collapsed states in the quotient to the value l that is to be multiplied to the transitions of those states + std::map collapsedStateToValueMap; + for (auto const& mec : mecs) { + collapsedStateToValueMap.emplace(quotientData.oldToNewStateMapping.at(mec.begin()->first), computeMecTraversalLowerBound(mec, transitions)); + } + + // Create a modified quotient with scaled transitions + storm::storage::SparseMatrixBuilder modifiedQuotientBuilder(quotientData.matrix.getRowCount(), quotientData.matrix.getColumnCount(), 0, true, + true, quotientData.matrix.getRowGroupCount()); + std::vector toSinkProbabilities; + toSinkProbabilities.reserve(quotientData.matrix.getRowGroupCount()); + for (uint64_t state = 0; state < quotientData.matrix.getRowGroupCount(); ++state) { + modifiedQuotientBuilder.newRowGroup(quotientData.matrix.getRowGroupIndices()[state]); + ValueType scalingFactor = storm::utility::one(); + if (auto findRes = collapsedStateToValueMap.find(state); findRes != collapsedStateToValueMap.end()) { + scalingFactor = findRes->second; + } + for (auto const choice : quotientData.matrix.getRowGroupIndices(state)) { + if (!storm::utility::isOne(scalingFactor)) { + modifiedQuotientBuilder.addDiagonalEntry(choice, storm::utility::one() - scalingFactor); + } + for (auto const& entry : quotientData.matrix.getRow(choice)) { + modifiedQuotientBuilder.addNextValue(choice, entry.getColumn(), scalingFactor * entry.getValue()); + } + if (quotientData.sinkRows.get(choice)) { + toSinkProbabilities.push_back(scalingFactor); + } else { + toSinkProbabilities.push_back(scalingFactor * transitions.getConstrainedRowSum(quotientData.newToOldRowMapping.at(choice), notSubsystem)); + } + } + } + auto modifiedQuotient = modifiedQuotientBuilder.build(); + STORM_LOG_ASSERT(modifiedQuotient.getRowCount() == toSinkProbabilities.size(), "Unexpected row count"); + + // compute upper visiting time bounds + auto quotientUpperBounds = + storm::modelchecker::helper::BaierUpperRewardBoundsComputer::computeUpperBoundOnExpectedVisitingTimes(modifiedQuotient, toSinkProbabilities); + + std::vector result; + result.reserve(subsystem.size()); + for (uint64_t state = 0; state < subsystem.size(); ++state) { + auto const& quotientState = quotientData.oldToNewStateMapping[state]; + if (quotientState < quotientData.matrix.getRowGroupCount()) { + result.push_back(quotientUpperBounds.at(quotientState)); + } else { + result.push_back(-storm::utility::one()); // not in given subsystem; + } + } + return result; +} + +template class VisitingTimesHelper; +template class VisitingTimesHelper; + +} // namespace storm::modelchecker::multiobjective \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h b/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h new file mode 100644 index 0000000000..9de4353237 --- /dev/null +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/VisitingTimesHelper.h @@ -0,0 +1,48 @@ +#pragma once + +#include + +namespace storm::storage { +class BitVector; +class MaximalEndComponent; +template +class SparseMatrix; +} // namespace storm::storage + +namespace storm::modelchecker::multiobjective { + +/*! + * Provides helper functions for computing bounds on the expected visiting times + * @see http://doi.org/10.18154/RWTH-2023-09669 Chapter 8.2 + */ +template +class VisitingTimesHelper { + public: + /*! + * Computes value l in (0,1] such that for any pair of distinct states s,s' and deterministic, memoryless scheduler under which s' is reachable from s we + * have that l is a lower bound for the probability that we reach s' from s without visiting s in between. + * @param assumeOptimalTransitionProbabilities if true, we assume transitions with the same graph structure (except for selfloops) with uniform distribution + */ + static ValueType computeMecTraversalLowerBound(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix const& transitions, + bool assumeOptimalTransitionProbabilities = false); + + /*! + * Computes an upper bound for the largest finite expected number of times a state s in the given MEC is visited without leaving the MEC (under all + * deterministic, memoryless scheduler) + * @param assumeOptimalTransitionProbabilities if true, we assume transitions with the same graph structure (except for selfloops and MEC exiting choices) + with uniform distribution + + */ + static ValueType computeMecVisitsUpperBound(storm::storage::MaximalEndComponent const& mec, storm::storage::SparseMatrix const& transitions, + bool assumeOptimalTransitionProbabilities = false); + + /*! + * Computes for each state in the given subsystem an upper bound for the maximal finite expected number of visits. Assumes that there is a strategy that + * yields finite expected number of visits for each state in the subsystem + * @return a vector with the upper bounds. The vector has size subsystem.size() with a -1 in it wherever subsystem is false. + */ + static std::vector computeUpperBoundsOnExpectedVisitingTimes(storm::storage::BitVector const& subsystem, + storm::storage::SparseMatrix const& transitions, + storm::storage::SparseMatrix const& backwardTransitions); +}; +} // namespace storm::modelchecker::multiobjective \ No newline at end of file diff --git a/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp b/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp index a8065df454..e3e71f30f1 100644 --- a/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp +++ b/src/storm/modelchecker/multiobjective/multiObjectiveModelChecking.cpp @@ -2,6 +2,7 @@ #include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.h" +#include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsAchievabilityChecker.h" #include "storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsParetoExplorer.h" #include "storm/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.h" #include "storm/modelchecker/multiobjective/pcaa/SparsePcaaParetoQuery.h" @@ -54,12 +55,18 @@ std::unique_ptr performMultiObjectiveModelChecking(Environment cons switch (method) { case MultiObjectiveMethod::Pcaa: { 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 = DeterministicSchedsParetoExplorer(preprocessorResult); - result = explorer.check(env); - if (env.modelchecker().multi().isExportPlotSet()) { - explorer.exportPlotOfCurrentApproximation(env); + if (preprocessorResult.queryType == preprocessing::SparseMultiObjectivePreprocessorResult::QueryType::Achievability || + preprocessorResult.queryType == preprocessing::SparseMultiObjectivePreprocessorResult::QueryType::Quantitative) { + auto achChecker = DeterministicSchedsAchievabilityChecker(preprocessorResult); + result = achChecker.check(env); + } else { + STORM_LOG_THROW(preprocessorResult.queryType == preprocessing::SparseMultiObjectivePreprocessorResult::QueryType::Pareto, + storm::exceptions::NotImplementedException, "The query type is not implemented with scheduler restrictions."); + auto explorer = DeterministicSchedsParetoExplorer(preprocessorResult); + result = explorer.check(env); + if (env.modelchecker().multi().isExportPlotSet()) { + explorer.exportPlotOfCurrentApproximation(env); + } } } else { std::unique_ptr> query; diff --git a/src/storm/settings/modules/MultiObjectiveSettings.cpp b/src/storm/settings/modules/MultiObjectiveSettings.cpp index fe08ad09af..f2c5cb4ff5 100644 --- a/src/storm/settings/modules/MultiObjectiveSettings.cpp +++ b/src/storm/settings/modules/MultiObjectiveSettings.cpp @@ -76,12 +76,25 @@ MultiObjectiveSettings::MultiObjectiveSettings() : ModuleSettings(moduleName) { storm::settings::OptionBuilder(moduleName, printResultsOptionName, true, "Prints intermediate results of the computation to standard output.") .setIsAdvanced() .build()); - std::vector encodingTypes = {"auto", "classic", "flow"}; - this->addOption(storm::settings::OptionBuilder(moduleName, encodingOptionName, true, "The preferred type of encoding for constraint-based methods.") + this->addOption(storm::settings::OptionBuilder(moduleName, encodingOptionName, true, "The preferred for encoding for constraint-based methods.") .setIsAdvanced() - .addArgument(storm::settings::ArgumentBuilder::createStringArgument("type", "The type.") + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("style", "The main type of encoding.") .setDefaultValueString("auto") - .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(encodingTypes)) + .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator({"auto", "classic", "flow"})) + .build()) + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("constraints", "type of constraints.") + .setDefaultValueString("bigm") + .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator({"bigm", "indicator"})) + .makeOptional() + .build()) + .addArgument(storm::settings::ArgumentBuilder::createStringArgument("bscc", "bscc encoding.") + .setDefaultValueString("flow") + .addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator({"flow", "order"})) + .makeOptional() + .build()) + .addArgument(storm::settings::ArgumentBuilder::createBooleanArgument("redundant", "enable redundant bscc constraints.") + .setDefaultValueBoolean(false) + .makeOptional() .build()) .build()); @@ -171,15 +184,35 @@ bool MultiObjectiveSettings::isPrintResultsSet() const { } bool MultiObjectiveSettings::isClassicEncodingSet() const { - return this->getOption(encodingOptionName).getArgumentByName("type").getValueAsString() == "classic"; + return this->getOption(encodingOptionName).getArgumentByName("style").getValueAsString() == "classic"; } bool MultiObjectiveSettings::isFlowEncodingSet() const { - return this->getOption(encodingOptionName).getArgumentByName("type").getValueAsString() == "flow"; + return this->getOption(encodingOptionName).getArgumentByName("style").getValueAsString() == "flow"; } bool MultiObjectiveSettings::isAutoEncodingSet() const { - return this->getOption(encodingOptionName).getArgumentByName("type").getValueAsString() == "auto"; + return this->getOption(encodingOptionName).getArgumentByName("style").getValueAsString() == "auto"; +} + +bool MultiObjectiveSettings::isBsccDetectionViaFlowConstraintsSet() const { + return this->getOption(encodingOptionName).getArgumentByName("bscc").getValueAsString() == "flow"; +} + +bool MultiObjectiveSettings::isBsccDetectionViaOrderConstraintsSet() const { + return this->getOption(encodingOptionName).getArgumentByName("bscc").getValueAsString() == "order"; +} + +bool MultiObjectiveSettings::isBigMConstraintsSet() const { + return this->getOption(encodingOptionName).getArgumentByName("constraints").getValueAsString() == "bigm"; +} + +bool MultiObjectiveSettings::isIndicatorConstraintsSet() const { + return this->getOption(encodingOptionName).getArgumentByName("constraints").getValueAsString() == "indicator"; +} + +bool MultiObjectiveSettings::isRedundantBsccConstraintsSet() const { + return this->getOption(encodingOptionName).getArgumentByName("redundant").getValueAsBoolean(); } bool MultiObjectiveSettings::isLexicographicModelCheckingSet() const { diff --git a/src/storm/settings/modules/MultiObjectiveSettings.h b/src/storm/settings/modules/MultiObjectiveSettings.h index 4aa753ee30..0f966164bf 100644 --- a/src/storm/settings/modules/MultiObjectiveSettings.h +++ b/src/storm/settings/modules/MultiObjectiveSettings.h @@ -96,11 +96,36 @@ class MultiObjectiveSettings : public ModuleSettings { */ bool isAutoEncodingSet() const; + /*! + * Retrieves whether the encoding for constraint-based methods should use flow constraints for BSCC detection. + */ + bool isBsccDetectionViaFlowConstraintsSet() const; + + /*! + * Retrieves whether the encoding for constraint-based methods should use order constraints for BSCC detection. + */ + bool isBsccDetectionViaOrderConstraintsSet() const; + + /*! + * Retrieves whether the encoding for constraint-based methods should use BigM constraints + */ + bool isBigMConstraintsSet() const; + + /*! + * Retrieves whether the encoding for constraint-based methods should use indicator constraints + */ + bool isIndicatorConstraintsSet() const; + /*! * Retrieves whether lexicographic model checking has been set */ bool isLexicographicModelCheckingSet() const; + /*! + * Retrieves whether redundant BSCC constraints are to be added + */ + bool isRedundantBsccConstraintsSet() const; + /*! * Checks whether the settings are consistent. If they are inconsistent, an exception is thrown. * diff --git a/src/test/storm/modelchecker/multiobjective/MultiObjectiveSchedRestModelCheckerTest.cpp b/src/test/storm/modelchecker/multiobjective/MultiObjectiveSchedRestModelCheckerTest.cpp index e6a323ca99..755a60f403 100644 --- a/src/test/storm/modelchecker/multiobjective/MultiObjectiveSchedRestModelCheckerTest.cpp +++ b/src/test/storm/modelchecker/multiobjective/MultiObjectiveSchedRestModelCheckerTest.cpp @@ -1,76 +1,181 @@ #include "storm-config.h" #include "test/storm_gtest.h" -#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" -#include "storm/modelchecker/multiobjective/multiObjectiveModelChecking.h" - #include "storm-parsers/api/storm-parsers.h" #include "storm/api/storm.h" #include "storm/environment/Environment.h" +#include "storm/environment/modelchecker/MultiObjectiveModelCheckerEnvironment.h" #include "storm/exceptions/InvalidOperationException.h" +#include "storm/modelchecker/multiobjective/multiObjectiveModelChecking.h" #include "storm/modelchecker/results/ExplicitParetoCurveCheckResult.h" -#include "storm/models/sparse/MarkovAutomaton.h" +#include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/models/sparse/Mdp.h" -#include "storm/settings/SettingsManager.h" -#include "storm/settings/modules/GeneralSettings.h" #include "storm/storage/SchedulerClass.h" -#include "storm/storage/jani/Property.h" #ifdef STORM_HAVE_Z3_OPTIMIZE namespace { -storm::Environment getPositionalDeterministicEnvironment() { - storm::Environment env; - env.modelchecker().multi().setSchedulerRestriction(storm::storage::SchedulerClass().setPositional().setIsDeterministic()); - return env; -} +class FlowBigMEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); + env.modelchecker().multi().setUseIndicatorConstraints(false); + env.modelchecker().multi().setUseBsccOrderEncoding(false); // not relevant for Flow + env.modelchecker().multi().setUseRedundantBsccConstraints(false); // not relevant for Flow + return env; + } +}; + +class FlowIndicatorEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); + env.modelchecker().multi().setUseIndicatorConstraints(true); + env.modelchecker().multi().setUseBsccOrderEncoding(false); // not relevant for Flow + env.modelchecker().multi().setUseRedundantBsccConstraints(false); // not relevant for Flow + return env; + } +}; + +class BigMEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(false); + env.modelchecker().multi().setUseBsccOrderEncoding(false); + env.modelchecker().multi().setUseRedundantBsccConstraints(false); + return env; + } +}; + +class IndicatorEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(true); + env.modelchecker().multi().setUseBsccOrderEncoding(false); + env.modelchecker().multi().setUseRedundantBsccConstraints(false); + return env; + } +}; + +class BigMOrderEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(false); + env.modelchecker().multi().setUseBsccOrderEncoding(true); + env.modelchecker().multi().setUseRedundantBsccConstraints(false); + return env; + } +}; + +class IndicatorOrderEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(true); + env.modelchecker().multi().setUseBsccOrderEncoding(true); + env.modelchecker().multi().setUseRedundantBsccConstraints(false); + return env; + } +}; + +class RedundantEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(false); + env.modelchecker().multi().setUseBsccOrderEncoding(false); + env.modelchecker().multi().setUseRedundantBsccConstraints(true); + return env; + } +}; + +class RedundantOrderEnvironment { + public: + static storm::Environment getEnv() { + storm::Environment env; + env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); + env.modelchecker().multi().setUseIndicatorConstraints(false); + env.modelchecker().multi().setUseBsccOrderEncoding(true); + env.modelchecker().multi().setUseRedundantBsccConstraints(true); + return env; + } +}; -storm::Environment getGoalDeterministicEnvironment() { - storm::Environment env; - env.modelchecker().multi().setSchedulerRestriction( - storm::storage::SchedulerClass().setMemoryPattern(storm::storage::SchedulerClass::MemoryPattern::GoalMemory).setIsDeterministic()); - return env; -} +template +class MultiObjectiveSchedRestModelCheckerTest : public ::testing::Test { + public: + typedef storm::RationalNumber ValueType; + + bool isFlowEncoding() const { + return TestType::getEnv().modelchecker().multi().getEncodingType() == storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow; + } + + storm::Environment getPositionalDeterministicEnvironment() { + auto env = TestType::getEnv(); + env.modelchecker().multi().setSchedulerRestriction(storm::storage::SchedulerClass().setPositional().setIsDeterministic()); + return env; + } + + storm::Environment getGoalDeterministicEnvironment() { + auto env = TestType::getEnv(); + env.modelchecker().multi().setSchedulerRestriction( + storm::storage::SchedulerClass().setMemoryPattern(storm::storage::SchedulerClass::MemoryPattern::GoalMemory).setIsDeterministic()); + return env; + } + + typedef std::vector Point; -typedef std::vector Point; - -std::vector parsePoints(std::vector const& input) { - std::vector result; - for (auto const& i : input) { - Point currPoint; - std::size_t pos1 = 0; - std::size_t pos2 = i.find(","); - while (pos2 != std::string::npos) { - currPoint.push_back(storm::utility::convertNumber(i.substr(pos1, pos2 - pos1))); - pos1 = pos2 + 1; - pos2 = i.find(",", pos1); + ValueType parseNumber(std::string const& input) const { + return storm::utility::convertNumber(input); + } + + std::vector parsePoints(std::vector const& input) { + std::vector result; + for (auto const& i : input) { + Point currPoint; + std::size_t pos1 = 0; + std::size_t pos2 = i.find(","); + while (pos2 != std::string::npos) { + currPoint.push_back(parseNumber(i.substr(pos1, pos2 - pos1))); + pos1 = pos2 + 1; + pos2 = i.find(",", pos1); + } + currPoint.push_back(parseNumber(i.substr(pos1))); + result.push_back(currPoint); } - currPoint.push_back(storm::utility::convertNumber(i.substr(pos1))); - result.push_back(currPoint); + return result; } - return result; -} -std::set setMinus(std::vector const& lhs, std::vector const& rhs) { - std::set result(lhs.begin(), lhs.end()); - for (auto const& r : rhs) { - for (auto lIt = result.begin(); lIt != result.end(); ++lIt) { - if (*lIt == r) { - result.erase(lIt); - break; + std::set setMinus(std::vector const& lhs, std::vector const& rhs) { + std::set result(lhs.begin(), lhs.end()); + for (auto const& r : rhs) { + for (auto lIt = result.begin(); lIt != result.end(); ++lIt) { + if (*lIt == r) { + result.erase(lIt); + break; + } } } + return result; } - return result; -} -std::string toString(std::vector pointset, bool asDouble) { - std::stringstream s; - for (auto const& p : pointset) { + std::string toString(Point point, bool asDouble) { + std::stringstream s; s << "["; bool first = true; - for (auto const& pi : p) { + for (auto const& pi : point) { if (first) { first = false; } else { @@ -82,187 +187,198 @@ std::string toString(std::vector pointset, bool asDouble) { s << pi; } } - s << "]\n"; + s << "]"; + return s.str(); } - return s.str(); -} -TEST(MultiObjectiveSchedRestModelCheckerTest, steps) { + std::string getDiffString(std::vector const& expected, std::vector const& actual) { + std::stringstream stream; + stream << "Unexpected set of Points:\n"; + stream << " Expected | Actual | Point\n"; + for (auto const& p : expected) { + if (std::find(actual.begin(), actual.end(), p) != actual.end()) { + stream << " yes | yes | " << toString(p, true) << "\n"; + } else { + stream << " --> yes | no | " << toString(p, true) << "\n"; + } + } + for (auto const& p : actual) { + if (std::find(expected.begin(), expected.end(), p) == expected.end()) { + stream << " --> no | yes | " << toString(p, true) << "\n"; + } + } + return stream.str(); + } + + bool isSame(std::vector const& expected, std::vector const& actual, std::string& diffString) { + if (expected.size() != actual.size()) { + diffString = getDiffString(expected, actual); + return false; + } + for (auto const& p : expected) { + if (std::find(actual.begin(), actual.end(), p) == actual.end()) { + diffString = getDiffString(expected, actual); + return false; + } + } + diffString = ""; + return true; + } + + template + bool testParetoFormula(storm::Environment const& env, SparseModelType const& model, storm::logic::Formula const& formula, + std::vector const& expected, std::string& errorString) { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *model, formula.asMultiObjectiveFormula()); + if (!result->isParetoCurveCheckResult()) { + errorString = "Not a ParetoCurveCheckResult"; + return false; + } + return isSame(expected, result->template asExplicitParetoCurveCheckResult().getPoints(), errorString); + } +}; + +typedef ::testing::Types + TestingTypes; + +TYPED_TEST_SUITE(MultiObjectiveSchedRestModelCheckerTest, TestingTypes, ); + +TYPED_TEST(MultiObjectiveSchedRestModelCheckerTest, steps) { + typedef typename TestFixture::ValueType ValueType; + typedef typename TestFixture::Point Point; + std::string programFile = STORM_TEST_RESOURCES_DIR "/mdp/multiobj_stairs.nm"; std::string constantsString = "N=3"; - std::string formulasAsString = "multi(Pmax=? [ F y=1], Pmax=? [ F y=2 ])"; + std::string formulasAsString = "multi(Pmax=? [ F y=1], Pmax=? [ F y=2 ]);"; + formulasAsString += "multi(P>=0.375 [ F y=1], Pmax>=0.5 [ F y=2 ]);"; + formulasAsString += "multi(P>=0.4 [ F y=1], Pmax>=0.4 [ F y=2 ]);"; + formulasAsString += "multi(Pmax=? [ F y=1], Pmax>=0.4 [ F y=2 ]);"; + formulasAsString += "multi(Pmax=? [ F y=1], Pmax>=0.9 [ F y=2 ]);"; // programm, model, formula storm::prism::Program program = storm::api::parseProgram(programFile); program = storm::utility::prism::preprocess(program, constantsString); std::vector> formulas = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program)); - std::shared_ptr> mdp = - storm::api::buildSparseModel(program, formulas)->as>(); - std::vector expected, actual; - std::set incorrectPoints, missingPoints; - std::unique_ptr result; - - storm::Environment env = getPositionalDeterministicEnvironment(); - - expected = parsePoints({"0.875,0", "0,0.875", "0.125,0.75", "0.25,0.625", "0.375,0.5", "0.5,0.375", "0.625,0.25", "0.75,0.125"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[0]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[0]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + std::shared_ptr> mdp = + storm::api::buildSparseModel(program, formulas)->template as>(); + std::string errorString; // storage for error reporting + + storm::Environment env = this->getPositionalDeterministicEnvironment(); + uint64_t formulaIndex = 0; + { + auto expected = this->parsePoints({"0.875,0", "0,0.875", "0.125,0.75", "0.25,0.625", "0.375,0.5", "0.5,0.375", "0.625,0.25", "0.75,0.125"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } + ++formulaIndex; + { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQualitativeCheckResult()); + EXPECT_TRUE(result->asExplicitQualitativeCheckResult()[*mdp->getInitialStates().begin()]); + } + ++formulaIndex; + { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQualitativeCheckResult()); + EXPECT_FALSE(result->asExplicitQualitativeCheckResult()[*mdp->getInitialStates().begin()]); + } + ++formulaIndex; + { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQuantitativeCheckResult()); + auto expected = storm::utility::convertNumber("0.375"); + EXPECT_EQ(result->template asExplicitQuantitativeCheckResult()[*mdp->getInitialStates().begin()], expected); + } + ++formulaIndex; + { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQualitativeCheckResult()); + EXPECT_FALSE(result->asExplicitQualitativeCheckResult()[*mdp->getInitialStates().begin()]); + } } -TEST(MultiObjectiveSchedRestModelCheckerTest, mecs) { +TYPED_TEST(MultiObjectiveSchedRestModelCheckerTest, mecs) { + typedef typename TestFixture::ValueType ValueType; + typedef typename TestFixture::Point Point; std::string programFile = STORM_TEST_RESOURCES_DIR "/mdp/multiobj_mecs.nm"; std::string constantsString = ""; std::string formulasAsString = "multi(Pmin=? [ F x=3], Pmax=? [ F x=4 ]);"; formulasAsString += "\nmulti(R{\"test\"}min=? [C], Pmin=? [F \"t1\"]);"; formulasAsString += "\nmulti(R{\"test\"}min=? [C], Pmin=? [F \"t2\"]);"; formulasAsString += "\nmulti(R{\"test\"}min=? [C], Pmin=? [F \"t2\"], Pmax=? [F \"t1\"]);"; + formulasAsString += "\nmulti(R{\"test\"}<=0 [C], P<=1 [F \"t2\"], P>=0 [F \"t1\"]);"; + formulasAsString += "\nmulti(R{\"test\"}<=1 [C], P<=1 [F \"t2\"], P>=0.1 [F \"t1\"]);"; // programm, model, formula storm::prism::Program program = storm::api::parseProgram(programFile); program = storm::utility::prism::preprocess(program, constantsString); std::vector> formulas = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program)); - std::shared_ptr> mdp = - storm::api::buildSparseModel(program, formulas)->as>(); - std::vector expected, actual; - std::set incorrectPoints, missingPoints; - std::unique_ptr result; + std::shared_ptr> mdp = + storm::api::buildSparseModel(program, formulas)->template as>(); + std::string errorString; // storage for error reporting + std::vector expected; - storm::Environment env = getPositionalDeterministicEnvironment(); + storm::Environment env = this->getPositionalDeterministicEnvironment(); uint64_t formulaIndex = 0; - expected = parsePoints({"0.5,0.5", "1,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - STORM_SILENT_EXPECT_THROW( - storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), - storm::exceptions::InvalidOperationException); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0.5,0.5", "1,1"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } ++formulaIndex; - expected = parsePoints({"0,0"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,0"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; ++formulaIndex; - expected = parsePoints({"0,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - STORM_SILENT_EXPECT_THROW( - storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), - storm::exceptions::InvalidOperationException); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,1"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } ++formulaIndex; - expected = parsePoints({"0,1,0", "2,1,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - STORM_SILENT_EXPECT_THROW( - storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), - storm::exceptions::InvalidOperationException); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,1,0", "2,1,1"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } + ++formulaIndex; + { + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW( + storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), + storm::exceptions::InvalidOperationException); + } else { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQualitativeCheckResult()); + EXPECT_TRUE(result->asExplicitQualitativeCheckResult()[*mdp->getInitialStates().begin()]); + } + } + ++formulaIndex; + { + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW( + storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), + storm::exceptions::InvalidOperationException); + } else { + auto result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); + ASSERT_TRUE(result->isExplicitQualitativeCheckResult()); + EXPECT_FALSE(result->asExplicitQualitativeCheckResult()[*mdp->getInitialStates().begin()]); + } + } } -TEST(MultiObjectiveSchedRestModelCheckerTest, compromise) { +TYPED_TEST(MultiObjectiveSchedRestModelCheckerTest, compromise) { + typedef typename TestFixture::ValueType ValueType; + typedef typename TestFixture::Point Point; std::string programFile = STORM_TEST_RESOURCES_DIR "/mdp/multiobj_compromise.nm"; std::string constantsString = ""; std::string formulasAsString = "multi(Pmax=? [ F x=1], Pmax=? [ F x=2 ]);"; @@ -274,260 +390,92 @@ TEST(MultiObjectiveSchedRestModelCheckerTest, compromise) { program = storm::utility::prism::preprocess(program, constantsString); std::vector> formulas = storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program)); - std::shared_ptr> mdp = - storm::api::buildSparseModel(program, formulas)->as>(); - std::vector expected, actual; - std::set incorrectPoints, missingPoints; - std::unique_ptr result; + std::shared_ptr> mdp = + storm::api::buildSparseModel(program, formulas)->template as>(); + std::string errorString; // storage for error reporting + std::vector expected; - storm::Environment env = getPositionalDeterministicEnvironment(); + storm::Environment env = this->getPositionalDeterministicEnvironment(); uint64_t formulaIndex = 0; - expected = parsePoints({"1,0", "0,1", "0.3,3/7"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - STORM_SILENT_EXPECT_THROW( - storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), - storm::exceptions::InvalidOperationException); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"1,0", "0,1", "0.3,3/7"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } ++formulaIndex; - expected = parsePoints({"0,0.3", "2,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,0.3", "2,1"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; ++formulaIndex; - expected = parsePoints({"1,0,0", "0,1,0", "0.3,3/7,4/7"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - STORM_SILENT_EXPECT_THROW( - storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()), - storm::exceptions::InvalidOperationException); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"1,0,0", "0,1,0", "0.3,3/7,4/7"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } ++formulaIndex; - expected = parsePoints({"0,1,0", "0,3/7,4/7"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env = getGoalDeterministicEnvironment(); + expected = this->parsePoints({"0,1,0", "0,3/7,4/7"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + + env = this->getGoalDeterministicEnvironment(); formulaIndex = 0; - expected = parsePoints({"1,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"1,1"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; ++formulaIndex; - expected = parsePoints({"0,0.3", "2,1"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,0.3", "2,1"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; ++formulaIndex; - expected = parsePoints({"1,1,0", "1,3/7,4/7"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"1,1,0", "1,3/7,4/7"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; ++formulaIndex; - expected = parsePoints({"0,1,0", "0,3/7,4/7"}); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Flow); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - - env.modelchecker().multi().setEncodingType(storm::MultiObjectiveModelCheckerEnvironment::EncodingType::Classic); - result = storm::modelchecker::multiobjective::performMultiObjectiveModelChecking(env, *mdp, formulas[formulaIndex]->asMultiObjectiveFormula()); - ASSERT_TRUE(result->isParetoCurveCheckResult()); - actual = result->asExplicitParetoCurveCheckResult().getPoints(); - missingPoints = setMinus(expected, actual); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the expected solution are missing:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); - incorrectPoints = setMinus(actual, expected); - EXPECT_TRUE(incorrectPoints.empty()) << "Some points of the returned solution are not expected:\n" - << "Expected:\n" - << toString(expected, true) << "Actual:\n" - << toString(actual, true); + expected = this->parsePoints({"0,1,0", "0,3/7,4/7"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; } + +TYPED_TEST(MultiObjectiveSchedRestModelCheckerTest, infrew) { + typedef typename TestFixture::ValueType ValueType; + typedef typename TestFixture::Point Point; + std::string programFile = STORM_TEST_RESOURCES_DIR "/mdp/multiobj_infrew.nm"; + std::string constantsString = ""; + std::string formulasAsString = "multi(R{\"one\"}min=? [ F x=2], R{\"two\"}min=? [ C ]);"; + // programm, model, formula + storm::prism::Program program = storm::api::parseProgram(programFile); + program = storm::utility::prism::preprocess(program, constantsString); + std::vector> formulas = + storm::api::extractFormulasFromProperties(storm::api::parsePropertiesForPrismProgram(formulasAsString, program)); + std::shared_ptr> mdp = + storm::api::buildSparseModel(program, formulas)->template as>(); + std::string errorString; // storage for error reporting + std::vector expected; + + storm::Environment env = this->getPositionalDeterministicEnvironment(); + + uint64_t formulaIndex = 0; + expected = this->parsePoints({"10,1"}); + if (this->isFlowEncoding()) { + STORM_SILENT_EXPECT_THROW(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString), + storm::exceptions::InvalidOperationException); + } else { + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; + } + + env = this->getGoalDeterministicEnvironment(); + + formulaIndex = 0; + expected = this->parsePoints({"0,1"}); + EXPECT_TRUE(this->testParetoFormula(env, mdp, *formulas[formulaIndex], expected, errorString)) << errorString; +} + } // namespace -#endif /* STORM_HAVE_Z3_OPTIMIZE */ +#endif /* defined STORM_HAVE_Z3_OPTIMIZE */