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<ValueType>(input,verificationCallback, postprocessingCallback);
+            if (ioSettings.isComputeSteadyStateDistributionSet()) {
+                storm::utility::Stopwatch watch(true);
+                std::unique_ptr<storm::modelchecker::CheckResult> result;
+                try {
+                    result = storm::api::computeSteadyStateDistributionWithSparseEngine<ValueType>(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 <storm::dd::DdType DdType, typename ValueType>
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<typename ValueType>
+        std::unique_ptr<storm::modelchecker::CheckResult> computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr<storm::models::sparse::Dtmc<ValueType>> const& dtmc) {
+            std::unique_ptr<storm::modelchecker::CheckResult> result;
+            storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<ValueType>> modelchecker(*dtmc);
+            return modelchecker.computeSteadyStateDistribution(env);
+        }
+
+        template<typename ValueType>
+        std::unique_ptr<storm::modelchecker::CheckResult> computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr<storm::models::sparse::Ctmc<ValueType>> const& ctmc) {
+            std::unique_ptr<storm::modelchecker::CheckResult> result;
+            storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<ValueType>> modelchecker(*ctmc);
+            return modelchecker.computeSteadyStateDistribution(env);
+        }
 
+        template<typename ValueType>
+        std::unique_ptr<storm::modelchecker::CheckResult> computeSteadyStateDistributionWithSparseEngine(storm::Environment const& env, std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model) {
+            std::unique_ptr<storm::modelchecker::CheckResult> result;
+            if (model->getType() == storm::models::ModelType::Dtmc) {
+                result = computeSteadyStateDistributionWithSparseEngine(env, model->template as<storm::models::sparse::Dtmc<ValueType>>());
+            } else if (model->getType() == storm::models::ModelType::Ctmc) {
+                result = computeSteadyStateDistributionWithSparseEngine(env, model->template as<storm::models::sparse::Ctmc<ValueType>>());
+            } 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<ValueType> result = storm::modelchecker::helper::SparseCtmcCslHelper::computeAllTransientProbabilities(env, this->getModel().getTransitionMatrix(), this->getModel().getInitialStates(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), upperBound);
             return result;
         }
+        
+        template <typename SparseCtmcModelType>
+        std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<SparseCtmcModelType>::computeSteadyStateDistribution(Environment const& env) {
+            // Initialize helper
+            auto probabilisticTransitions = this->getModel().computeProbabilityMatrix();
+            storm::modelchecker::helper::SparseDeterministicInfiniteHorizonHelper<ValueType> helper(probabilisticTransitions, this->getModel().getExitRateVector());
+            
+            // Compute result
+            std::vector<ValueType> 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<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(numInitStates);
+                result = helper.computeLongRunAverageStateDistribution(env, [&initialStates,&initProb](uint64_t const& stateIndex) { return initialStates.get(stateIndex) ? initProb : storm::utility::zero<ValueType>(); });
+            }
+
+            // Return CheckResult
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(result)));
+        }
+
 
         // Explicitly instantiate the model checker.
         template class SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>>;
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<ValueType> computeAllTransientProbabilities(Environment const& env, CheckTask<storm::logic::BoundedUntilFormula, ValueType> const& checkTask);
 
+            /*!
+             * Computes the long run average (or: steady state) distribution over all states
+             * Assumes a uniform distribution over initial states.
+             */
+            std::unique_ptr<CheckResult> 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 <typename ValueType>
             std::pair<ValueType, std::vector<ValueType>> SparseDeterministicInfiniteHorizonHelper<ValueType>::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 <typename ValueType>
-            std::pair<ValueType, std::vector<ValueType>> SparseDeterministicInfiniteHorizonHelper<ValueType>::computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) {
-
+            std::vector<ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>()} };
+                    return { storm::utility::one<ValueType>() };
                 }
                 
                 // 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<ValueType> lraDistr(bscc.size(), storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size()));
-                solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide);
+                std::vector<ValueType> steadyStateDistr(bscc.size(), storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(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<ValueType>());
+                    storm::utility::vector::scaleVectorInPlace<ValueType, ValueType>(steadyStateDistr, storm::utility::one<ValueType>() / sum);
+                }
+                
+                return steadyStateDistr;
+            }
+            
+            template <typename ValueType>
+            std::pair<ValueType, std::vector<ValueType>> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>();
-                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<ValueType, std::vector<ValueType>>(std::move(result), std::move(lraDistr));
+                return std::pair<ValueType, std::vector<ValueType>>(std::move(result), std::move(steadyStateDistr));
             }
 
             template <typename ValueType>
@@ -423,6 +442,131 @@ namespace storm {
                 }
                 return result;
             }
+
+            template <typename ValueType>
+            std::vector<ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::computeLongRunAverageStateDistribution(Environment const& env) {
+                createDecomposition();
+                STORM_LOG_THROW(this->_longRunComponentDecomposition->size() <= 1, storm::exceptions::InvalidOperationException, "");
+                return computeLongRunAverageStateDistribution(env, [] (uint64_t) { return storm::utility::zero<ValueType>();} );
+            }
+            
+            template <typename ValueType>
+            std::vector<ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>() : storm::utility::zero<ValueType>();} );
+            }
+            
+            template <typename ValueType>
+            std::vector<ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType> steadyStateDistr(this->_transitionMatrix.getRowGroupCount(), storm::utility::zero<ValueType>());
+                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 <typename ValueType>
+            std::vector<ValueType> computeUpperBoundsForExpectedVisitingTimes(storm::storage::SparseMatrix<ValueType> const& nonBsccMatrix, std::vector<ValueType> const& toBsccProbabilities) {
+                return storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ValueType>::computeUpperBoundOnExpectedVisitingTimes(nonBsccMatrix, toBsccProbabilities);
+            }
+            
+            template <>
+            std::vector<storm::RationalFunction> computeUpperBoundsForExpectedVisitingTimes(storm::storage::SparseMatrix<storm::RationalFunction> const& nonBsccMatrix, std::vector<storm::RationalFunction> const& toBsccProbabilities) {
+                STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Computing upper bounds for expected visiting times over rational functions is not supported.");
+            }
+            
+            template <typename ValueType>
+            std::vector<ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType> bsccReachProbs(numBSCCs, storm::utility::zero<ValueType>());
+                if (numBSCCs == 1) {
+                    bsccReachProbs.front() = storm::utility::one<ValueType>();
+                } else {
+                    // Compute mapping from state indices to their BSCC and
+                    // the set of states that do not lie in a BSCC
+                    std::vector<uint64_t> stateToBsccMap(this->_transitionMatrix.getRowGroupCount(), std::numeric_limits<uint64_t>::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<ValueType> linearEquationSolverFactory;
+                    bool isEqSysFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
+                    auto eqSysMatrix = this->_transitionMatrix.getSubmatrix(false, statesNotInComponent, statesNotInComponent, isEqSysFormat);
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(env, eqSysMatrix);
+                    solver->setLowerBound(storm::utility::zero<ValueType>());
+                    // 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<ValueType> eqSysRhs;
+                    eqSysRhs.reserve(eqSysMatrix.getRowCount());
+                    for (auto state : statesNotInComponent) {
+                        eqSysRhs.push_back(initialDistributionGetter(state));
+                    }
+                    std::vector<ValueType> 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<ValueType>());
+                    storm::utility::vector::scaleVectorInPlace<ValueType, ValueType>(bsccReachProbs, storm::utility::one<ValueType>() / sum);
+                }
+                return bsccReachProbs;
+            }
+
             
             template class SparseDeterministicInfiniteHorizonHelper<double>;
             template class SparseDeterministicInfiniteHorizonHelper<storm::RationalNumber>;
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<ValueType> computeLongRunAverageStateDistribution(Environment const& env);
+                std::vector<ValueType> computeLongRunAverageStateDistribution(Environment const& env, uint64_t const& initialState);
+                std::vector<ValueType> 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<ValueType> computeBsccReachabilityProbabilities(Environment const& env, ValueGetter const& initialDistributionGetter);
+                
+                /*!
+                 * Computes the long run average (steady state) distribution for the given BSCC.
+                 */
+                std::vector<ValueType> computeSteadyStateDistrForBscc(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc);
+
+                
                 std::pair<bool, ValueType> 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<ValueType, std::vector<ValueType>> 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<ValueType, std::vector<ValueType>> computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc);
+                std::pair<ValueType, std::vector<ValueType>> computeLraForBsccSteadyStateDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc);
   
                 std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> buildSspMatrixVector(std::vector<ValueType> const& bsccLraValues, std::vector<uint64_t> 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 <typename SparseCtmcModelType>
+        std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<SparseCtmcModelType>::computeSteadyStateDistribution(Environment const& env) {
+            // Initialize helper
+            storm::modelchecker::helper::SparseDeterministicInfiniteHorizonHelper<ValueType> helper(this->getModel().getTransitionMatrix());
+            
+            // Compute result
+            std::vector<ValueType> 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<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(numInitStates);
+                result = helper.computeLongRunAverageStateDistribution(env, [&initialStates,&initProb](uint64_t const& stateIndex) { return initialStates.get(stateIndex) ? initProb : storm::utility::zero<ValueType>(); });
+            }
+
+            // Return CheckResult
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(result)));
+        }
+
+        
         template class SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>>;
 
 #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<CheckResult> computeLongRunAverageRewards(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::LongRunAverageRewardFormula, ValueType> const& checkTask) override;
             virtual std::unique_ptr<CheckResult> computeReachabilityTimes(Environment const& env, storm::logic::RewardMeasureType rewardMeasureType, CheckTask<storm::logic::EventuallyFormula, ValueType> const& checkTask) override;
             virtual std::unique_ptr<CheckResult> checkQuantileFormula(Environment const& env, CheckTask<storm::logic::QuantileFormula, ValueType> 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<CheckResult> 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;