diff --git a/src/storm/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/storm/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
index 010fe8b23..dd467f6d0 100644
--- a/src/storm/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
+++ b/src/storm/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
@@ -186,9 +186,9 @@ namespace storm {
             
             // Then, we convert the reduced matrix to a more flexible format to be able to perform state elimination more easily.
             storm::storage::FlexibleSparseMatrix<ValueType> flexibleMatrix(transitionMatrix);
-            flexibleMatrix.createSubmatrix(maybeStates, maybeStates);
+            flexibleMatrix.filterEntries(maybeStates, maybeStates);
             storm::storage::FlexibleSparseMatrix<ValueType> flexibleBackwardTransitions(backwardTransitions);
-            flexibleBackwardTransitions.createSubmatrix(maybeStates, maybeStates);
+            flexibleBackwardTransitions.filterEntries(maybeStates, maybeStates);
             auto conversionEnd = std::chrono::high_resolution_clock::now();
             
             std::chrono::high_resolution_clock::time_point modelCheckingStart = std::chrono::high_resolution_clock::now();
diff --git a/src/storm/modelchecker/region/ApproximationModel.cpp b/src/storm/modelchecker/region/ApproximationModel.cpp
index b79fee942..eabea80f1 100644
--- a/src/storm/modelchecker/region/ApproximationModel.cpp
+++ b/src/storm/modelchecker/region/ApproximationModel.cpp
@@ -33,7 +33,7 @@ namespace storm {
                     this->computeRewards=true;
                     STORM_LOG_THROW(this->typeOfParametricModel==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Approximation with rewards is only implemented for Dtmcs");
                     STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should be unique");
-                    STORM_LOG_THROW(parametricModel.getUniqueRewardModel().hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should have state rewards only");
+                    STORM_LOG_THROW(parametricModel.getUniqueRewardModel().hasStateActionRewards() && !parametricModel.getUniqueRewardModel().hasStateRewards() && !parametricModel.getUniqueRewardModel().hasTransitionRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should have state action rewards only");
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << formula << ". Approximation model only supports eventually or reachability reward formulae.");
                 }
@@ -205,13 +205,13 @@ namespace storm {
                 STORM_LOG_THROW(this->vectorData.vector.size()==this->matrixData.matrix.getRowCount(), storm::exceptions::UnexpectedException, "The size of the eq-sys vector does not match to the number of rows in the eq-sys matrix");
                 this->vectorData.assignment.reserve(vectorData.vector.size());
                 
-                // run through the state reward vector of the parametric model.
+                // run through the state action reward vector of the parametric model.
                 // Constant entries can be set directly.
                 // For Parametric entries we set a dummy value and insert the corresponding function and the assignment
                 ConstantType dummyNonZeroValue = storm::utility::one<ConstantType>();
                 auto vectorIt = this->vectorData.vector.begin();
                 for(auto oldState : this->maybeStates){
-                    if(storm::utility::isConstant(parametricModel.getUniqueRewardModel().getStateRewardVector()[oldState])){
+                    if(storm::utility::isConstant(parametricModel.getUniqueRewardModel().getStateActionRewardVector()[oldState])){
                         ConstantType reward = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel().getStateRewardVector()[oldState]));
                         //Add one of these entries for every row in the row group of oldState
                         for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[oldState+1]; ++matrixRow){
@@ -220,7 +220,7 @@ namespace storm {
                         }
                     } else {
                         std::set<VariableType> occurringRewVariables;
-                        storm::utility::region::gatherOccurringVariables(parametricModel.getUniqueRewardModel().getStateRewardVector()[oldState], occurringRewVariables);
+                        storm::utility::region::gatherOccurringVariables(parametricModel.getUniqueRewardModel().getStateActionRewardVector()[oldState], occurringRewVariables);
                         // For each row in the row group of oldState, we get the corresponding substitution and insert the FunctionSubstitution
                         for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[oldState+1]; ++matrixRow){
                             //Extend the substitution for the probabilities with the reward parameters
@@ -230,7 +230,7 @@ namespace storm {
                                 substitution.insert(typename std::map<VariableType, RegionBoundary>::value_type(rewardVar, RegionBoundary::UNSPECIFIED));
                             }
                             // insert the FunctionSubstitution
-                            auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel().getStateRewardVector()[oldState], this->matrixData.rowSubstitutions[matrixRow]), dummyNonZeroValue)).first;
+                            auto functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel().getStateActionRewardVector()[oldState], this->matrixData.rowSubstitutions[matrixRow]), dummyNonZeroValue)).first;
                             //insert assignment and dummy data
                             this->vectorData.assignment.emplace_back(vectorIt, functionsIt->second);
                             *vectorIt = dummyNonZeroValue;
diff --git a/src/storm/modelchecker/region/SamplingModel.cpp b/src/storm/modelchecker/region/SamplingModel.cpp
index 5c312d74d..ef4ff39b3 100644
--- a/src/storm/modelchecker/region/SamplingModel.cpp
+++ b/src/storm/modelchecker/region/SamplingModel.cpp
@@ -43,7 +43,7 @@ namespace storm {
                     this->computeRewards=true;
                     STORM_LOG_THROW(this->typeOfParametricModel==storm::models::ModelType::Dtmc, storm::exceptions::InvalidArgumentException, "Sampling with rewards is only implemented for Dtmcs");
                     STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should be unique");
-                    STORM_LOG_THROW(parametricModel.getUniqueRewardModel().hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should have state rewards only");
+                    STORM_LOG_THROW(parametricModel.getUniqueRewardModel().hasStateActionRewards() && !parametricModel.getUniqueRewardModel().hasStateRewards() && !parametricModel.getUniqueRewardModel().hasTransitionRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sempling model should have state action rewards only");
                     STORM_LOG_THROW(formula->getSubformula().isEventuallyFormula(), storm::exceptions::InvalidArgumentException, "The subformula should be a reachabilityreward formula");
                     STORM_LOG_THROW(formula->getSubformula().asEventuallyFormula().getSubformula().isInFragment(storm::logic::propositional()), storm::exceptions::InvalidArgumentException, "The subsubformula should be a propositional formula");
                 } else {
@@ -131,7 +131,7 @@ namespace storm {
                     std::vector<ConstantType> b;
                     if(this->computeRewards){
                         b.resize(submatrix.getRowCount());
-                        storm::utility::vector::selectVectorValues(b, this->maybeStates, instantiatedModel.getUniqueRewardModel().getStateRewardVector());
+                        storm::utility::vector::selectVectorValues(b, this->maybeStates, instantiatedModel.getUniqueRewardModel().getStateActionRewardVector());
                     } else {
                         b = instantiatedModel.getTransitionMatrix().getConstrainedRowSumVector(this->maybeStates, this->targetStates);
                     }
diff --git a/src/storm/models/sparse/StandardRewardModel.cpp b/src/storm/models/sparse/StandardRewardModel.cpp
index 4f43c62fb..78f4c39af 100644
--- a/src/storm/models/sparse/StandardRewardModel.cpp
+++ b/src/storm/models/sparse/StandardRewardModel.cpp
@@ -232,6 +232,34 @@ namespace storm {
                 return result;
             }
             
+            template<typename ValueType>
+            template<typename MatrixValueType>
+            storm::storage::BitVector StandardRewardModel<ValueType>::getStatesWithZeroReward(storm::storage::SparseMatrix<MatrixValueType> const& transitionMatrix) const {
+                storm::storage::BitVector result = this->hasStateRewards() ? storm::utility::vector::filterZero(this->getStateRewardVector()) : storm::storage::BitVector(transitionMatrix.getRowGroupCount(), true);
+                if (this->hasStateActionRewards()) {
+                    for (uint_fast64_t state = 0; state < transitionMatrix.getRowGroupCount(); ++state) {
+                        for (uint_fast64_t row = transitionMatrix.getRowGroupIndices()[state]; row < transitionMatrix.getRowGroupIndices()[state+1]; ++row) {
+                            if(!storm::utility::isZero(this->getStateActionRewardVector()[row])) {
+                                result.set(state, false);
+                                break;
+                            }
+                        }
+                    }
+                }
+                if (this->hasTransitionRewards()) {
+                    for (uint_fast64_t state = 0; state < transitionMatrix.getRowGroupCount(); ++state) {
+                        for (auto const& rewardMatrixEntry : this->getTransitionRewardMatrix().getRowGroup(state)) {
+                            if(!storm::utility::isZero(rewardMatrixEntry.getValue())) {
+                                result.set(state, false);
+                                break;
+                            }
+                        }
+                    }
+                }
+                return result;
+            }
+
+            
             template<typename ValueType>
             void StandardRewardModel<ValueType>::setStateActionRewardValue(uint_fast64_t row, ValueType const& value) {
                 this->optionalStateActionRewardVector.get()[row] = value;
@@ -304,6 +332,7 @@ namespace storm {
             template std::vector<double> StandardRewardModel<double>::getTotalRewardVector(uint_fast64_t numberOfRows, storm::storage::SparseMatrix<double> const& transitionMatrix, storm::storage::BitVector const& filter) const;
             template std::vector<double> StandardRewardModel<double>::getTotalRewardVector(storm::storage::SparseMatrix<double> const& transitionMatrix, std::vector<double> const& weights, bool scaleTransAndActions) const;
             template std::vector<double> StandardRewardModel<double>::getTotalActionRewardVector(storm::storage::SparseMatrix<double> const& transitionMatrix,  std::vector<double> const& stateRewardWeights) const;
+            template storm::storage::BitVector StandardRewardModel<double>::getStatesWithZeroReward(storm::storage::SparseMatrix<double> const& transitionMatrix) const;
             template void StandardRewardModel<double>::reduceToStateBasedRewards(storm::storage::SparseMatrix<double> const& transitionMatrix, bool reduceToStateRewards);
             template void StandardRewardModel<double>::setStateActionReward(uint_fast64_t choiceIndex, double const & newValue);
             template void StandardRewardModel<double>::setStateReward(uint_fast64_t state, double const & newValue);
@@ -325,6 +354,7 @@ namespace storm {
             template std::vector<storm::RationalNumber> StandardRewardModel<storm::RationalNumber>::getTotalRewardVector(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix) const;
             template std::vector<storm::RationalNumber> StandardRewardModel<storm::RationalNumber>::getTotalRewardVector(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, std::vector<storm::RationalNumber> const& weights, bool scaleTransAndActions) const;
             template std::vector<storm::RationalNumber> StandardRewardModel<storm::RationalNumber>::getTotalActionRewardVector(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix,  std::vector<storm::RationalNumber> const& stateRewardWeights) const;
+            template storm::storage::BitVector StandardRewardModel<storm::RationalNumber>::getStatesWithZeroReward(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix) const;
             template void StandardRewardModel<storm::RationalNumber>::reduceToStateBasedRewards(storm::storage::SparseMatrix<storm::RationalNumber> const& transitionMatrix, bool reduceToStateRewards);
             template void StandardRewardModel<storm::RationalNumber>::setStateActionReward(uint_fast64_t choiceIndex, storm::RationalNumber const & newValue);
             template void StandardRewardModel<storm::RationalNumber>::setStateReward(uint_fast64_t state, storm::RationalNumber const & newValue);
@@ -334,6 +364,8 @@ namespace storm {
             template std::vector<storm::RationalFunction> StandardRewardModel<storm::RationalFunction>::getTotalRewardVector(uint_fast64_t numberOfRows, storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, storm::storage::BitVector const& filter) const;
             template std::vector<storm::RationalFunction> StandardRewardModel<storm::RationalFunction>::getTotalRewardVector(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix) const;
             template std::vector<storm::RationalFunction> StandardRewardModel<storm::RationalFunction>::getTotalRewardVector(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, std::vector<storm::RationalFunction> const& weights, bool scaleTransAndActions) const;
+            template storm::storage::BitVector StandardRewardModel<storm::RationalFunction>::getStatesWithZeroReward(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix) const;
+
             template std::vector<storm::RationalFunction> StandardRewardModel<storm::RationalFunction>::getTotalActionRewardVector(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix,  std::vector<storm::RationalFunction> const& stateRewardWeights) const;
             template void StandardRewardModel<storm::RationalFunction>::reduceToStateBasedRewards(storm::storage::SparseMatrix<storm::RationalFunction> const& transitionMatrix, bool reduceToStateRewards);
             template void StandardRewardModel<storm::RationalFunction>::setStateActionReward(uint_fast64_t choiceIndex, storm::RationalFunction const & newValue);
diff --git a/src/storm/models/sparse/StandardRewardModel.h b/src/storm/models/sparse/StandardRewardModel.h
index 7251ac656..e98430d64 100644
--- a/src/storm/models/sparse/StandardRewardModel.h
+++ b/src/storm/models/sparse/StandardRewardModel.h
@@ -230,6 +230,15 @@ namespace storm {
                 template<typename MatrixValueType>
                 std::vector<ValueType> getTotalActionRewardVector(storm::storage::SparseMatrix<MatrixValueType> const& transitionMatrix,  std::vector<MatrixValueType> const& stateRewardWeights) const;
 
+                /*!
+                 * Returns the set of states at which a all rewards (state-, action- and transition-rewards) are zero.
+                 *
+                 * @param transitionMatrix the transition matrix of the model (used to determine the actions and transitions that belong to a state)
+                 * @ return a vector representing all states at which the reward is zero
+                 */
+                 template<typename MatrixValueType>
+                 storm::storage::BitVector getStatesWithZeroReward(storm::storage::SparseMatrix<MatrixValueType> const& transitionMatrix) const;
+                 
                 /*!
                  * Sets the given value in the state-action reward vector at the given row. This assumes that the reward
                  * model has state-action rewards.
diff --git a/src/storm/storage/FlexibleSparseMatrix.cpp b/src/storm/storage/FlexibleSparseMatrix.cpp
index d25450bd9..29308f410 100644
--- a/src/storm/storage/FlexibleSparseMatrix.cpp
+++ b/src/storm/storage/FlexibleSparseMatrix.cpp
@@ -109,16 +109,12 @@ namespace storm {
         
         template<typename ValueType>
         typename FlexibleSparseMatrix<ValueType>::index_type FlexibleSparseMatrix<ValueType>::getRowGroupCount() const {
-            return rowGroupIndices.size();
+            return rowGroupIndices.size() - 1;
         }
         
         template<typename ValueType>
         typename FlexibleSparseMatrix<ValueType>::index_type FlexibleSparseMatrix<ValueType>::getRowGroupSize(index_type group) const {
-            if (group == getRowGroupCount() - 1) {
-                return getRowCount() - rowGroupIndices[group];
-            } else {
-                return rowGroupIndices[group + 1] - rowGroupIndices[group];
-            }
+            return rowGroupIndices[group + 1] - rowGroupIndices[group];
         }
         
         template<typename ValueType>
@@ -157,9 +153,9 @@ namespace storm {
         bool FlexibleSparseMatrix<ValueType>::hasTrivialRowGrouping() const {
             return trivialRowGrouping;
         }
-
+        
         template<typename ValueType>
-        void FlexibleSparseMatrix<ValueType>::createSubmatrix(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) {
+        void FlexibleSparseMatrix<ValueType>::filterEntries(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) {
             for (uint_fast64_t rowIndex = 0; rowIndex < this->data.size(); ++rowIndex) {
                 auto& row = this->data[rowIndex];
                 if (!rowConstraint.get(rowIndex)) {
@@ -179,12 +175,83 @@ namespace storm {
         
         template<typename ValueType>
         storm::storage::SparseMatrix<ValueType> FlexibleSparseMatrix<ValueType>::createSparseMatrix() {
-            storm::storage::SparseMatrixBuilder<ValueType> matrixBuilder(getRowCount(), getColumnCount());
-            for (uint_fast64_t rowIndex = 0; rowIndex < getRowCount(); ++rowIndex) {
-                auto& row = this->data[rowIndex];
+            uint_fast64_t numEntries = 0;
+            for (auto const& row : this->data) {
+                numEntries += row.size();
+            }
+            
+            storm::storage::SparseMatrixBuilder<ValueType> matrixBuilder(getRowCount(), getColumnCount(), numEntries, hasTrivialRowGrouping(), hasTrivialRowGrouping() ? 0 : getRowGroupCount());
+            uint_fast64_t currRowIndex = 0;
+            auto rowGroupIndexIt = getRowGroupIndices().begin();
+            for (auto const& row : this->data) {
+                if(!hasTrivialRowGrouping()) {
+                    while (currRowIndex >= *rowGroupIndexIt) {
+                        matrixBuilder.newRowGroup(currRowIndex);
+                        ++rowGroupIndexIt;
+                    }
+                }
                 for (auto const& entry : row) {
-                    matrixBuilder.addNextValue(rowIndex, entry.getColumn(), entry.getValue());
+                    matrixBuilder.addNextValue(currRowIndex, entry.getColumn(), entry.getValue());
+                }
+                ++currRowIndex;
+            }
+            // The matrix might end with one or more empty row groups
+            if(!hasTrivialRowGrouping()) {
+                while (currRowIndex >= *rowGroupIndexIt) {
+                    matrixBuilder.newRowGroup(currRowIndex);
+                    ++rowGroupIndexIt;
+                }
+            }
+            return matrixBuilder.build();
+        }
+        
+        template<typename ValueType>
+        storm::storage::SparseMatrix<ValueType> FlexibleSparseMatrix<ValueType>::createSparseMatrix(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint) {
+            uint_fast64_t numEntries = 0;
+            for (auto const& rowIndex : rowConstraint) {
+                auto const& row = data[rowIndex];
+                for(auto const& entry : row) {
+                    if (columnConstraint.get(entry.getColumn())) {
+                        ++numEntries;
+                    }
+                }
+            }
+            uint_fast64_t numRowGroups = 0;
+            if (!hasTrivialRowGrouping()) {
+                auto lastRowGroupIndexIt = getRowGroupIndices().end() - 1;
+                auto rowGroupIndexIt = getRowGroupIndices().begin();
+                while (rowGroupIndexIt != lastRowGroupIndexIt) {
+                    // Check whether the rowGroup will be nonempty
+                    if(rowConstraint.getNextSetIndex(*rowGroupIndexIt) < *(++rowGroupIndexIt)) {
+                        ++numRowGroups;
+                    }
+                }
+            }
+            
+            std::vector<uint_fast64_t> oldToNewColumnIndexMapping(getColumnCount(), getColumnCount());
+            uint_fast64_t newColumnIndex = 0;
+            for (auto const& oldColumnIndex : columnConstraint) {
+                oldToNewColumnIndexMapping[oldColumnIndex] = newColumnIndex++;
+            }
+            
+            storm::storage::SparseMatrixBuilder<ValueType> matrixBuilder(rowConstraint.getNumberOfSetBits(), newColumnIndex, numEntries, hasTrivialRowGrouping(), numRowGroups);
+            uint_fast64_t currRowIndex = 0;
+            auto rowGroupIndexIt = getRowGroupIndices().begin();
+            for (auto const& oldRowIndex : rowConstraint) {
+                if(!hasTrivialRowGrouping() && oldRowIndex >= *rowGroupIndexIt) {
+                    matrixBuilder.newRowGroup(currRowIndex);
+                    // Skip empty row groups
+                    do {
+                        ++rowGroupIndexIt;
+                    } while (currRowIndex >= *rowGroupIndexIt);
+                }
+                auto const& row = data[oldRowIndex];
+                for (auto const& entry : row) {
+                    if(columnConstraint.get(entry.getColumn())) {
+                        matrixBuilder.addNextValue(currRowIndex, oldToNewColumnIndexMapping[entry.getColumn()], entry.getValue());
+                    }
                 }
+                ++currRowIndex;
             }
             return matrixBuilder.build();
         }
@@ -237,7 +304,7 @@ namespace storm {
                 FlexibleIndex rowGroupCount = matrix.getRowGroupCount();
                 for (FlexibleIndex rowGroup = 0; rowGroup < rowGroupCount; ++rowGroup) {
                     out << "\t---- group " << rowGroup << "/" << (rowGroupCount - 1) << " ---- " << std::endl;
-                    FlexibleIndex endRow = rowGroup+1 < rowGroupCount ? matrix.rowGroupIndices[rowGroup + 1] : matrix.getRowCount();
+                    FlexibleIndex endRow = matrix.rowGroupIndices[rowGroup + 1];
                     // Iterate over all rows.
                     for (FlexibleIndex row = matrix.rowGroupIndices[rowGroup]; row < endRow; ++row) {
                         // Print the actual row.
diff --git a/src/storm/storage/FlexibleSparseMatrix.h b/src/storm/storage/FlexibleSparseMatrix.h
index af2b0b7ec..92d385c81 100644
--- a/src/storm/storage/FlexibleSparseMatrix.h
+++ b/src/storm/storage/FlexibleSparseMatrix.h
@@ -158,19 +158,31 @@ namespace storm {
             bool hasTrivialRowGrouping() const;
 
             /*!
-             * Creates a submatrix of the current matrix in place by dropping all rows and columns whose bits are not
-             * set to one in the given bit vector.
+             * Erases all entries whose row and column does not satisfy the given rowConstraint and the given columnConstraint
              *
-             * @param rowConstraint A bit vector indicating which rows to keep.
-             * @param columnConstraint A bit vector indicating which columns to keep.
+             * @param rowConstraint A bit vector indicating which row entries to keep.
+             * @param columnConstraint A bit vector indicating which column entries to keep.
              */
-            void createSubmatrix(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint);
-
+            void filterEntries(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint);
+            
             /*!
              * Creates a sparse matrix from the flexible sparse matrix.
              * @return The sparse matrix.
              */
             storm::storage::SparseMatrix<ValueType> createSparseMatrix();
+                        
+            /*!
+             * Creates a sparse matrix from the flexible sparse matrix.
+             * Only the selected rows and columns will be considered.
+             * Empty rowGroups will be ignored
+             *
+             * @param rowConstraint A bit vector indicating which rows to keep.
+             * @param columnConstraint A bit vector indicating which columns to keep.
+
+             *
+             * @return The sparse matrix.
+             */
+            storm::storage::SparseMatrix<ValueType> createSparseMatrix(storm::storage::BitVector const& rowConstraint, storm::storage::BitVector const& columnConstraint);
 
             /*!
              * Checks whether the given state has a self-loop with an arbitrary probability in the probability matrix.
diff --git a/src/storm/transformer/GoalStateMerger.cpp b/src/storm/transformer/GoalStateMerger.cpp
new file mode 100644
index 000000000..69ac6cd16
--- /dev/null
+++ b/src/storm/transformer/GoalStateMerger.cpp
@@ -0,0 +1,190 @@
+#include "storm/transformer/GoalStateMerger.h"
+
+#include <limits>
+#include <memory>
+
+#include "storm/utility/constants.h"
+#include "storm/utility/macros.h"
+#include "storm/utility/vector.h"
+#include "storm/models/sparse/Dtmc.h"
+#include "storm/models/sparse/Mdp.h"
+#include "storm/models/sparse/StandardRewardModel.h"
+
+#include "storm/exceptions/InvalidArgumentException.h"
+
+
+namespace storm {
+    namespace transformer {
+        
+        template <typename SparseModelType>
+        GoalStateMerger<SparseModelType>::GoalStateMerger(SparseModelType const& model) : originalModel(model) {
+            // Intentionally left empty
+        }
+            
+        template <typename SparseModelType>
+        std::shared_ptr<SparseModelType> GoalStateMerger<SparseModelType>::mergeTargetAndSinkStates(storm::storage::BitVector const& maybeStates, storm::storage::BitVector& targetStates, storm::storage::BitVector& sinkStates, std::vector<std::string> const& selectedRewardModels) {
+            STORM_LOG_THROW(maybeStates.isDisjointFrom(targetStates) && targetStates.isDisjointFrom(sinkStates) && sinkStates.isDisjointFrom(maybeStates), storm::exceptions::InvalidArgumentException, "maybestates, targetstates, and sinkstates are assumed to be disjoint when creating the submodel. However, this is not the case.");
+  
+            boost::optional<uint_fast64_t> targetState, sinkState;
+            auto builder = initializeTransitionMatrixBuilder(maybeStates, targetStates, sinkStates, targetState, sinkState);
+            auto transitionMatrix = buildTransitionMatrix(maybeStates, targetStates, sinkStates, targetState, sinkState, builder);
+            
+            uint_fast64_t resNumStates = transitionMatrix.getRowGroupCount();
+                
+            // Get the labeling for the initial states
+            storm::storage::BitVector initialStates = originalModel.getInitialStates() % maybeStates;
+            initialStates.resize(resNumStates, false);
+            if(!originalModel.getInitialStates().isDisjointFrom(targetStates)) {
+                initialStates.set(*targetState, true);
+            }
+            if(!originalModel.getInitialStates().isDisjointFrom(sinkStates)) {
+                initialStates.set(*sinkState, true);
+            }
+            storm::models::sparse::StateLabeling labeling(resNumStates);
+            labeling.addLabel("init", std::move(initialStates));
+                
+            // Get the reward models
+            std::unordered_map<std::string, typename SparseModelType::RewardModelType> rewardModels;
+            for (auto rewardModelName : selectedRewardModels) {
+                auto origTotalRewards = originalModel.getRewardModel(rewardModelName).getTotalRewardVector(originalModel.getTransitionMatrix());
+                auto resTotalRewards = storm::utility::vector::filterVector(origTotalRewards, maybeStates);
+                resTotalRewards.resize(resNumStates, storm::utility::zero<typename SparseModelType::RewardModelType::ValueType>());
+                rewardModels.insert(std::make_pair(rewardModelName, typename SparseModelType::RewardModelType(boost::none, resTotalRewards)));
+            }
+                
+            // modify the given target and sink states
+            targetStates = storm::storage::BitVector(resNumStates, false);
+            if(targetState) {
+                targetStates.set(*targetState, true);
+            }
+            sinkStates = storm::storage::BitVector(resNumStates, false);
+            if(sinkState) {
+                sinkStates.set(*sinkState, true);
+            }
+                
+            // Return the result
+            return std::make_shared<SparseModelType>(std::move(transitionMatrix), std::move(labeling), std::move(rewardModels));
+        }
+        
+        template <typename SparseModelType>
+        storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType> GoalStateMerger<SparseModelType>::initializeTransitionMatrixBuilder(storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& sinkStates, boost::optional<uint_fast64_t>& newTargetState, boost::optional<uint_fast64_t>& newSinkState) {
+            
+            storm::storage::SparseMatrix<typename SparseModelType::ValueType> const& origMatrix = originalModel.getTransitionMatrix();
+                
+            // Get the number of rows, cols and entries that the resulting transition matrix will have.
+            uint_fast64_t resNumStates(maybeStates.getNumberOfSetBits()), resNumActions(0), resNumTransitions(0);
+            bool targetStateRequired = !originalModel.getInitialStates().isDisjointFrom(targetStates);
+            bool sinkStateRequired = !originalModel.getInitialStates().isDisjointFrom(targetStates);
+            for( auto state : maybeStates) {
+                resNumActions += origMatrix.getRowGroupSize(state);
+                bool hasTransitionToTarget(false), hasTransitionToSink(false);
+                auto const& endOfRowGroup = origMatrix.getRowGroupIndices()[state+1];
+                for (uint_fast64_t row = origMatrix.getRowGroupIndices()[state]; row < endOfRowGroup; ++row) {
+                    for (auto const& entry : origMatrix.getRow(row)) {
+                        if(maybeStates.get(entry.getColumn())) {
+                            ++resNumTransitions;
+                        } else if (targetStates.get(entry.getColumn())) {
+                            hasTransitionToTarget = true;
+                        } else if (sinkStates.get(entry.getColumn())) {
+                            hasTransitionToSink = true;
+                        } else {
+                            STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "There is a transition originating from a maybestate that does not lead to a maybe-, target-, or sinkstate.");
+                        }
+                    }
+                    if(hasTransitionToTarget) {
+                        ++resNumTransitions;
+                        targetStateRequired = true;
+                    }
+                    if(hasTransitionToSink) {
+                        ++resNumTransitions;
+                        sinkStateRequired = true;
+                    }
+                }
+            }
+            
+            // Get the index of the target/ sink state in the resulting model (if these states will exist)
+            if(targetStateRequired) {
+                newTargetState = resNumStates;
+                ++resNumStates;
+                ++resNumActions;
+                ++resNumTransitions;
+            }
+            if(sinkStateRequired) {
+                newSinkState = resNumStates;
+                ++resNumStates;
+                ++resNumActions;
+                ++resNumTransitions;
+            }
+            
+            return storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType>(resNumActions, resNumStates, resNumTransitions, true, true, resNumStates);
+
+        }
+        
+        
+        
+        template <typename SparseModelType>
+        storm::storage::SparseMatrix<typename SparseModelType::ValueType> GoalStateMerger<SparseModelType>::buildTransitionMatrix(storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& sinkStates, boost::optional<uint_fast64_t> const& newTargetState, boost::optional<uint_fast64_t>const& newSinkState, storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType>& builder) {
+        
+            // Get a Mapping that yields for each column in the old matrix the corresponding column in the new matrix
+            std::vector<uint_fast64_t> oldToNewIndexMap(maybeStates.size(), std::numeric_limits<uint_fast64_t>::max()); // init with some invalid state
+            uint_fast64_t newStateIndex = 0;
+            for (auto maybeState : maybeStates) {
+                oldToNewIndexMap[maybeState] = newStateIndex;
+                ++newStateIndex;
+            }
+                
+            // Build the transition matrix
+            storm::storage::SparseMatrix<typename SparseModelType::ValueType> const& origMatrix = originalModel.getTransitionMatrix();
+            uint_fast64_t currRow = 0;
+            for (auto state : maybeStates) {
+                builder.newRowGroup(currRow);
+                boost::optional<typename SparseModelType::ValueType> targetProbability, sinkProbability;
+                auto const& endOfRowGroup = origMatrix.getRowGroupIndices()[state+1];
+                for (uint_fast64_t row = origMatrix.getRowGroupIndices()[state]; row < endOfRowGroup; ++row) {
+                    for (auto const& entry : origMatrix.getRow(row)) {
+                        if(maybeStates.get(entry.getColumn())) {
+                            builder.addNextValue(currRow, oldToNewIndexMap[entry.getColumn()], entry.getValue());
+                        } else if (targetStates.get(entry.getColumn())) {
+                            targetProbability = targetProbability.is_initialized() ? *targetProbability + entry.getValue() : entry.getValue();
+                        } else if (sinkStates.get(entry.getColumn())) {
+                            sinkProbability = sinkProbability.is_initialized() ? *sinkProbability + entry.getValue() : entry.getValue();
+                        } else {
+                            STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "There is a transition originating from a maybestate that does not lead to a maybe-, target-, or sinkstate.");
+                        }
+                    }
+                    if(targetProbability) {
+                        assert(newTargetState);
+                        builder.addNextValue(currRow, *newTargetState, storm::utility::simplify(*targetProbability));
+                    }
+                    if(sinkProbability) {
+                        assert(newSinkState);
+                        builder.addNextValue(currRow, *newSinkState, storm::utility::simplify(*sinkProbability));
+                    }
+                    ++currRow;
+                }
+            }
+            // Add the selfloops at target and sink
+            if(newTargetState) {
+                builder.newRowGroup(currRow);
+                builder.addNextValue(currRow, *newTargetState, storm::utility::one<typename SparseModelType::ValueType>());
+                ++currRow;
+            }
+            if(newSinkState) {
+                builder.newRowGroup(currRow);
+                builder.addNextValue(currRow, *newSinkState, storm::utility::one<typename SparseModelType::ValueType>());
+                ++currRow;
+            }
+            
+            return builder.build();
+        }
+        
+        template class GoalStateMerger<storm::models::sparse::Dtmc<double>>;
+        template class GoalStateMerger<storm::models::sparse::Mdp<double>>;
+        template class GoalStateMerger<storm::models::sparse::Dtmc<storm::RationalNumber>>;
+        template class GoalStateMerger<storm::models::sparse::Mdp<storm::RationalNumber>>;
+        template class GoalStateMerger<storm::models::sparse::Dtmc<storm::RationalFunction>>;
+        template class GoalStateMerger<storm::models::sparse::Mdp<storm::RationalFunction>>;
+
+
+    }
+}
diff --git a/src/storm/transformer/GoalStateMerger.h b/src/storm/transformer/GoalStateMerger.h
index 08d8cb684..712eb9af6 100644
--- a/src/storm/transformer/GoalStateMerger.h
+++ b/src/storm/transformer/GoalStateMerger.h
@@ -1,15 +1,12 @@
 #pragma once
 
-#include <limits>
 #include <memory>
+#include <string>
+#include <vector>
 #include <boost/optional.hpp>
 
-#include "storm/utility/constants.h"
-#include "storm/utility/macros.h"
-#include "storm/models/sparse/Dtmc.h"
-#include "storm/models/sparse/Mdp.h"
-
-#include "storm/exceptions/InvalidArgumentException.h"
+#include "storm/storage/BitVector.h"
+#include "storm/storage/SparseMatrix.h"
 
 
 namespace storm {
@@ -21,7 +18,9 @@ namespace storm {
         template <typename SparseModelType>
         class GoalStateMerger {
         public:
-                                    
+                      
+            GoalStateMerger(SparseModelType const& model);
+            
             /* Computes a submodel of the specified model that only considers the states given by maybeStates as well as
              *  * one target state to which all transitions to a state selected by targetStates are redirected and
              *  * one sink state to which all transitions to a state selected by sinkStates are redirected.
@@ -30,133 +29,31 @@ namespace storm {
              *  The two given bitvectors targetStates and sinkStates are modified such that they represent the corresponding state in the obtained submodel.
              *
              *  Notes:
-             *  * The resulting model will not have any rewardmodels, labels (other then "init") etc.
+             *  * The resulting model will not have any labels (other then "init") and only the selected reward models.
+             *  * Rewards are reduced to stateActionRewards.
+             *  * The target and sink states will not get any reward
              *  * It is assumed that the given set of maybeStates can only be left via either a target or a sink state. Otherwise an exception is thrown.
              *  * It is assumed that maybeStates, targetStates, and sinkStates are all disjoint. Otherwise an exception is thrown.
              */
-            static std::shared_ptr<SparseModelType> mergeTargetAndSinkStates(SparseModelType const& model, storm::storage::BitVector const& maybeStates, storm::storage::BitVector& targetStates, storm::storage::BitVector& sinkStates) {
-                STORM_LOG_THROW(maybeStates.isDisjointFrom(targetStates) && targetStates.isDisjointFrom(sinkStates) && sinkStates.isDisjointFrom(maybeStates), storm::exceptions::InvalidArgumentException, "maybestates, targetstates, and sinkstates are assumed to be disjoint when creating the submodel. However, this is not the case.");
-                storm::storage::SparseMatrix<typename SparseModelType::ValueType> const& origMatrix = model.getTransitionMatrix();
-                
-                // Get the number of rows, cols and entries that the resulting transition matrix will have.
-                uint_fast64_t resNumStates(maybeStates.getNumberOfSetBits()), resNumActions(0), resNumTransitions(0);
-                bool targetStateRequired = !model.getInitialStates().isDisjointFrom(targetStates);
-                bool sinkStateRequired = !model.getInitialStates().isDisjointFrom(targetStates);
-                for( auto state : maybeStates) {
-                    resNumActions += origMatrix.getRowGroupSize(state);
-                    bool hasTransitionToTarget(false), hasTransitionToSink(false);
-                    auto const& endOfRowGroup = origMatrix.getRowGroupIndices()[state+1];
-                    for (uint_fast64_t row = origMatrix.getRowGroupIndices()[state]; row < endOfRowGroup; ++row) {
-                        for (auto const& entry : origMatrix.getRow(row)) {
-                            if(maybeStates.get(entry.getColumn())) {
-                                ++resNumTransitions;
-                            } else if (targetStates.get(entry.getColumn())) {
-                                hasTransitionToTarget = true;
-                            } else if (sinkStates.get(entry.getColumn())) {
-                                hasTransitionToSink = true;
-                            } else {
-                                STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "There is a transition originating from a maybestate that does not lead to a maybe-, target-, or sinkstate.");
-                            }
-                        }
-                        if(hasTransitionToTarget) {
-                            ++resNumTransitions;
-                            targetStateRequired = true;
-                        }
-                        if(hasTransitionToSink) {
-                            ++resNumTransitions;
-                            sinkStateRequired = true;
-                        }
-                    }
-                }
-                
-                // Get the index of the target/ sink state in the resulting model (if these states will exist)
-                uint_fast64_t targetState(std::numeric_limits<uint_fast64_t>::max()), sinkState(std::numeric_limits<uint_fast64_t>::max()); // init with some invalid state
-                if(targetStateRequired) {
-                    targetState = resNumStates;
-                    ++resNumStates;
-                    ++resNumActions;
-                    ++resNumTransitions;
-                }
-                if(sinkStateRequired) {
-                    sinkState = resNumStates;
-                    ++resNumStates;
-                    ++resNumActions;
-                    ++resNumTransitions;
-                }
-                
-                // Get a Mapping that yields for each column in the old matrix the corresponding column in the new matrix
-                std::vector<uint_fast64_t> oldToNewIndexMap(maybeStates.size(), std::numeric_limits<uint_fast64_t>::max()); // init with some invalid state
-                uint_fast64_t newStateIndex = 0;
-                for (auto maybeState : maybeStates) {
-                    oldToNewIndexMap[maybeState] = newStateIndex;
-                    ++newStateIndex;
-                }
-                
-                // Build the transition matrix
-                storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType> builder(resNumActions, resNumStates, resNumTransitions, true, true, resNumStates);
-                uint_fast64_t currRow = 0;
-                for (auto state : maybeStates) {
-                    builder.newRowGroup(currRow);
-                    boost::optional<typename SparseModelType::ValueType> targetProbability, sinkProbability;
-                    auto const& endOfRowGroup = origMatrix.getRowGroupIndices()[state+1];
-                    for (uint_fast64_t row = origMatrix.getRowGroupIndices()[state]; row < endOfRowGroup; ++row) {
-                        for (auto const& entry : origMatrix.getRow(row)) {
-                            if(maybeStates.get(entry.getColumn())) {
-                                builder.addNextValue(currRow, oldToNewIndexMap[entry.getColumn()], entry.getValue());
-                            } else if (targetStates.get(entry.getColumn())) {
-                                targetProbability = targetProbability.is_initialized() ? *targetProbability + entry.getValue() : entry.getValue();
-                            } else if (sinkStates.get(entry.getColumn())) {
-                                sinkProbability = sinkProbability.is_initialized() ? *sinkProbability + entry.getValue() : entry.getValue();
-                            } else {
-                                STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "There is a transition originating from a maybestate that does not lead to a maybe-, target-, or sinkstate.");
-                            }
-                        }
-                        if(targetProbability) {
-                            builder.addNextValue(currRow, targetState, storm::utility::simplify(*targetProbability));
-                        }
-                        if(sinkProbability) {
-                            builder.addNextValue(currRow, sinkState, storm::utility::simplify(*sinkProbability));
-                        }
-                        ++currRow;
-                    }
-                }
-                // Add the selfloops at target and sink
-                if(targetStateRequired) {
-                    builder.newRowGroup(currRow);
-                    builder.addNextValue(currRow, targetState, storm::utility::one<typename SparseModelType::ValueType>());
-                    ++currRow;
-                }
-                if(sinkStateRequired) {
-                    builder.newRowGroup(currRow);
-                    builder.addNextValue(currRow, sinkState, storm::utility::one<typename SparseModelType::ValueType>());
-                    ++currRow;
-                }
-                
-                // Get the labeling for the initial states
-                storm::storage::BitVector initialStates = model.getInitialStates() % maybeStates;
-                initialStates.resize(resNumStates, false);
-                if(!model.getInitialStates().isDisjointFrom(targetStates)) {
-                    initialStates.set(targetState, true);
-                }
-                if(!model.getInitialStates().isDisjointFrom(sinkStates)) {
-                    initialStates.set(sinkState, true);
-                }
-                storm::models::sparse::StateLabeling labeling(resNumStates);
-                labeling.addLabel("init", std::move(initialStates));
-                
-                // modify the given target and sink states
-                targetStates = storm::storage::BitVector(resNumStates, false);
-                if(targetStateRequired) {
-                    targetStates.set(targetState, true);
-                }
-                sinkStates = storm::storage::BitVector(resNumStates, false);
-                if(sinkStateRequired) {
-                    sinkStates.set(sinkState, true);
-                }
-                
-                // Return the result
-                return std::make_shared<SparseModelType>(builder.build(), std::move(labeling));
-            }
+            std::shared_ptr<SparseModelType> mergeTargetAndSinkStates(storm::storage::BitVector const& maybeStates, storm::storage::BitVector& targetStates, storm::storage::BitVector& sinkStates, std::vector<std::string> const& selectedRewardModels = std::vector<std::string>());
+            
+        private:
+            SparseModelType const& originalModel;
+            
+            /*!
+             * Initializes the matrix builder for the transition matrix of the resulting model
+             *
+             * @param newTargetState Will be set to the index of the target state of the resulting model (Only if it has a target state)
+             * @param newSinkState Will be set to the index of the sink state of the resulting model (Only if it has a sink state)
+             */
+            storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType> initializeTransitionMatrixBuilder(storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& sinkStates, boost::optional<uint_fast64_t>& newTargetState, boost::optional<uint_fast64_t>& newSinkState);
+            
+            /*!
+             * Builds the transition matrix of the resulting model
+             */
+            storm::storage::SparseMatrix<typename SparseModelType::ValueType> buildTransitionMatrix(storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& sinkStates, boost::optional<uint_fast64_t> const& newTargetState, boost::optional<uint_fast64_t> const& newSinkState, storm::storage::SparseMatrixBuilder<typename SparseModelType::ValueType>& builder);
+
+            
         };
     }
 }
diff --git a/src/storm/transformer/SparseParametricDtmcSimplifier.cpp b/src/storm/transformer/SparseParametricDtmcSimplifier.cpp
index 72ef29b61..9ea8aecc3 100644
--- a/src/storm/transformer/SparseParametricDtmcSimplifier.cpp
+++ b/src/storm/transformer/SparseParametricDtmcSimplifier.cpp
@@ -9,6 +9,8 @@
 #include "storm/transformer/GoalStateMerger.h"
 #include "storm/utility/graph.h"
 
+#include <storm/exceptions/NotSupportedException.h>
+#include <storm/exceptions/UnexpectedException.h>
 
 namespace storm {
     namespace transformer {
@@ -30,11 +32,12 @@ namespace storm {
             storm::storage::BitVector psiStates = std::move(propositionalChecker.check(formula.getSubformula().asUntilFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
             std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(this->originalModel.getBackwardTransitions(), phiStates, psiStates);
             // Only consider the maybestates that are reachable from one initial state without hopping over a target (i.e., prob1) state
-            storm::storage::BitVector reachableGreater0States = storm::utility::graph::getReachableStates(this->originalModel.getTransitionMatrix(), this->originalModel.getInitialStates(), ~statesWithProbability01.first, statesWithProbability01.second);
+            storm::storage::BitVector reachableGreater0States = storm::utility::graph::getReachableStates(this->originalModel.getTransitionMatrix(), this->originalModel.getInitialStates() & ~statesWithProbability01.first, ~statesWithProbability01.first, statesWithProbability01.second);
             storm::storage::BitVector maybeStates = reachableGreater0States & ~statesWithProbability01.second;
             
             // obtain the resulting subsystem
-            this->simplifiedModel = storm::transformer::GoalStateMerger<SparseModelType>::mergeTargetAndSinkStates(this->originalModel, maybeStates, statesWithProbability01.second, statesWithProbability01.first);
+            storm::transformer::GoalStateMerger<SparseModelType> goalStateMerger(this->originalModel);
+            this->simplifiedModel = goalStateMerger.mergeTargetAndSinkStates(maybeStates, statesWithProbability01.second, statesWithProbability01.first);
             this->simplifiedModel->getStateLabeling().addLabel("target", statesWithProbability01.second);
             this->simplifiedModel->getStateLabeling().addLabel("sink", statesWithProbability01.first);
             
@@ -51,14 +54,76 @@ namespace storm {
         
         template<typename SparseModelType>
         bool SparseParametricDtmcSimplifier<SparseModelType>::simplifyForBoundedUntilProbabilities(storm::logic::ProbabilityOperatorFormula const& formula) {
-            // If this method was not overriden by any subclass, simplification is not possible
-            return false;
+            STORM_LOG_THROW(!formula.getSubformula().asBoundedUntilFormula().hasLowerBound(), storm::exceptions::NotSupportedException, "Lower step bounds are not supported.");
+            STORM_LOG_THROW(formula.getSubformula().asBoundedUntilFormula().hasUpperBound(), storm::exceptions::UnexpectedException, "Expected a bounded until formula with an upper bound.");
+            STORM_LOG_THROW(formula.getSubformula().asBoundedUntilFormula().hasIntegerUpperBound(), storm::exceptions::UnexpectedException, "Expected a bounded until formula with an integral upper bound.");
+            
+            // Get the prob0, target, and the maybeStates
+            storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> propositionalChecker(this->originalModel);
+            if(!propositionalChecker.canHandle(formula.getSubformula().asBoundedUntilFormula().getLeftSubformula()) || !propositionalChecker.canHandle(formula.getSubformula().asBoundedUntilFormula().getRightSubformula())) {
+                STORM_LOG_DEBUG("Can not simplify when Until-formula has non-propositional subformula(s). Formula: " << formula);
+                return false;
+            }
+            storm::storage::BitVector phiStates = std::move(propositionalChecker.check(formula.getSubformula().asBoundedUntilFormula().getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+            storm::storage::BitVector psiStates = std::move(propositionalChecker.check(formula.getSubformula().asBoundedUntilFormula().getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+            storm::storage::BitVector probGreater0States = storm::utility::graph::performProbGreater0(this->originalModel.getBackwardTransitions(), phiStates, psiStates, true, formula.getSubformula().asBoundedUntilFormula().getUpperBound().evaluateAsInt());
+            storm::storage::BitVector prob0States = ~probGreater0States;
+            
+            
+            // Only consider the maybestates that are reachable from one initial probGreater0 state within the given amount of steps and without hopping over a target state
+            storm::storage::BitVector reachableGreater0States = storm::utility::graph::getReachableStates(this->originalModel.getTransitionMatrix(), this->originalModel.getInitialStates() & probGreater0States, probGreater0States, psiStates, true, formula.getSubformula().asBoundedUntilFormula().getUpperBound().evaluateAsInt());
+            storm::storage::BitVector maybeStates = reachableGreater0States & ~psiStates;
+            
+            // obtain the resulting subsystem
+            storm::transformer::GoalStateMerger<SparseModelType> goalStateMerger(this->originalModel);
+            this->simplifiedModel = goalStateMerger.mergeTargetAndSinkStates(maybeStates, prob0States, psiStates);
+            this->simplifiedModel->getStateLabeling().addLabel("target", psiStates);
+            this->simplifiedModel->getStateLabeling().addLabel("sink", prob0States);
+            
+            // obtain the simplified formula for the simplified model
+            auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula const> ("target");
+            auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula const>(storm::logic::Formula::getTrueFormula(), labelFormula, boost::none, storm::logic::);
+            this->simplifiedFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula const>(eventuallyFormula, formula.getOperatorInformation());
+            
+            return true;
         }
         
         template<typename SparseModelType>
         bool SparseParametricDtmcSimplifier<SparseModelType>::simplifyForReachabilityRewards(storm::logic::RewardOperatorFormula const& formula) {
-            // If this method was not overriden by any subclass, simplification is not possible
-            return false;
+            typename SparseModelType::RewardModelType const& originalRewardModel = formula.hasRewardModelName() ? this->originalModel.getRewardModel(formula.getRewardModelName()) : this->originalModel.getUniqueRewardModel();
+            
+            // Get the prob1 and the maybeStates
+            storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> propositionalChecker(this->originalModel);
+            if(!propositionalChecker.canHandle(formula.getSubformula().asEventuallyFormula().getSubformula())) {
+                STORM_LOG_DEBUG("Can not simplify when reachability reward formula has non-propositional subformula(s). Formula: " << formula);
+                return false;
+            }
+            storm::storage::BitVector targetStates = std::move(propositionalChecker.check(formula.getSubformula().asEventuallyFormula().getSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector());
+            // The set of target states can be extended by the states that reach target with probability 1 without collecting any reward
+            targetStates = storm::utility::graph::performProb1(this->originalModel.getBackwardTransitions(), originalRewardModel.getStatesWithZeroReward(this->originalModel.getTransitionMatrix()), targetStates);
+            storm::storage::BitVector statesWithProb1 = storm::utility::graph::performProb1(this->originalModel.getBackwardTransitions(), storm::storage::BitVector(this->originalModel.getNumberOfStates(), true), targetStates);
+            storm::storage::BitVector infinityStates = ~statesWithProb1;
+            // Only consider the states that are reachable from an initial state without hopping over a target state
+            storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->originalModel.getTransitionMatrix(), this->originalModel.getInitialStates() & statesWithProb1, statesWithProb1, targetStates);
+            storm::storage::BitVector maybeStates = reachableStates & ~targetStates;
+            targetStates = targetStates & reachableStates;
+            
+            // obtain the resulting subsystem
+            std::vector<std::string> rewardModelNameAsVector(1, formula.hasRewardModelName() ? formula.getRewardModelName() : this->originalModel.getRewardModels().begin()->first);
+            storm::transformer::GoalStateMerger<SparseModelType> goalStateMerger(this->originalModel);
+            this->simplifiedModel = goalStateMerger.mergeTargetAndSinkStates(maybeStates, targetStates, infinityStates, rewardModelNameAsVector);
+            this->simplifiedModel->getStateLabeling().addLabel("target", targetStates);
+            this->simplifiedModel->getStateLabeling().addLabel("sink", infinityStates);
+            
+            // obtain the simplified formula for the simplified model
+            auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula const> ("target");
+            auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula const>(labelFormula, storm::logic::FormulaContext::Reward);
+            this->simplifiedFormula = std::make_shared<storm::logic::RewardOperatorFormula const>(eventuallyFormula, rewardModelNameAsVector.front(), formula.getOperatorInformation(), storm::logic::RewardMeasureType::Expectation);
+            
+            // Eliminate all states for which all outgoing transitions are constant
+            this->simplifiedModel = this->eliminateConstantDeterministicStates(*this->simplifiedModel);
+                
+            return true;
         }
         
         template<typename SparseModelType>
diff --git a/src/storm/transformer/SparseParametricModelSimplifier.cpp b/src/storm/transformer/SparseParametricModelSimplifier.cpp
index b972e20cd..58883468c 100644
--- a/src/storm/transformer/SparseParametricModelSimplifier.cpp
+++ b/src/storm/transformer/SparseParametricModelSimplifier.cpp
@@ -10,8 +10,10 @@
 #include "storm/solver/stateelimination/PrioritizedStateEliminator.h"
 #include "storm/storage/FlexibleSparseMatrix.h"
 #include "storm/utility/vector.h"
+
 #include "storm/exceptions/InvalidStateException.h"
 #include "storm/exceptions/NotImplementedException.h"
+#include "storm/exceptions/InvalidPropertyException.h"
 
 namespace storm {
     namespace transformer {
@@ -37,6 +39,7 @@ namespace storm {
                 }
             } else if (formula.isRewardOperatorFormula()) {
                 storm::logic::RewardOperatorFormula rewOpForm = formula.asRewardOperatorFormula();
+                STORM_LOG_THROW((rewOpForm.hasRewardModelName() && originalModel.hasRewardModel(rewOpForm.getRewardModelName())) || (!rewOpForm.hasRewardModelName() && originalModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model specified by formula " << formula << " is not available in the given model.");
                 if (rewOpForm.getSubformula().isReachabilityRewardFormula()) {
                     return simplifyForReachabilityRewards(rewOpForm);
                 } else if (rewOpForm.getSubformula().isCumulativeRewardFormula()) {
@@ -141,7 +144,6 @@ namespace storm {
             storm::solver::stateelimination::PrioritizedStateEliminator<typename SparseModelType::ValueType> stateEliminator(flexibleMatrix, flexibleBackwardTransitions, statesToEliminate, stateValues);
             stateEliminator.eliminateAll();
             
-            flexibleMatrix.createSubmatrix(keptRows, keptStates);
             stateValues = storm::utility::vector::filterVector(stateValues, keptRows);
             
             std::unordered_map<std::string, typename SparseModelType::RewardModelType> rewardModels;
@@ -149,7 +151,7 @@ namespace storm {
                 rewardModels.insert(std::make_pair(model.getRewardModels().begin()->first, typename SparseModelType::RewardModelType(boost::none, stateValues)));
             }
             
-            return std::make_shared<SparseModelType>(flexibleMatrix.createSparseMatrix(), model.getStateLabeling().getSubLabeling(keptStates), std::move(rewardModels));
+            return std::make_shared<SparseModelType>(flexibleMatrix.createSparseMatrix(keptRows, keptStates), model.getStateLabeling().getSubLabeling(keptStates), std::move(rewardModels));
         }