From 19925ac74d6b2dfe2e49cba85a2433442f48e3e2 Mon Sep 17 00:00:00 2001 From: TimQu Date: Fri, 30 Jun 2017 19:05:02 +0200 Subject: [PATCH] implemented value iteration based Long run average rewards for Markov automata by Butkova et al. (TACAS 2017) --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 212 ++++++++++++++---- .../helper/SparseMarkovAutomatonCslHelper.h | 27 ++- .../models/sparse/StandardRewardModel.cpp | 28 +++ src/storm/models/sparse/StandardRewardModel.h | 14 ++ 4 files changed, 231 insertions(+), 50 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index d48c38a6c..411c0fecd 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -213,7 +213,7 @@ namespace storm { } template - std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount(); @@ -228,32 +228,18 @@ namespace storm { } // Otherwise, reduce the long run average probabilities to long run average rewards. - // Every Markovian goal state s gets 1/E(s) reward for its (unique) action. - std::vector totalActionRewards(transitionMatrix.getRowCount(), storm::utility::zero()); - storm::storage::BitVector markovianGoalStates = markovianStates & psiStates; - for (auto const& state : markovianGoalStates) { - totalActionRewards[transitionMatrix.getRowGroupIndices()[state]] = storm::utility::one() / exitRateVector[state]; - } - return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, totalActionRewards, minMaxLinearEquationSolverFactory); + // Every Markovian goal state gets reward one. + std::vector stateRewards(transitionMatrix.getRowGroupCount(), storm::utility::zero()); + storm::utility::vector::setVectorValues(stateRewards, markovianStates & psiStates, storm::utility::one()); + storm::models::sparse::StandardRewardModel rewardModel(std::move(stateRewards)); + + return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, rewardModel, minMaxLinearEquationSolverFactory, useLpBasedTechnique); } template - std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { - // Obtain the total action reward vector where the state rewards are scaled accordingly - std::vector stateRewardWeights(transitionMatrix.getRowGroupCount(), storm::utility::zero()); - for (auto const markovianState : markovianStates) { - stateRewardWeights[markovianState] = storm::utility::one() / exitRateVector[markovianState]; - } - std::vector totalRewardVector = rewardModel.getTotalActionRewardVector(transitionMatrix, stateRewardWeights); - - return computeLongRunAverageRewards(dir, transitionMatrix, backwardTransitions, exitRateVector, markovianStates, totalRewardVector, minMaxLinearEquationSolverFactory); - } - - template - std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { - uint_fast64_t numberOfStates = transitionMatrix.getRowGroupCount(); // Start by decomposing the Markov automaton into its MECs. @@ -281,7 +267,7 @@ namespace storm { } // Compute the LRA value for the current MEC. - lraValuesForEndComponents.push_back(computeLraForMaximalEndComponent(dir, transitionMatrix, exitRateVector, markovianStates, totalActionRewards, mec)); + lraValuesForEndComponents.push_back(computeLraForMaximalEndComponent(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec, minMaxLinearEquationSolverFactory, useLpBasedTechnique)); } // For fast transition rewriting, we build some auxiliary data structures. @@ -410,8 +396,31 @@ namespace storm { return SparseMdpPrctlHelper::computeReachabilityRewards(dir, transitionMatrix, backwardTransitions, rewardModel, psiStates, false, false, minMaxLinearEquationSolverFactory).values; } - template - ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::storage::MaximalEndComponent const& mec) { + template + ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique) { + + // If the mec only consists of a single state, we compute the LRA value directly + if (++mec.begin() == mec.end()) { + uint_fast64_t state = mec.begin()->first; + STORM_LOG_THROW(markovianStates.get(state), storm::exceptions::InvalidOperationException, "Markov Automaton has Zeno behavior. Computation of Long Run Average values not supported."); + ValueType result = rewardModel.hasStateRewards() ? rewardModel.getStateReward(state) : storm::utility::zero(); + if (rewardModel.hasStateActionRewards() || rewardModel.hasTransitionRewards()) { + STORM_LOG_ASSERT(mec.begin()->second.size() == 1, "Markovian state has nondeterministic behavior."); + uint_fast64_t choice = *mec.begin()->second.begin(); + result += exitRateVector[state] * rewardModel.getTotalStateActionReward(state, choice, transitionMatrix, storm::utility::zero()); + } + return result; + } + + if (useLpBasedTechnique) { + return computeLraForMaximalEndComponentLP(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec); + } else { + return computeLraForMaximalEndComponentVI(dir, transitionMatrix, exitRateVector, markovianStates, rewardModel, mec, minMaxLinearEquationSolverFactory); + } + } + + template + ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec) { std::unique_ptr lpSolverFactory(new storm::utility::solver::LpSolverFactory()); std::unique_ptr solver = lpSolverFactory->create("LRA for MEC"); solver->setOptimizationDirection(invert(dir)); @@ -442,7 +451,8 @@ namespace storm { } constraint = constraint + solver->getManager().rational(storm::utility::one() / exitRateVector[state]) * k; - storm::expressions::Expression rightHandSide = solver->getManager().rational(totalActionRewards[choice]); + + storm::expressions::Expression rightHandSide = solver->getManager().rational(rewardModel.getTotalStateActionReward(state, choice, transitionMatrix, (ValueType) (storm::utility::one() / exitRateVector[state]))); if (dir == OptimizationDirection::Minimize) { constraint = constraint <= rightHandSide; } else { @@ -458,8 +468,8 @@ namespace storm { for (auto element : transitionMatrix.getRow(choice)) { constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getManager().rational(element.getValue()); } - - storm::expressions::Expression rightHandSide = solver->getManager().rational(totalActionRewards[choice]); + + storm::expressions::Expression rightHandSide = solver->getManager().rational(rewardModel.getTotalStateActionReward(state, choice, transitionMatrix, storm::utility::zero())); if (dir == OptimizationDirection::Minimize) { constraint = constraint <= rightHandSide; } else { @@ -473,8 +483,128 @@ namespace storm { solver->optimize(); return storm::utility::convertNumber(solver->getContinuousValue(k)); } - + template + ValueType SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory) { + + // Initialize data about the mec + + storm::storage::BitVector mecStates(transitionMatrix.getRowGroupCount(), false); + storm::storage::BitVector mecChoices(transitionMatrix.getRowCount(), false); + for (auto const& stateChoicesPair : mec) { + mecStates.set(stateChoicesPair.first); + for (auto const& choice : stateChoicesPair.second) { + mecChoices.set(choice); + } + } + storm::storage::BitVector markovianMecStates = mecStates & markovianStates; + storm::storage::BitVector probabilisticMecStates = mecStates & ~markovianStates; + storm::storage::BitVector probabilisticMecChoices = transitionMatrix.getRowFilter(probabilisticMecStates) & mecChoices; + STORM_LOG_THROW(!markovianMecStates.empty(), storm::exceptions::InvalidOperationException, "Markov Automaton has Zeno behavior. Computation of Long Run Average values not supported."); + + // Get the uniformization rate + + ValueType uniformizationRate = storm::utility::vector::max_if(exitRateVector, markovianMecStates); + // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop. + // Hence, we increase the uniformization rate a little. + uniformizationRate += storm::utility::one(); // Todo: try other values such as *=1.01 + + // Get the transitions of the submodel, that is + // * a matrix aMarkovian with all (uniformized) transitions from Markovian mec states to all Markovian mec states. + // * a matrix aMarkovianToProbabilistic with all (uniformized) transitions from Markovian mec states to all probabilistic mec states. + // * a matrix aProbabilistic with all transitions from probabilistic mec states to other probabilistic mec states. + // * a matrix aProbabilisticToMarkovian with all transitions from probabilistic mec states to all Markovian mec states. + typename storm::storage::SparseMatrix aMarkovian = transitionMatrix.getSubmatrix(true, markovianMecStates, markovianMecStates, true); + typename storm::storage::SparseMatrix aMarkovianToProbabilistic = transitionMatrix.getSubmatrix(true, markovianMecStates, probabilisticMecStates); + typename storm::storage::SparseMatrix aProbabilistic = transitionMatrix.getSubmatrix(false, probabilisticMecChoices, probabilisticMecStates); + typename storm::storage::SparseMatrix aProbabilisticToMarkovian = transitionMatrix.getSubmatrix(false, probabilisticMecChoices, markovianMecStates); + // The matrices with transitions from Markovian states need to be uniformized. + uint_fast64_t subState = 0; + for (auto state : markovianMecStates) { + ValueType uniformizationFactor = exitRateVector[state] / uniformizationRate; + for (auto& entry : aMarkovianToProbabilistic.getRow(subState)) { + entry.setValue(entry.getValue() * uniformizationFactor); + } + for (auto& entry : aMarkovian.getRow(subState)) { + if (entry.getColumn() == subState) { + entry.setValue(storm::utility::one() - uniformizationFactor * (storm::utility::one() - entry.getValue())); + } else { + entry.setValue(entry.getValue() * uniformizationFactor); + } + } + ++subState; + } + + // Compute the rewards obtained in a single uniformization step + + std::vector markovianChoiceRewards; + markovianChoiceRewards.reserve(aMarkovian.getRowCount()); + for (auto const& state : markovianMecStates) { + ValueType stateRewardScalingFactor = storm::utility::one() / uniformizationRate; + ValueType actionRewardScalingFactor = exitRateVector[state] / uniformizationRate; + assert(transitionMatrix.getRowGroupSize(state) == 1); + uint_fast64_t choice = transitionMatrix.getRowGroupIndices()[state]; + markovianChoiceRewards.push_back(rewardModel.getTotalStateActionReward(state, choice, transitionMatrix, stateRewardScalingFactor, actionRewardScalingFactor)); + } + + std::vector probabilisticChoiceRewards; + probabilisticChoiceRewards.reserve(aProbabilistic.getRowCount()); + for (auto const& state : probabilisticMecStates) { + uint_fast64_t groupStart = transitionMatrix.getRowGroupIndices()[state]; + uint_fast64_t groupEnd = transitionMatrix.getRowGroupIndices()[state + 1]; + for (uint_fast64_t choice = probabilisticMecChoices.getNextSetIndex(groupStart); choice < groupEnd; choice = probabilisticMecChoices.getNextSetIndex(choice + 1)) { + probabilisticChoiceRewards.push_back(rewardModel.getTotalStateActionReward(state, choice, transitionMatrix, storm::utility::zero())); + } + } + + // start the iterations + + ValueType precision = storm::utility::convertNumber(storm::settings::getModule().getPrecision()) / uniformizationRate; + std::vector v(aMarkovian.getRowCount(), storm::utility::zero()); + std::vector vOld = v; + std::vector w = v; + std::vector x(aProbabilistic.getRowGroupCount(), storm::utility::zero()); + std::vector xPrime = x; + std::vector b = probabilisticChoiceRewards; + auto solver = minMaxLinearEquationSolverFactory.create(std::move(aProbabilistic)); + solver->setCachingEnabled(true); + + while (true) { + + // Compute the expected total rewards for the probabilistic states + solver->solveEquations(dir, x, b); + + // now compute the values for the markovian states + for (uint_fast64_t row = 0; row < aMarkovian.getRowCount(); ++row) { + v[row] = markovianChoiceRewards[row] + aMarkovianToProbabilistic.multiplyRowWithVector(row, x) + aMarkovian.multiplyRowWithVector(row, w); + } + + // Check for convergence + auto vIt = v.begin(); + auto vEndIt = v.end(); + auto vOldIt = vOld.begin(); + ValueType maxDiff = *vIt - *vOldIt; + ValueType minDiff = maxDiff; + for (++vIt, ++vOldIt; vIt != vEndIt; ++vIt, ++vOldIt) { + ValueType diff = *vIt - *vOldIt; + maxDiff = std::max(maxDiff, diff); + minDiff = std::min(minDiff, diff); + } + if (maxDiff - minDiff < precision) { + break; + } + + // update the rhs of the MinMax equation system + ValueType referenceValue = v.front(); + storm::utility::vector::applyPointwise(v, w, [&referenceValue] (ValueType const& v_i) -> ValueType { return v_i - referenceValue; }); + aProbabilisticToMarkovian.multiplyWithVector(w, b); + storm::utility::vector::addVectors(b, probabilisticChoiceRewards, b); + + vOld = v; + } + return v.front() * uniformizationRate; + + } template std::vector SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair const& boundsPair, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); @@ -482,35 +612,39 @@ namespace storm { template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template void SparseMarkovAutomatonCslHelper::computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, double delta, uint_fast64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::storage::MaximalEndComponent const& mec); + template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); + + template double SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair const& boundsPair, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMarkovAutomatonCslHelper::computeUntilProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); - template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template std::vector SparseMarkovAutomatonCslHelper::computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); template std::vector SparseMarkovAutomatonCslHelper::computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); template void SparseMarkovAutomatonCslHelper::computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, storm::RationalNumber delta, uint_fast64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::storage::MaximalEndComponent const& mec); + template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique); + + template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec); + + template storm::RationalNumber SparseMarkovAutomatonCslHelper::computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::models::sparse::StandardRewardModel const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); } } diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h index 613f81d91..8faddede5 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h @@ -29,13 +29,10 @@ namespace storm { template - static std::vector computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::vector computeLongRunAverageProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); template - static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - - template - static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + static std::vector computeLongRunAverageRewards(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); template static std::vector computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); @@ -44,24 +41,32 @@ namespace storm { template ::SupportsExponential, int>::type = 0> static void computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); - template ::SupportsExponential, int>::type = 0> - static void computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + template ::SupportsExponential, int>::type = 0> + static void computeBoundedReachabilityProbabilities(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector& markovianNonGoalValues, std::vector& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); /*! * Computes the long-run average value for the given maximal end component of a Markov automaton. * + * Implementations are based on Linear Programming (LP) and Value Iteration (VI). + * * @param dir Sets whether the long-run average is to be minimized or maximized. * @param transitionMatrix The transition matrix of the underlying Markov automaton. * @param markovianStates A bit vector storing all markovian states. * @param exitRateVector A vector with exit rates for all states. Exit rates of probabilistic states are * assumed to be zero. - * @param totalActionRewards A vector indicating the total rewards obtained for each action + * @param rewardModel The considered reward model + * @param actionRewards The action rewards (earned instantaneously). * @param mec The maximal end component to consider for computing the long-run average. + * @param minMaxLinearEquationSolverFactory The factory for creating MinMaxLinearEquationSolvers (if needed by the performed method * @return The long-run average of being in a goal state for the given MEC. */ - template - static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, std::vector const& totalActionRewards, storm::storage::MaximalEndComponent const& mec); - + template + static ValueType computeLraForMaximalEndComponent(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, bool useLpBasedTechnique = false); + template + static ValueType computeLraForMaximalEndComponentLP(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec); + template + static ValueType computeLraForMaximalEndComponentVI(OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, RewardModelType const& rewardModel, storm::storage::MaximalEndComponent const& mec, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory); + }; } diff --git a/src/storm/models/sparse/StandardRewardModel.cpp b/src/storm/models/sparse/StandardRewardModel.cpp index bffd8a7e9..721608047 100644 --- a/src/storm/models/sparse/StandardRewardModel.cpp +++ b/src/storm/models/sparse/StandardRewardModel.cpp @@ -139,6 +139,30 @@ namespace storm { return StandardRewardModel(std::move(newStateRewardVector), std::move(newStateActionRewardVector), std::move(newTransitionRewardMatrix)); } + template + template + ValueType StandardRewardModel::getTotalStateActionReward(uint_fast64_t stateIndex, uint_fast64_t choiceIndex, storm::storage::SparseMatrix const& transitionMatrix, MatrixValueType const& stateRewardWeight, MatrixValueType const& actionRewardWeight) const { + ValueType result = this->hasStateRewards() ? (this->hasStateActionRewards() ? (ValueType) (this->getStateReward(stateIndex) * stateRewardWeight + this->getStateActionReward(choiceIndex) * actionRewardWeight) + : (ValueType) (this->getStateReward(stateIndex) * stateRewardWeight)) + : (this->hasStateActionRewards() ? (ValueType) (this->getStateActionReward(choiceIndex) * actionRewardWeight) + : storm::utility::zero()); + if (this->hasTransitionRewards()) { + auto rewMatrixEntryIt = this->getTransitionRewardMatrix().begin(choiceIndex); + for (auto const& transitionEntry : transitionMatrix.getRow(choiceIndex)) { + assert(rewMatrixEntryIt != this->getTransitionRewardMatrix().end(choiceIndex)); + if (transitionEntry.getColumn() < rewMatrixEntryIt->getColumn()) { + continue; + } else { + // We assume that the transition reward matrix is a submatrix of the given transition matrix. Hence, the following must hold + assert(transitionEntry.getColumn() == rewMatrixEntryIt->getColumn()); + result += actionRewardWeight * rewMatrixEntryIt->getValue() * storm::utility::convertNumber(transitionEntry.getValue()); + ++rewMatrixEntryIt; + } + } + } + return result; + } + template template void StandardRewardModel::reduceToStateBasedRewards(storm::storage::SparseMatrix const& transitionMatrix, bool reduceToStateRewards) { @@ -362,6 +386,8 @@ namespace storm { template std::vector StandardRewardModel::getTotalActionRewardVector(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& stateRewardWeights) const; template storm::storage::BitVector StandardRewardModel::getStatesWithZeroReward(storm::storage::SparseMatrix const& transitionMatrix) const; template storm::storage::BitVector StandardRewardModel::getChoicesWithZeroReward(storm::storage::SparseMatrix const& transitionMatrix) const; + template double StandardRewardModel::getTotalStateActionReward(uint_fast64_t stateIndex, uint_fast64_t choiceIndex, storm::storage::SparseMatrix const& transitionMatrix, double const& stateRewardWeight, double const& actionRewardWeight) const; + template void StandardRewardModel::reduceToStateBasedRewards(storm::storage::SparseMatrix const& transitionMatrix, bool reduceToStateRewards); template void StandardRewardModel::setStateActionReward(uint_fast64_t choiceIndex, double const & newValue); template void StandardRewardModel::setStateReward(uint_fast64_t state, double const & newValue); @@ -385,6 +411,7 @@ namespace storm { template std::vector StandardRewardModel::getTotalActionRewardVector(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& stateRewardWeights) const; template storm::storage::BitVector StandardRewardModel::getStatesWithZeroReward(storm::storage::SparseMatrix const& transitionMatrix) const; template storm::storage::BitVector StandardRewardModel::getChoicesWithZeroReward(storm::storage::SparseMatrix const& transitionMatrix) const; + template storm::RationalNumber StandardRewardModel::getTotalStateActionReward(uint_fast64_t stateIndex, uint_fast64_t choiceIndex, storm::storage::SparseMatrix const& transitionMatrix, storm::RationalNumber const& stateRewardWeight, storm::RationalNumber const& actionRewardWeight) const; template void StandardRewardModel::reduceToStateBasedRewards(storm::storage::SparseMatrix const& transitionMatrix, bool reduceToStateRewards); template void StandardRewardModel::setStateActionReward(uint_fast64_t choiceIndex, storm::RationalNumber const & newValue); template void StandardRewardModel::setStateReward(uint_fast64_t state, storm::RationalNumber const & newValue); @@ -398,6 +425,7 @@ namespace storm { template storm::storage::BitVector StandardRewardModel::getChoicesWithZeroReward(storm::storage::SparseMatrix const& transitionMatrix) const; template std::vector StandardRewardModel::getTotalActionRewardVector(storm::storage::SparseMatrix const& transitionMatrix, std::vector const& stateRewardWeights) const; + template storm::RationalFunction StandardRewardModel::getTotalStateActionReward(uint_fast64_t stateIndex, uint_fast64_t choiceIndex, storm::storage::SparseMatrix const& transitionMatrix, storm::RationalFunction const& stateRewardWeight, storm::RationalFunction const& actionRewardWeight) const; template void StandardRewardModel::reduceToStateBasedRewards(storm::storage::SparseMatrix const& transitionMatrix, bool reduceToStateRewards); template void StandardRewardModel::setStateActionReward(uint_fast64_t choiceIndex, storm::RationalFunction const & newValue); template void StandardRewardModel::setStateReward(uint_fast64_t state, storm::RationalFunction const & newValue); diff --git a/src/storm/models/sparse/StandardRewardModel.h b/src/storm/models/sparse/StandardRewardModel.h index 2a3033ea7..c58f38312 100644 --- a/src/storm/models/sparse/StandardRewardModel.h +++ b/src/storm/models/sparse/StandardRewardModel.h @@ -4,6 +4,7 @@ #include #include "storm/storage/SparseMatrix.h" +#include "storm/utility/constants.h" #include "storm/utility/OsDetection.h" #include "storm/adapters/RationalFunctionAdapter.h" @@ -161,6 +162,19 @@ namespace storm { * @return The transition reward matrix if there is one. */ boost::optional> const& getOptionalTransitionRewardMatrix() const; + + /*! + * Retrieves the total reward for the given state action pair (including (scaled) state rewards, action rewards and transition rewards + * + * @param stateIndex The index of the considered state + * @param choiceIndex The index of the considered choice + * @param transitionMatrix The matrix that is used to weight the values of the transition reward matrix. + * @param stateRewardWeight The weight that is used to scale the state reward + * @param actionRewardWeight The weight that is used to scale the action based rewards (includes action and transition rewards). + * @return + */ + template + ValueType getTotalStateActionReward(uint_fast64_t stateIndex, uint_fast64_t choiceIndex, storm::storage::SparseMatrix const& transitionMatrix, MatrixValueType const& stateRewardWeight = storm::utility::one(), MatrixValueType const& actionRewardWeight = storm::utility::one()) const; /*! * Creates a new reward model by restricting the actions of the action-based rewards.