From 2e593dc014084b6071db4a3ac8be0abbb2bf8794 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Fri, 11 Dec 2020 17:08:55 +0100 Subject: [PATCH] Added computation of steady state probabilities for DTMC/CTMC in the sparse engine. --- CHANGELOG.md | 1 + src/storm-cli-utilities/model-handling.h | 12 ++ src/storm/api/verification.h | 27 +++ .../csl/SparseCtmcCslModelChecker.cpp | 23 +++ .../csl/SparseCtmcCslModelChecker.h | 6 + ...arseDeterministicInfiniteHorizonHelper.cpp | 172 ++++++++++++++++-- ...SparseDeterministicInfiniteHorizonHelper.h | 29 ++- .../prctl/SparseDtmcPrctlModelChecker.cpp | 22 +++ .../prctl/SparseDtmcPrctlModelChecker.h | 7 + .../helper/BaierUpperRewardBoundsComputer.h | 2 +- .../helper/DsMpiUpperRewardBoundsComputer.h | 4 +- src/storm/settings/modules/IOSettings.cpp | 6 + src/storm/settings/modules/IOSettings.h | 8 +- 13 files changed, 299 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c9f74db7..ed86c3df0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ Version 1.6.x ------------- ## Version 1.6.4 (20xx/xx) - Added an export of check results to json. Use `--exportresult` in the command line interface. +- Added computation of steady state probabilities for DTMC/CTMC in the sparse engine. Use `--steadystate` in the command line interface. ## Version 1.6.3 (2020/11) - Added support for multi-objective model checking of long-run average objectives including mixtures with other kinds of objectives. diff --git a/src/storm-cli-utilities/model-handling.h b/src/storm-cli-utilities/model-handling.h index fda64be44..b5103cae8 100644 --- a/src/storm-cli-utilities/model-handling.h +++ b/src/storm-cli-utilities/model-handling.h @@ -1043,6 +1043,18 @@ namespace storm { ++exportCount; }; verifyProperties(input,verificationCallback, postprocessingCallback); + if (ioSettings.isComputeSteadyStateDistributionSet()) { + storm::utility::Stopwatch watch(true); + std::unique_ptr result; + try { + result = storm::api::computeSteadyStateDistributionWithSparseEngine(mpi.env, sparseModel); + } catch (storm::exceptions::BaseException const& ex) { + STORM_LOG_WARN("Cannot compute steady-state probabilities: " << ex.what()); + } + watch.stop(); + postprocessingCallback(result); + STORM_PRINT((storm::utility::resources::isTerminate() ? "Result till abort:" : "Result:") << *result); + } } template diff --git a/src/storm/api/verification.h b/src/storm/api/verification.h index 3e5f916a5..f9cbda29f 100644 --- a/src/storm/api/verification.h +++ b/src/storm/api/verification.h @@ -273,7 +273,34 @@ namespace storm { Environment env; return verifyWithSparseEngine(env, model, task); } + + template + std::unique_ptr computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr> const& dtmc) { + std::unique_ptr result; + storm::modelchecker::SparseDtmcPrctlModelChecker> modelchecker(*dtmc); + return modelchecker.computeSteadyStateDistribution(env); + } + + template + std::unique_ptr computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr> const& ctmc) { + std::unique_ptr result; + storm::modelchecker::SparseCtmcCslModelChecker> modelchecker(*ctmc); + return modelchecker.computeSteadyStateDistribution(env); + } + template + std::unique_ptr computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr> const& model) { + std::unique_ptr result; + if (model->getType() == storm::models::ModelType::Dtmc) { + result = computeSteadyStateDistributionWithSparseEngine(env, model->template as>()); + } else if (model->getType() == storm::models::ModelType::Ctmc) { + result = computeSteadyStateDistributionWithSparseEngine(env, model->template as>()); + } else { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing the long run average distribution for the model type " << model->getType() << " is not supported."); + } + return result; + } + // // Verifying with Hybrid engine // diff --git a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 11530673b..513f9d8cb 100644 --- a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -172,6 +172,29 @@ namespace storm { std::vector result = storm::modelchecker::helper::SparseCtmcCslHelper::computeAllTransientProbabilities(env, this->getModel().getTransitionMatrix(), this->getModel().getInitialStates(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), upperBound); return result; } + + template + std::unique_ptr SparseCtmcCslModelChecker::computeSteadyStateDistribution(Environment const& env) { + // Initialize helper + auto probabilisticTransitions = this->getModel().computeProbabilityMatrix(); + storm::modelchecker::helper::SparseDeterministicInfiniteHorizonHelper helper(probabilisticTransitions, this->getModel().getExitRateVector()); + + // Compute result + std::vector result; + auto const& initialStates = this->getModel().getInitialStates(); + uint64_t numInitStates = initialStates.getNumberOfSetBits(); + if (numInitStates == 1) { + result = helper.computeLongRunAverageStateDistribution(env, *initialStates.begin()); + } else { + STORM_LOG_WARN("Multiple initial states found. A uniform distribution over initial states is assumed."); + ValueType initProb = storm::utility::one() / storm::utility::convertNumber(numInitStates); + result = helper.computeLongRunAverageStateDistribution(env, [&initialStates,&initProb](uint64_t const& stateIndex) { return initialStates.get(stateIndex) ? initProb : storm::utility::zero(); }); + } + + // Return CheckResult + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + } + // Explicitly instantiate the model checker. template class SparseCtmcCslModelChecker>; diff --git a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.h index bdb896a7b..fa43b6f9e 100644 --- a/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/storm/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -43,6 +43,12 @@ namespace storm { */ std::vector computeAllTransientProbabilities(Environment const& env, CheckTask const& checkTask); + /*! + * Computes the long run average (or: steady state) distribution over all states + * Assumes a uniform distribution over initial states. + */ + std::unique_ptr computeSteadyStateDistribution(Environment const& env); + }; } // namespace modelchecker diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp index 23ceaf9cb..a1db02872 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp @@ -2,15 +2,13 @@ #include "storm/modelchecker/helper/infinitehorizon/internal/ComponentUtility.h" #include "storm/modelchecker/helper/infinitehorizon/internal/LraViHelper.h" - +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" #include "storm/storage/SparseMatrix.h" #include "storm/storage/StronglyConnectedComponentDecomposition.h" #include "storm/storage/Scheduler.h" #include "storm/solver/LinearEquationSolver.h" -#include "storm/solver/Multiplier.h" -#include "storm/solver/LpSolver.h" #include "storm/utility/SignalHandler.h" #include "storm/utility/solver.h" @@ -71,7 +69,7 @@ namespace storm { return computeLraForBsccVi(env, stateValueGetter, actionValueGetter, component); } else if (method == storm::solver::LraMethod::LraDistributionEquations) { // We only need the first element of the pair as the lra distribution is not relevant at this point. - return computeLraForBsccLraDistr(env, stateValueGetter, actionValueGetter, component).first; + return computeLraForBsccSteadyStateDistr(env, stateValueGetter, actionValueGetter, component).first; } STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to " << storm::solver::toString(storm::solver::LraMethod::GainBiasEquations) << "."); // We don't need the bias values @@ -124,6 +122,9 @@ namespace storm { template std::pair> SparseDeterministicInfiniteHorizonHelper::computeLraForBsccGainBias(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) { + // we want that the returned vector is sorted as the bscc. So let's assert that the bscc is sorted ascendingly. + STORM_LOG_ASSERT(std::is_sorted(bscc.begin(), bscc.end()), "Expected that bsccs are sorted."); + // We build the equation system as in Line 3 of Algorithm 3 from // Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017) // https://doi.org/10.1007/978-3-319-68167-2_25 @@ -212,17 +213,19 @@ namespace storm { } template - std::pair> SparseDeterministicInfiniteHorizonHelper::computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) { - + std::vector SparseDeterministicInfiniteHorizonHelper::computeSteadyStateDistrForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc) { + // We want that the returned values are sorted properly. Let's assert that strongly connected components use a sorted container type + STORM_LOG_ASSERT(std::is_sorted(bscc.begin(), bscc.end()), "Expected that bsccs are sorted."); + // Let A be ab auxiliary Matrix with A[s,s] = R(s,s) - r(s) & A[s,s'] = R(s,s') for s,s' in BSCC and s!=s'. // We build and solve the equation system for // x*A=0 & x_0+...+x_n=1 <=> A^t*x=0=x-x & x_0+...+x_n=1 <=> (1+A^t)*x = x & 1-x_0-...-x_n-1=x_n // Then, x[i] will be the fraction of the time we are in state i. - // This method assumes that this BSCC consist of more than one state + // This method assumes that this BSCC consist of more than one state. + // We therefore catch the (easy) case where the BSCC is a singleton. if (bscc.size() == 1) { - ValueType lraValue = stateValuesGetter(*bscc.begin()) + (this->isContinuousTime() ? (*this->_exitRates)[*bscc.begin()] * actionValuesGetter(*bscc.begin()) : actionValuesGetter(*bscc.begin())); - return { lraValue, {storm::utility::one()} }; + return { storm::utility::one() }; } // Prepare an environment for the underlying linear equation solver @@ -318,12 +321,28 @@ namespace storm { requirements.clearUpperBounds(); STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); - std::vector lraDistr(bscc.size(), storm::utility::one() / storm::utility::convertNumber(bscc.size())); - solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide); + std::vector steadyStateDistr(bscc.size(), storm::utility::one() / storm::utility::convertNumber(bscc.size())); + solver->solveEquations(subEnv, steadyStateDistr, bsccEquationSystemRightSide); + + // As a last step, we normalize these values to counter numerical inaccuracies a bit. + // This is only reasonable in non-exact mode. + if (!subEnv.solver().isForceExact()) { + ValueType sum = std::accumulate(steadyStateDistr.begin(), steadyStateDistr.end(), storm::utility::zero()); + storm::utility::vector::scaleVectorInPlace(steadyStateDistr, storm::utility::one() / sum); + } + + return steadyStateDistr; + } + + template + std::pair> SparseDeterministicInfiniteHorizonHelper::computeLraForBsccSteadyStateDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) { + + // Get the LRA distribution for the provided bscc. + auto steadyStateDistr = computeSteadyStateDistrForBscc(env, bscc); // Calculate final LRA Value ValueType result = storm::utility::zero(); - auto solIt = lraDistr.begin(); + auto solIt = steadyStateDistr.begin(); for (auto const& globalState : bscc) { if (this->isContinuousTime()) { result += (*solIt) * (stateValuesGetter(globalState) + (*this->_exitRates)[globalState] * actionValuesGetter(globalState)); @@ -332,9 +351,9 @@ namespace storm { } ++solIt; } - assert(solIt == lraDistr.end()); + assert(solIt == steadyStateDistr.end()); - return std::pair>(std::move(result), std::move(lraDistr)); + return std::pair>(std::move(result), std::move(steadyStateDistr)); } template @@ -423,6 +442,131 @@ namespace storm { } return result; } + + template + std::vector SparseDeterministicInfiniteHorizonHelper::computeLongRunAverageStateDistribution(Environment const& env) { + createDecomposition(); + STORM_LOG_THROW(this->_longRunComponentDecomposition->size() <= 1, storm::exceptions::InvalidOperationException, ""); + return computeLongRunAverageStateDistribution(env, [] (uint64_t) { return storm::utility::zero();} ); + } + + template + std::vector SparseDeterministicInfiniteHorizonHelper::computeLongRunAverageStateDistribution(Environment const& env, uint64_t const& initialState) { + STORM_LOG_ASSERT(initialState < this->_transitionMatrix.getRowGroupCount(), "Invlid initial state index: " << initialState << ". Have only " << this->_transitionMatrix.getRowGroupCount() << " states."); + return computeLongRunAverageStateDistribution(env, [&initialState] (uint64_t stateIndex ) { return initialState == stateIndex ? storm::utility::one() : storm::utility::zero();} ); + } + + template + std::vector SparseDeterministicInfiniteHorizonHelper::computeLongRunAverageStateDistribution(Environment const& env, ValueGetter const& initialDistributionGetter) { + createDecomposition(); + + // Compute for each BSCC get the probability with which we reach that BSCC + auto bsccReachProbs = computeBsccReachabilityProbabilities(env, initialDistributionGetter); + // We are now ready to compute the resulting lra distribution + std::vector steadyStateDistr(this->_transitionMatrix.getRowGroupCount(), storm::utility::zero()); + for (uint64_t currentComponentIndex = 0; currentComponentIndex < this->_longRunComponentDecomposition->size(); ++currentComponentIndex) { + auto const& component = (*this->_longRunComponentDecomposition)[currentComponentIndex]; + // Compute distribution for current bscc + auto bsccDistr = this->computeSteadyStateDistrForBscc(env, component); + // Scale with probability to reach that bscc + auto const& scalingFactor = bsccReachProbs[currentComponentIndex]; + if (!storm::utility::isOne(scalingFactor)) { + storm::utility::vector::scaleVectorInPlace(bsccDistr, scalingFactor); + } + // Set the values in the result vector + auto bsccDistrIt = bsccDistr.begin(); + for (auto const& element : component) { + uint64_t state = internal::getComponentElementState(element); + steadyStateDistr[state] = *bsccDistrIt; + ++bsccDistrIt; + } + STORM_LOG_ASSERT(bsccDistrIt == bsccDistr.end(), "Unexpected number of entries in bscc distribution"); + } + return steadyStateDistr; + } + + template + std::vector computeUpperBoundsForExpectedVisitingTimes(storm::storage::SparseMatrix const& nonBsccMatrix, std::vector const& toBsccProbabilities) { + return storm::modelchecker::helper::BaierUpperRewardBoundsComputer::computeUpperBoundOnExpectedVisitingTimes(nonBsccMatrix, toBsccProbabilities); + } + + template <> + std::vector computeUpperBoundsForExpectedVisitingTimes(storm::storage::SparseMatrix const& nonBsccMatrix, std::vector const& toBsccProbabilities) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing upper bounds for expected visiting times over rational functions is not supported."); + } + + template + std::vector SparseDeterministicInfiniteHorizonHelper::computeBsccReachabilityProbabilities(Environment const& env, ValueGetter const& initialDistributionGetter) { + STORM_LOG_ASSERT(this->_longRunComponentDecomposition != nullptr, "Decomposition not computed, yet."); + uint64_t numBSCCs = this->_longRunComponentDecomposition->size(); + STORM_LOG_ASSERT(numBSCCs > 0, "Found 0 BSCCs in a Markov chain. This should not be possible."); + + // Compute for each BSCC get the probability with which we reach that BSCC + std::vector bsccReachProbs(numBSCCs, storm::utility::zero()); + if (numBSCCs == 1) { + bsccReachProbs.front() = storm::utility::one(); + } else { + // Compute mapping from state indices to their BSCC and + // the set of states that do not lie in a BSCC + std::vector stateToBsccMap(this->_transitionMatrix.getRowGroupCount(), std::numeric_limits::max()); + storm::storage::BitVector statesNotInComponent(this->_transitionMatrix.getRowGroupCount(), true); + for (uint64_t currentComponentIndex = 0; currentComponentIndex < this->_longRunComponentDecomposition->size(); ++currentComponentIndex) { + for (auto const& element : (*this->_longRunComponentDecomposition)[currentComponentIndex]) { + uint64_t state = internal::getComponentElementState(element); + statesNotInComponent.set(state, false); + stateToBsccMap[state] = currentComponentIndex; + } + } + // For each non-BSCC state, compute the expected number of times we visit that state by solving an equation system + // The equation system basically considers an equation `init(s) + sum_s' P(s',s) * x_s' = x_s` for each state s. + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + bool isEqSysFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + auto eqSysMatrix = this->_transitionMatrix.getSubmatrix(false, statesNotInComponent, statesNotInComponent, isEqSysFormat); + std::unique_ptr> solver = linearEquationSolverFactory.create(env, eqSysMatrix); + solver->setLowerBound(storm::utility::zero()); + // Check solver requirements + auto requirements = solver->getRequirements(env); + requirements.clearLowerBounds(); + if (requirements.upperBounds()) { + auto toBsccProbabilities = this->_transitionMatrix.getConstrainedRowSumVector(statesNotInComponent, ~statesNotInComponent); + solver->setUpperBounds(computeUpperBoundsForExpectedVisitingTimes(eqSysMatrix, toBsccProbabilities)); + requirements.clearUpperBounds(); + } + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + eqSysMatrix.transpose(false, true); // Transpose so that each row considers the predecessors + if (isEqSysFormat) { + eqSysMatrix.convertToEquationSystem(); + } + std::vector eqSysRhs; + eqSysRhs.reserve(eqSysMatrix.getRowCount()); + for (auto state : statesNotInComponent) { + eqSysRhs.push_back(initialDistributionGetter(state)); + } + std::vector eqSysValues(eqSysMatrix.getRowCount()); + solver->solveEquations(env, eqSysValues, eqSysRhs); + + // Now that we have the expected number of times we visit each non-BSCC state, we can compute the probabilities to reach each bscc. + // Run over all transitions that enter a BSCC + uint64_t eqSysStateIndex = 0; + for (auto state : statesNotInComponent) { + for (auto const& entry : this->_transitionMatrix.getRow(state)) { + if (!statesNotInComponent.get(entry.getColumn())) { + // Add the probability that the BSCC is reached from this state. + bsccReachProbs[stateToBsccMap[entry.getColumn()]] += entry.getValue() * eqSysValues[eqSysStateIndex]; + } + } + ++eqSysStateIndex; + } + } + // As a last step, we normalize these values to counter numerical inaccuracies a bit. + // This is only reasonable in non-exact mode. + if (!env.solver().isForceExact()) { + ValueType sum = std::accumulate(bsccReachProbs.begin(), bsccReachProbs.end(), storm::utility::zero()); + storm::utility::vector::scaleVectorInPlace(bsccReachProbs, storm::utility::one() / sum); + } + return bsccReachProbs; + } + template class SparseDeterministicInfiniteHorizonHelper; template class SparseDeterministicInfiniteHorizonHelper; diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h index b5a78fad1..800adddb2 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h @@ -39,10 +39,34 @@ namespace storm { */ virtual ValueType computeLraForComponent(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& component) override; + /*! + * Computes the long run average state distribution, i.e., a vector that assigns for each state s the average fraction of the time we spent in s, + * assuming an infinite run. + * Note that all the probability muss will be in bottom SCCs. + * If there are multiple bottom SCCs, the value will depend on the initial state. + * It is also possible to provide a distribution over initial states. + * However, the standard treatment of multiple initial states is not possible here, as this is not defined then. + * @throws InvalidOperationException if no initial state or distribution is provided and the model has multiple BSCCs. + */ + std::vector computeLongRunAverageStateDistribution(Environment const& env); + std::vector computeLongRunAverageStateDistribution(Environment const& env, uint64_t const& initialState); + std::vector computeLongRunAverageStateDistribution(Environment const& env, ValueGetter const& initialDistributionGetter); + protected: virtual void createDecomposition() override; + /*! + * Computes for each BSCC the probability to reach that SCC assuming the given distribution over initial states. + */ + std::vector computeBsccReachabilityProbabilities(Environment const& env, ValueGetter const& initialDistributionGetter); + + /*! + * Computes the long run average (steady state) distribution for the given BSCC. + */ + std::vector computeSteadyStateDistrForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc); + + std::pair computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); /*! @@ -55,10 +79,11 @@ namespace storm { * @see Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017), https://doi.org/10.1007/978-3-319-68167-2_25 */ std::pair> computeLraForBsccGainBias(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + /*! - * As computeLraForComponent but solves a linear equation system consisting encoding the long run average (steady state) distribution (independent of what is set in env) + * As computeLraForComponent but does the computation by computing the long run average (steady state) distribution (independent of what is set in env) */ - std::pair> computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + std::pair> computeLraForBsccSteadyStateDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); std::pair, std::vector> buildSspMatrixVector(std::vector const& bsccLraValues, std::vector const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem); diff --git a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 831c3ac96..30858ae3a 100644 --- a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -244,6 +244,28 @@ namespace storm { } } + template + std::unique_ptr SparseDtmcPrctlModelChecker::computeSteadyStateDistribution(Environment const& env) { + // Initialize helper + storm::modelchecker::helper::SparseDeterministicInfiniteHorizonHelper helper(this->getModel().getTransitionMatrix()); + + // Compute result + std::vector result; + auto const& initialStates = this->getModel().getInitialStates(); + uint64_t numInitStates = initialStates.getNumberOfSetBits(); + if (numInitStates == 1) { + result = helper.computeLongRunAverageStateDistribution(env, *initialStates.begin()); + } else { + STORM_LOG_WARN("Multiple initial states found. A uniform distribution over initial states is assumed."); + ValueType initProb = storm::utility::one() / storm::utility::convertNumber(numInitStates); + result = helper.computeLongRunAverageStateDistribution(env, [&initialStates,&initProb](uint64_t const& stateIndex) { return initialStates.get(stateIndex) ? initProb : storm::utility::zero(); }); + } + + // Return CheckResult + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(result))); + } + + template class SparseDtmcPrctlModelChecker>; #ifdef STORM_HAVE_CARL diff --git a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index 8e72833d6..93ae8dbfb 100644 --- a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -41,6 +41,13 @@ namespace storm { virtual std::unique_ptr computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr computeReachabilityTimes(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask const& checkTask) override; virtual std::unique_ptr checkQuantileFormula(Environment const& env, CheckTask const& checkTask) override; + + /*! + * Computes the long run average (or: steady state) distribution over all states + * Assumes a uniform distribution over initial states. + */ + std::unique_ptr computeSteadyStateDistribution(Environment const& env); + }; } // namespace modelchecker diff --git a/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h b/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h index cf47c4c5a..864414c23 100644 --- a/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h +++ b/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h @@ -16,7 +16,7 @@ namespace storm { public: /*! * Creates an object that can compute upper bounds on the *maximal* expected rewards for the provided MDP. - * + * @see http://doi.org/10.1007/978-3-319-63387-9_8 * @param transitionMatrix The matrix defining the transitions of the system without the transitions * that lead directly to the goal state. * @param rewards The rewards of each choice. diff --git a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h index f49470dcf..d691ecb18 100644 --- a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h +++ b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h @@ -20,7 +20,7 @@ namespace storm { public: /*! * Creates an object that can compute upper bounds on the expected rewards for the provided DTMC. - * + * @see http://doi.org/10.1145/1102351.1102423 * @param transitionMatrix The matrix defining the transitions of the system without the transitions * that lead directly to the goal state. * @param rewards The rewards of each state. @@ -81,7 +81,7 @@ namespace storm { public: /*! * Creates an object that can compute upper bounds on the *minimal* expected rewards for the provided MDP. - * + * @see http://doi.org/10.1145/1102351.1102423 * @param transitionMatrix The matrix defining the transitions of the system without the transitions * that lead directly to the goal state. * @param rewards The rewards of each choice. diff --git a/src/storm/settings/modules/IOSettings.cpp b/src/storm/settings/modules/IOSettings.cpp index 82d2de9c8..40c1a006c 100644 --- a/src/storm/settings/modules/IOSettings.cpp +++ b/src/storm/settings/modules/IOSettings.cpp @@ -47,6 +47,7 @@ namespace storm { const std::string IOSettings::janiPropertyOptionShortName = "jprop"; const std::string IOSettings::propertyOptionName = "prop"; const std::string IOSettings::propertyOptionShortName = "prop"; + const std::string IOSettings::steadyStateDistrOptionName = "steadystate"; const std::string IOSettings::qvbsInputOptionName = "qvbs"; const std::string IOSettings::qvbsInputOptionShortName = "qvbs"; @@ -100,6 +101,7 @@ namespace storm { .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of constants and their value, e.g. a=1,b=2,c=3.").setDefaultValueString("").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, janiPropertyOptionName, false, "Specifies the properties from the jani model (given by --" + janiInputOptionName + ") to be checked.").setShortName(janiPropertyOptionShortName) .addArgument(storm::settings::ArgumentBuilder::createStringArgument("values", "A comma separated list of properties to be checked").setDefaultValueString("").makeOptional().build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, steadyStateDistrOptionName, false, "Computes the steady state distribution. Result can be exported using --" + exportCheckResultOptionName +".").setIsAdvanced().build()); this->addOption(storm::settings::OptionBuilder(moduleName, qvbsInputOptionName, false, "Selects a model from the Quantitative Verification Benchmark Set.").setShortName(qvbsInputOptionShortName) .addArgument(storm::settings::ArgumentBuilder::createStringArgument("model", "The short model name as in the benchmark set.").build()) @@ -302,6 +304,10 @@ namespace storm { return this->getOption(propertyOptionName).getArgumentByName("filter").getValueAsString(); } + bool IOSettings::isComputeSteadyStateDistributionSet() const { + return this->getOption(steadyStateDistrOptionName).getHasOptionBeenSet(); + } + bool IOSettings::isQvbsInputSet() const { return this->getOption(qvbsInputOptionName).getHasOptionBeenSet(); } diff --git a/src/storm/settings/modules/IOSettings.h b/src/storm/settings/modules/IOSettings.h index 3c7f0bb44..34875441f 100644 --- a/src/storm/settings/modules/IOSettings.h +++ b/src/storm/settings/modules/IOSettings.h @@ -323,7 +323,12 @@ namespace storm { * @return The property filter. */ std::string getPropertyFilter() const; - + + /*! + * Retrieves whether the steady-state distribution is to be computed. + */ + bool isComputeSteadyStateDistributionSet() const; + /*! * Retrieves whether the input model is to be read from the quantitative verification benchmark set (QVBS) */ @@ -390,6 +395,7 @@ namespace storm { static const std::string janiPropertyOptionShortName; static const std::string propertyOptionName; static const std::string propertyOptionShortName; + static const std::string steadyStateDistrOptionName; static const std::string qvbsInputOptionName; static const std::string qvbsInputOptionShortName; static const std::string qvbsRootOptionName;