diff --git a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp index 9656f28f6..5519e5b55 100644 --- a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp +++ b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp @@ -530,46 +530,6 @@ namespace storm { auto const& rewModel = result.preprocessedModel->getRewardModel(result.objectives[objIndex].formula->asRewardOperatorFormula().getRewardModelName()); auto actionRewards = rewModel.getTotalRewardVector(transitions); - // TODO: Consider cumulative reward formulas less naively - if (result.objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) { - auto const& cumulativeRewardFormula = result.objectives[objIndex].formula->getSubformula().asCumulativeRewardFormula(); - ValueType rewardBound = cumulativeRewardFormula.template getBound(); - if (cumulativeRewardFormula.getTimeBoundReference().isRewardBound()) { - auto const& costModel = result.preprocessedModel->getRewardModel(cumulativeRewardFormula.getTimeBoundReference().getRewardName()); - if (!costModel.hasTransitionRewards()) { - auto actionCosts = costModel.getTotalRewardVector(transitions); - typename SparseModelType::ValueType largestRewardPerCost = storm::utility::zero(); - bool isFinite = true; - for (auto rewIt = actionRewards.begin(), costIt = actionCosts.begin(); rewIt != actionRewards.end(); ++rewIt, ++costIt) { - if (!storm::utility::isZero(*rewIt)) { - if (storm::utility::isZero(*costIt)) { - isFinite = false; - break; - } - ValueType rewardPerCost = *rewIt / *costIt; - largestRewardPerCost = std::max(largestRewardPerCost, rewardPerCost); - } - } - - if (isFinite) { - ValueType newResultBound = largestRewardPerCost * rewardBound; - if (upperBound) { - upperBound = std::min(upperBound.get(), newResultBound); - } else { - upperBound = newResultBound; - } - } - } - } else { - ValueType newResultBound = (*std::max_element(actionRewards.begin(), actionRewards.end())) * rewardBound; - if (upperBound) { - upperBound = std::min(upperBound.get(), newResultBound); - } else { - upperBound = newResultBound; - } - } - } - if (result.objectives[objIndex].formula->getSubformula().isTotalRewardFormula() || result.objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) { // We have to eliminate ECs here to treat zero-reward ECs @@ -590,11 +550,21 @@ namespace storm { rew0StateProbs.reserve(ecElimRes.matrix.getRowCount()); for (uint64_t state = 0; state < allStates.size(); ++ state) { for (uint64_t choice = ecElimRes.matrix.getRowGroupIndices()[state]; choice < ecElimRes.matrix.getRowGroupIndices()[state + 1]; ++choice) { - ValueType outProb = storm::utility::one() - ecElimRes.matrix.getRowSum(choice); - if (!storm::utility::isZero(outProb)) { - outStates.set(state, true); + // Check whether the choice lead to a state with expRew 0 in the original model + bool isOutChoice = false; + uint64_t originalModelChoice = ecElimRes.newToOldRowMapping[choice]; + for (auto const& entry : transitions.getRow(originalModelChoice)) { + if (!expRewGreater0EStates.get(entry.getColumn())) { + isOutChoice = true; + outStates.set(state, true); + rew0StateProbs.push_back(storm::utility::one() - ecElimRes.matrix.getRowSum(choice)); + assert (!storm::utility::isZero(rew0StateProbs.back())); + break; + } + } + if (!isOutChoice) { + rew0StateProbs.push_back(storm::utility::zero()); } - rew0StateProbs.push_back(outProb); } } diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp index 2f76e3dae..1887e86ce 100644 --- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp +++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpRewardBoundedPcaaWeightVectorChecker.cpp @@ -34,6 +34,12 @@ namespace storm { STORM_LOG_THROW(preprocessorResult.rewardFinitenessType == SparseMultiObjectivePreprocessorResult::RewardFinitenessType::AllFinite, storm::exceptions::NotSupportedException, "There is a scheduler that yields infinite reward for one objective. This is not supported."); STORM_LOG_THROW(preprocessorResult.preprocessedModel->getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::NotSupportedException, "The model has multiple initial states."); + // Update the objective bounds with what the reward unfolding can compute + for (uint64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) { + this->objectives[objIndex].lowerResultBound = rewardUnfolding.getLowerObjectiveBound(objIndex); + this->objectives[objIndex].upperResultBound = rewardUnfolding.getUpperObjectiveBound(objIndex); + } + numCheckedEpochs = 0; numChecks = 0; } diff --git a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 8712debb7..bd1ef4f4b 100644 --- a/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -67,7 +67,7 @@ namespace storm { } auto formula = std::make_shared(checkTask.getFormula().asSharedPointer(), opInfo); helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); - auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(helper::SolutionType::UntilProbabilities, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); + auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } else { STORM_LOG_THROW(!pathFormula.hasLowerBound() && pathFormula.hasUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have single upper time bound."); @@ -146,7 +146,7 @@ namespace storm { } auto formula = std::make_shared(checkTask.getFormula().asSharedPointer(), checkTask.getRewardModel(), opInfo); helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding(this->getModel(), formula); - auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(helper::SolutionType::ExpectedRewards, checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); + auto numericResult = storm::modelchecker::helper::SparseMdpPrctlHelper::computeRewardBoundedValues(checkTask.getOptimizationDirection(), rewardUnfolding, this->getModel().getInitialStates(), *minMaxLinearEquationSolverFactory); return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); } else { diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index e63fb7374..204590e6d 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -134,7 +134,7 @@ namespace storm { } template - std::vector analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ValueType const& precision, SolutionType const& type) { + std::vector analyzeNonTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ValueType const& precision, boost::optional const& lowerBound, boost::optional const& upperBound) { // Update some data for the case that the Matrix has changed if (epochModel.epochMatrixChanged) { @@ -146,14 +146,13 @@ namespace storm { minMaxSolver->setCachingEnabled(true); minMaxSolver->setTrackScheduler(true); auto req = minMaxSolver->getRequirements(dir); - minMaxSolver->setLowerBound(storm::utility::zero()); - req.clearLowerBounds(); - if (type == SolutionType::UntilProbabilities) { - minMaxSolver->setUpperBound(storm::utility::one()); + if (lowerBound) { + minMaxSolver->setLowerBound(lowerBound.get()); + req.clearLowerBounds(); + } + if (upperBound) { + minMaxSolver->setUpperBound(upperBound.get()); req.clearUpperBounds(); - } else if (type == SolutionType::ExpectedRewards) { - // TODO - STORM_LOG_WARN_COND(!req.requiresUpperBounds(), "Upper bounds for expected reward are not specified."); } STORM_LOG_THROW(req.empty(), storm::exceptions::UncheckedRequirementException, "At least one requirement was not checked."); minMaxSolver->setRequirementsChecked(); @@ -182,8 +181,14 @@ namespace storm { } template - std::map SparseMdpPrctlHelper::computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::map SparseMdpPrctlHelper::computeRewardBoundedValues(OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { storm::utility::Stopwatch swAll(true), swBuild, swCheck; + + // Get lower and upper bounds for the solution. + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + + // Initialize epoch models auto initEpoch = rewardUnfolding.getStartEpoch(); auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch); @@ -205,7 +210,7 @@ namespace storm { if (epochModel.epochMatrix.getEntryCount() == 0) { rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialEpochModel(dir, epochModel)); } else { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialEpochModel(dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, precision, type)); + rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialEpochModel(dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, precision, lowerBound, upperBound)); } swCheck.stop(); ++numCheckedEpochs; diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h index 1990a4873..96243c555 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.h @@ -38,7 +38,7 @@ namespace storm { static std::vector computeStepBoundedUntilProbabilities(storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint()); - static std::map computeRewardBoundedValues(SolutionType const& type, OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::map computeRewardBoundedValues(OptimizationDirection dir, rewardbounded::MultiDimensionalRewardUnfolding& rewardUnfolding, storm::storage::BitVector const& initialStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); static std::vector computeNextProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index ee09ba970..93a9afd02 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -9,6 +9,8 @@ #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" + #include "storm/transformer/EndComponentEliminator.h" #include "storm/exceptions/UnexpectedException.h" @@ -548,6 +550,103 @@ namespace storm { return precision / storm::utility::convertNumber(sumOfDimensions); } + template + boost::optional MultiDimensionalRewardUnfolding::getUpperObjectiveBound(uint64_t objectiveIndex) { + auto& objective = this->objectives[objectiveIndex]; + if (!objective.upperResultBound) { + if (objective.formula->isProbabilityOperatorFormula()) { + objective.upperResultBound = storm::utility::one(); + } else if (objective.formula->isRewardOperatorFormula()) { + auto const& rewModel = this->model.getRewardModel(objective.formula->asRewardOperatorFormula().getRewardModelName()); + auto actionRewards = rewModel.getTotalRewardVector(this->model.getTransitionMatrix()); + if (objective.formula->getSubformula().isCumulativeRewardFormula()) { + // Try to get an upper bound by computing the maximal reward achievable within one epoch step + auto const& cumulativeRewardFormula = objective.formula->getSubformula().asCumulativeRewardFormula(); + ValueType rewardBound = cumulativeRewardFormula.template getBound(); + if (cumulativeRewardFormula.getTimeBoundReference().isRewardBound()) { + auto const& costModel = this->model.getRewardModel(cumulativeRewardFormula.getTimeBoundReference().getRewardName()); + if (!costModel.hasTransitionRewards()) { + auto actionCosts = costModel.getTotalRewardVector(this->model.getTransitionMatrix()); + ValueType largestRewardPerCost = storm::utility::zero(); + bool isFinite = true; + for (auto rewIt = actionRewards.begin(), costIt = actionCosts.begin(); rewIt != actionRewards.end(); ++rewIt, ++costIt) { + if (!storm::utility::isZero(*rewIt)) { + if (storm::utility::isZero(*costIt)) { + isFinite = false; + break; + } + ValueType rewardPerCost = *rewIt / *costIt; + largestRewardPerCost = std::max(largestRewardPerCost, rewardPerCost); + } + } + if (isFinite) { + objective.upperResultBound = largestRewardPerCost * rewardBound; + } + } + } else { + objective.upperResultBound = (*std::max_element(actionRewards.begin(), actionRewards.end())) * rewardBound; + } + + // If we could not find an upper bound, try to get an upper bound for the unbounded case + if (!objective.upperResultBound) { + storm::storage::BitVector allStates(model.getNumberOfStates(), true); + // Get the set of states from which reward is reachable + auto nonZeroRewardStates = rewModel.getStatesWithZeroReward(model.getTransitionMatrix()); + nonZeroRewardStates.complement(); + auto expRewGreater0EStates = storm::utility::graph::performProbGreater0E(model.getBackwardTransitions(), allStates, nonZeroRewardStates); + // Eliminate zero-reward ECs + auto zeroRewardChoices = rewModel.getChoicesWithZeroReward(model.getTransitionMatrix()); + auto ecElimRes = storm::transformer::EndComponentEliminator::transform(model.getTransitionMatrix(), expRewGreater0EStates, zeroRewardChoices, ~allStates); + allStates.resize(ecElimRes.matrix.getRowGroupCount()); + storm::storage::BitVector outStates(allStates.size(), false); + std::vector rew0StateProbs; + rew0StateProbs.reserve(ecElimRes.matrix.getRowCount()); + for (uint64_t state = 0; state < allStates.size(); ++ state) { + for (uint64_t choice = ecElimRes.matrix.getRowGroupIndices()[state]; choice < ecElimRes.matrix.getRowGroupIndices()[state + 1]; ++choice) { + // Check whether the choice lead to a state with expRew 0 in the original model + bool isOutChoice = false; + uint64_t originalModelChoice = ecElimRes.newToOldRowMapping[choice]; + for (auto const& entry : model.getTransitionMatrix().getRow(originalModelChoice)) { + if (!expRewGreater0EStates.get(entry.getColumn())) { + isOutChoice = true; + outStates.set(state, true); + rew0StateProbs.push_back(storm::utility::one() - ecElimRes.matrix.getRowSum(choice)); + assert (!storm::utility::isZero(rew0StateProbs.back())); + break; + } + } + if (!isOutChoice) { + rew0StateProbs.push_back(storm::utility::zero()); + } + } + } + // An upper reward bound can only be computed if it is below infinity + if (storm::utility::graph::performProb1A(ecElimRes.matrix, ecElimRes.matrix.getRowGroupIndices(), ecElimRes.matrix.transpose(true), allStates, outStates).full()) { + std::vector rewards; + rewards.reserve(ecElimRes.matrix.getRowCount()); + for (auto row : ecElimRes.newToOldRowMapping) { + rewards.push_back(actionRewards[row]); + } + storm::modelchecker::helper::BaierUpperRewardBoundsComputer baier(ecElimRes.matrix, rewards, rew0StateProbs); + objective.upperResultBound = baier.computeUpperBound(); + } + } + } + } + } + return objective.upperResultBound; + } + + template + boost::optional MultiDimensionalRewardUnfolding::getLowerObjectiveBound(uint64_t objectiveIndex) { + auto& objective = this->objectives[objectiveIndex]; + if (!objective.lowerResultBound) { + objective.lowerResultBound = storm::utility::zero(); + } + return objective.lowerResultBound; + } + + template void MultiDimensionalRewardUnfolding::setSolutionForCurrentEpoch(std::vector&& inStateSolutions) { STORM_LOG_ASSERT(currentEpoch, "Tried to set a solution for the current epoch, but no epoch was specified before."); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h index ea6c23080..11a15a82a 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -58,6 +58,12 @@ namespace storm { */ ValueType getRequiredEpochModelPrecision(Epoch const& startEpoch, ValueType const& precision); + /*! + * Returns an upper/lower bound for the objective result in every state (if this bound could be computed) + */ + boost::optional getUpperObjectiveBound(uint64_t objectiveIndex = 0); + boost::optional getLowerObjectiveBound(uint64_t objectiveIndex = 0); + void setSolutionForCurrentEpoch(std::vector&& inStateSolutions); SolutionType const& getInitialStateResult(Epoch const& epoch); // Assumes that the initial state is unique SolutionType const& getInitialStateResult(Epoch const& epoch, uint64_t initialStateIndex);