diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index 71d9f16b4..80e99295b 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -167,7 +167,6 @@ namespace storm { template std::map getSubEndComponents(storm::storage::SparseMatrix const& mecTransitions) { auto backwardTransitions = mecTransitions.transpose(true); - auto const& groups = mecTransitions.getRowGroupIndices(); std::map unprocessed, processed; storm::storage::BitVector allStates(mecTransitions.getRowGroupCount()); storm::storage::BitVector allChoices(mecTransitions.getRowCount()); @@ -237,7 +236,7 @@ namespace storm { uint64_t subEcChoice = toGlobalChoiceIndexMapping[localSubEcChoiceIndex]; for (auto const& objHelper : objectiveHelper) { if (objHelper.getChoiceValueOffsets().count(subEcChoice) > 0) { - STORM_LOG_ASSERT(!storm::utility::isZero(objHelper.getChoiceValueOffsets().at(subEcChoice)), "Expectet non-zero choice-value offset."); + STORM_LOG_ASSERT(!storm::utility::isZero(objHelper.getChoiceValueOffsets().at(subEcChoice)), "Expected non-zero choice-value offset."); subEcChoicesWithValueZero.set(localSubEcChoiceIndex, false); break; } @@ -275,12 +274,16 @@ namespace storm { } } // Assert that the ecVar is one iff the sum over the zero-value-choice variables equals the number of states in this ec - lpModel->addConstraint("", ecVar >= one + storm::expressions::sum(ecChoiceVars) - lpModel->getConstant(storm::utility::convertNumber(numSubEcStatesWithMultipleChoices))); + storm::expressions::Expression ecVarLowerBound = one - lpModel->getConstant(storm::utility::convertNumber(numSubEcStatesWithMultipleChoices)).simplify(); + if (!ecChoiceVars.empty()) { + ecVarLowerBound = ecVarLowerBound + storm::expressions::sum(ecChoiceVars); + } + lpModel->addConstraint("", ecVar >= ecVarLowerBound); } } } - STORM_LOG_TRACE("Found " << ecCounter << " End components."); + STORM_LOG_TRACE("Found " << ecCounter << " end components."); return result; } @@ -353,24 +356,8 @@ namespace storm { STORM_LOG_ERROR_COND(storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), nonBottomStates, bottomStates).full(), "End components not yet treated correctly."); // Compute upper bounds for each state - std::vector visitingTimesUpperBounds; - { - // TODO: this doesn't work if there are end components. - storm::storage::SparseMatrix transitions = model.getTransitionMatrix().getSubmatrix(true, nonBottomStates, nonBottomStates); - std::vector probabilitiesToBottomStates = model.getTransitionMatrix().getConstrainedRowGroupSumVector(nonBottomStates, bottomStates); - auto subsystemBounds = storm::modelchecker::helper::BaierUpperRewardBoundsComputer::computeUpperBoundOnExpectedVisitingTimes(transitions, probabilitiesToBottomStates); - uint64_t subsystemState = 0; - visitingTimesUpperBounds.reserve(bottomStates.size()); - for (uint64_t state = 0; state < bottomStates.size(); ++state) { - if (bottomStates.get(state)) { - visitingTimesUpperBounds.push_back(storm::utility::zero()); - } else { - visitingTimesUpperBounds.push_back(subsystemBounds[subsystemState]); - ++subsystemState; - } - } - assert(subsystemState == subsystemBounds.size()); - } + std::vector visitingTimesUpperBounds = DeterministicSchedsObjectiveHelper::computeUpperBoundOnExpectedVisitingTimes(model.getTransitionMatrix(), bottomStates, nonBottomStates, hasEndComponents); + // create choiceValue variables and assert deterministic ones. std::vector choiceValVars(model.getNumberOfChoices()); for (auto const& state : nonBottomStates) { @@ -561,12 +548,12 @@ namespace storm { if (objectiveHelper[objIndex].minimizing()) { // TODO: these are optional if (isMinNegativeEncoding()) { - lpModel->addConstraint("", stateVars[state] >= (ecVar.getExpression() - one) * objectiveHelper[objIndex].getUpperValueBoundAtState(env, state)); + lpModel->addConstraint("", stateVars[state] >= (ecVar.getExpression() - one) * lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state))); } else { - lpModel->addConstraint("", stateVars[state] <= (one - ecVar.getExpression()) * objectiveHelper[objIndex].getUpperValueBoundAtState(env, state)); + lpModel->addConstraint("", stateVars[state] <= (one - ecVar.getExpression()) * lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state))); } } else { - lpModel->addConstraint("", stateVars[state] <= (one - ecVar.getExpression()) * objectiveHelper[objIndex].getUpperValueBoundAtState(env, state)); + lpModel->addConstraint("", stateVars[state] <= (one - ecVar.getExpression()) * lpModel->getConstant(objectiveHelper[objIndex].getUpperValueBoundAtState(env, state))); } } } diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp index b1f9ece62..756963a07 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp @@ -7,14 +7,17 @@ #include "storm/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h" #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/storage/BitVector.h" +#include "storm/storage/MaximalEndComponentDecomposition.h" #include "storm/utility/graph.h" #include "storm/utility/FilteredRewardModel.h" #include "storm/utility/vector.h" #include "storm/logic/Formulas.h" #include "storm/logic/CloneVisitor.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/transformer/EndComponentEliminator.h" #include "storm/exceptions/UnexpectedException.h" @@ -192,21 +195,53 @@ namespace storm { } } + 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); + } + 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()); + } + template typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getUpperValueBoundAtState(Environment const& env, uint64_t state) const{ - //return objective.upperResultBound.get(); if (!upperResultBounds) { upperResultBounds = computeValueBounds(env, false, model, *objective.formula); - storm::utility::vector::clip(upperResultBounds.get(), objective.lowerResultBound, objective.upperResultBound); - STORM_LOG_THROW(!storm::utility::vector::hasInfinityEntry(upperResultBounds.get()), storm::exceptions::NotSupportedException, "The upper bound for objective " << *objective.originalFormula << " is infinity at some state. This is not supported."); + auto upperResultBound = objective.upperResultBound; + if (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; + } + } + storm::utility::vector::clip(upperResultBounds.get(), objective.lowerResultBound, upperResultBound); } return upperResultBounds.get()[state]; - } template typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper::getLowerValueBoundAtState(Environment const& env, uint64_t state) const{ - // return objective.lowerResultBound.get(); if (!lowerResultBounds) { lowerResultBounds = computeValueBounds(env, true, model, *objective.formula); storm::utility::vector::clip(lowerResultBounds.get(), objective.lowerResultBound, objective.upperResultBound); @@ -224,6 +259,79 @@ namespace storm { bool DeterministicSchedsObjectiveHelper::isTotalRewardObjective() const { return objective.formula->isRewardOperatorFormula() && objective.formula->getSubformula().isTotalRewardFormula(); } + + 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, true); + transitions = std::move(mecElimination.matrix); + modelToSubsystemStateMapping = std::move(mecElimination.oldToNewStateMapping); + probabilitiesToBottomStates.reserve(transitions.getRowCount()); + for (uint64_t row = 0; row < transitions.getRowCount(); ++row) { + probabilitiesToBottomStates.push_back(modelTransitions.getConstrainedRowSum(mecElimination.newToOldRowMapping[row], bottomStates)); + } + // replace 'selfloop probability' for mec states by 1-lpath + for (auto const& mec : mecs) { + ValueType lpath = storm::utility::one(); + for (auto const& stateChoices : mec) { + ValueType minProb = storm::utility::one(); + for (auto const& choice : stateChoices.second) { + for (auto const& transition : modelTransitions.getRow(choice)) { + if (!storm::utility::isZero(transition.getValue())) { + minProb = std::min(minProb, transition.getValue()); + } + } + } + lpath *= minProb; + } + STORM_LOG_ASSERT(!storm::utility::isZero(lpath), "unexpected value of lpath"); + uint64_t mecState = mecElimination.oldToNewStateMapping[mec.begin()->first]; + bool foundEntry = false; + for (uint64_t mecChoice = transitions.getRowGroupIndices()[mecState]; mecChoice < transitions.getRowGroupIndices()[mecState + 1]; ++mecChoice) { + if (transitions.getRow(mecChoice).getNumberOfEntries() == 1) { + auto& entry = *transitions.getRow(mecChoice).begin(); + if (entry.getColumn() == mecState && storm::utility::isOne(entry.getValue())) { + entry.setValue(storm::utility::one() - lpath); + foundEntry = true; + probabilitiesToBottomStates[mecChoice] = lpath; + break; + } + } + } + STORM_LOG_THROW(foundEntry, storm::exceptions::UnexpectedException, "Unable to find self loop entry at mec state."); + } + } else { + transitions = modelTransitions.getSubmatrix(true, nonBottomStates, nonBottomStates); + probabilitiesToBottomStates = modelTransitions.getConstrainedRowGroupSumVector(nonBottomStates, bottomStates); + } + + auto subsystemBounds = storm::modelchecker::helper::BaierUpperRewardBoundsComputer::computeUpperBoundOnExpectedVisitingTimes(transitions, probabilitiesToBottomStates); + uint64_t subsystemState = 0; + + 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()); + } else { + if (modelToSubsystemStateMapping) { + visitingTimesUpperBounds.push_back(subsystemBounds[modelToSubsystemStateMapping.get()[state]]); + } else { + visitingTimesUpperBounds.push_back(subsystemBounds[subsystemState]); + } + ++subsystemState; + } + } + assert(subsystemState == subsystemBounds.size()); + } + template class DeterministicSchedsObjectiveHelper>; template class DeterministicSchedsObjectiveHelper>; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h index c93ef6086..4250c4ebe 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.h @@ -4,6 +4,9 @@ #include #include "storm/modelchecker/multiobjective/Objective.h" +#include "storm/storage/SparseMatrix.h" +#include "storm/storage/BitVector.h" + namespace storm { @@ -41,6 +44,8 @@ namespace storm { */ bool isTotalRewardObjective() const; + static std::vector computeUpperBoundOnExpectedVisitingTimes(storm::storage::SparseMatrix const& modelTransitions, storm::storage::BitVector const& bottomStates, storm::storage::BitVector const& nonBottomStates, bool hasEndComponents); + private: mutable boost::optional> schedulerIndependentStateValues;