diff --git a/src/storm/logic/TimeBoundType.h b/src/storm/logic/TimeBoundType.h
index 7200cc89d..0bd6aba0c 100644
--- a/src/storm/logic/TimeBoundType.h
+++ b/src/storm/logic/TimeBoundType.h
@@ -24,7 +24,10 @@ namespace storm {
                 assert(rewardName != ""); // Empty reward name is reserved.
             }
 
-
+            TimeBoundType getType() const {
+                return type;
+            }
+            
             bool isStepBound() const {
                 return type ==  TimeBoundType::Steps;
             }
diff --git a/src/storm/modelchecker/multiobjective/Objective.h b/src/storm/modelchecker/multiobjective/Objective.h
index d5d269775..6b1f44f1f 100644
--- a/src/storm/modelchecker/multiobjective/Objective.h
+++ b/src/storm/modelchecker/multiobjective/Objective.h
@@ -13,57 +13,27 @@ namespace storm {
         namespace multiobjective {
             template <typename ValueType>
             struct Objective {
+                
                 // the original input formula
                 std::shared_ptr<storm::logic::Formula const> originalFormula;
                 
-                // the name of the considered reward model in the preprocessedModel
-                boost::optional<std::string> rewardModelName;
+                // The preprocessed (simplified) formula
+                std::shared_ptr<storm::logic::OperatorFormula const> formula;
                 
                 // True iff the complementary event is considered.
                 // E.g. if we consider P<1-t [F !"safe"] instead of P>=t [ G "safe"]
                 bool considersComplementaryEvent;
                 
-                // The probability/reward threshold for the preprocessed model (if originalFormula specifies one).
-                boost::optional<storm::logic::Bound> bound;
-                // The optimization direction for the preprocessed model
-                // if originalFormula does ot specifies one, the direction is derived from the bound.
-                storm::solver::OptimizationDirection optimizationDirection;
-                
-                // Lower and upper time/step/reward bouds
-                boost::optional<storm::logic::TimeBound> lowerTimeBound, upperTimeBound;
-                boost::optional<storm::logic::TimeBoundReference> timeBoundReference;
-                
+                // Limitations for the quantitative objective value (e.g. 0 <= value <= 1 for probabilities).
+                // Can be used to guide the underlying solver
                 boost::optional<ValueType> lowerResultBound, upperResultBound;
                 
                 void printToStream(std::ostream& out) const {
-                    out  << originalFormula->toString();
-                    out << " \t";
-                    out << "direction: ";
-                    out << optimizationDirection;
-                    out << " \t";
-                    out << "intern bound: ";
-                    if (bound){
-                        out << *bound;
-                    } else {
-                        out << " -none-    ";
-                    }
+                    out  << "Original: " << *originalFormula;
                     out << " \t";
-                    out << "time bounds: ";
-                    if (lowerTimeBound && upperTimeBound) {
-                        out << (lowerTimeBound->isStrict() ? "(" : "[") << lowerTimeBound->getBound() << "," << upperTimeBound->getBound() << (upperTimeBound->isStrict() ? ")" : "]");
-                    } else if (lowerTimeBound) {
-                        out << (lowerTimeBound->isStrict() ? ">" : ">=") << lowerTimeBound->getBound();
-                    } else if (upperTimeBound) {
-                        out << (upperTimeBound->isStrict() ? "<" : "<=") << upperTimeBound->getBound();
-                    } else {
-                        out << " -none-    ";
-                    }
-                    out << " \t";
-                    out << "intern reward model: ";
-                    if (rewardModelName) {
-                        out << *rewardModelName;
-                    } else {
-                        out << " -none-    ";
+                    out << "Preprocessed: " << *formula;
+                    if (considersComplementaryEvent) {
+                        out << " (Complementary event)";
                     }
                     out << " \t";
                     out << "result bounds: ";
diff --git a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp
index 0596303c8..7d359a5ab 100644
--- a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp
+++ b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.cpp
@@ -69,9 +69,7 @@ namespace storm {
             
             template <typename SparseModelType>
             SparseMultiObjectivePreprocessor<SparseModelType>::PreprocessorData::PreprocessorData(SparseModelType const& model) : originalModel(model) {
-                storm::storage::MemoryStructureBuilder<ValueType, RewardModelType> memoryBuilder(1, model);
-                memoryBuilder.setTransition(0,0, storm::storage::BitVector(model.getNumberOfStates(), true));
-                memory = std::make_shared<storm::storage::MemoryStructure>(memoryBuilder.build());
+                memory = std::make_shared<storm::storage::MemoryStructure>(storm::storage::MemoryStructureBuilder<ValueType, RewardModelType>::buildTrivialMemoryStructure(model));
                 
                 // The memoryLabelPrefix should not be a prefix of a state label of the given model to ensure uniqueness of label names
                 memoryLabelPrefix = "mem";
@@ -115,80 +113,81 @@ namespace storm {
                 
                 Objective<ValueType>& objective = *data.objectives.back();
                 
-                objective.considersComplementaryEvent = false;
+                // Check whether the complementary event is considered
+                objective.considersComplementaryEvent = formula.isProbabilityOperatorFormula() && formula.getSubformula().isGloballyFormula();
                 
+                storm::logic::OperatorInformation opInfo;
                 if (formula.hasBound()) {
-                    STORM_LOG_THROW(!formula.getBound().threshold.containsVariables(), storm::exceptions::InvalidPropertyException, "The formula " << formula << "considers a non-constant threshold");
-                    objective.bound = formula.getBound();
-                    if (storm::logic::isLowerBound(formula.getBound().comparisonType)) {
-                        objective.optimizationDirection = storm::solver::OptimizationDirection::Maximize;
-                    } else {
-                        objective.optimizationDirection = storm::solver::OptimizationDirection::Minimize;
-                    }
-                    STORM_LOG_WARN_COND(!formula.hasOptimalityType() || formula.getOptimalityType() == objective.optimizationDirection, "Optimization direction of formula " << formula << " ignored as the formula also specifies a threshold.");
-                } else if (formula.hasOptimalityType()){
-                    objective.optimizationDirection = formula.getOptimalityType();
-                } else {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Current objective " << formula << " does not specify whether to minimize or maximize");
-                }
-                
-                if (formula.isProbabilityOperatorFormula()){
-                    preprocessProbabilityOperatorFormula(formula.asProbabilityOperatorFormula(), data);
-                } else if (formula.isRewardOperatorFormula()){
-                    preprocessRewardOperatorFormula(formula.asRewardOperatorFormula(), data);
-                } else if (formula.isTimeOperatorFormula()){
-                    preprocessTimeOperatorFormula(formula.asTimeOperatorFormula(), data);
-                } else {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Could not preprocess the objective " << formula << " because it is not supported");
-                }
-                
-                // Invert the bound and optimization direction (if necessary)
-                if (objective.considersComplementaryEvent) {
-                    if (objective.bound) {
-                        objective.bound->threshold = objective.bound->threshold.getManager().rational(storm::utility::one<storm::RationalNumber>()) - objective.bound->threshold;
-                        switch (objective.bound->comparisonType) {
+                    opInfo.bound = formula.getBound();
+                    // Invert the bound (if necessary)
+                    if (objective.considersComplementaryEvent) {
+                        opInfo.bound->threshold = opInfo.bound->threshold.getManager().rational(storm::utility::one<storm::RationalNumber>()) - opInfo.bound->threshold;
+                        switch (opInfo.bound->comparisonType) {
                             case storm::logic::ComparisonType::Greater:
-                                objective.bound->comparisonType = storm::logic::ComparisonType::Less;
+                                opInfo.bound->comparisonType = storm::logic::ComparisonType::Less;
                                 break;
                             case storm::logic::ComparisonType::GreaterEqual:
-                                objective.bound->comparisonType = storm::logic::ComparisonType::LessEqual;
+                                opInfo.bound->comparisonType = storm::logic::ComparisonType::LessEqual;
                                 break;
                             case storm::logic::ComparisonType::Less:
-                                objective.bound->comparisonType = storm::logic::ComparisonType::Greater;
+                                opInfo.bound->comparisonType = storm::logic::ComparisonType::Greater;
                                 break;
                             case storm::logic::ComparisonType::LessEqual:
-                                objective.bound->comparisonType = storm::logic::ComparisonType::GreaterEqual;
+                                opInfo.bound->comparisonType = storm::logic::ComparisonType::GreaterEqual;
                                 break;
                             default:
                                 STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Current objective " << formula << " has unexpected comparison type");
                         }
                     }
-                    objective.optimizationDirection = storm::solver::invert(objective.optimizationDirection);
+                    if (storm::logic::isLowerBound(opInfo.bound->comparisonType)) {
+                        opInfo.optimalityType = storm::solver::OptimizationDirection::Maximize;
+                    } else {
+                        opInfo.optimalityType = storm::solver::OptimizationDirection::Minimize;
+                    }
+                    STORM_LOG_WARN_COND(!formula.hasOptimalityType(), "Optimization direction of formula " << formula << " ignored as the formula also specifies a threshold.");
+                } else if (formula.hasOptimalityType()){
+                    opInfo.optimalityType = formula.getOptimalityType();
+                    // Invert the optimality type (if necessary)
+                    if (objective.considersComplementaryEvent) {
+                        opInfo.optimalityType = storm::solver::invert(opInfo.optimalityType.get());
+                    }
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Objective " << formula << " does not specify whether to minimize or maximize");
+                }
+                
+                if (formula.isProbabilityOperatorFormula()){
+                    preprocessProbabilityOperatorFormula(formula.asProbabilityOperatorFormula(), opInfo, data);
+                } else if (formula.isRewardOperatorFormula()){
+                    preprocessRewardOperatorFormula(formula.asRewardOperatorFormula(), opInfo, data);
+                } else if (formula.isTimeOperatorFormula()){
+                    preprocessTimeOperatorFormula(formula.asTimeOperatorFormula(), opInfo, data);
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Could not preprocess the objective " << formula << " because it is not supported");
                 }
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessProbabilityOperatorFormula(storm::logic::ProbabilityOperatorFormula const& formula, PreprocessorData& data) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessProbabilityOperatorFormula(storm::logic::ProbabilityOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
                 
                 // Probabilities are between zero and one
                 data.objectives.back()->lowerResultBound = storm::utility::zero<ValueType>();
                 data.objectives.back()->upperResultBound = storm::utility::one<ValueType>();
                 
                 if (formula.getSubformula().isUntilFormula()){
-                    preprocessUntilFormula(formula.getSubformula().asUntilFormula(), data);
+                    preprocessUntilFormula(formula.getSubformula().asUntilFormula(), opInfo, data);
                 } else if (formula.getSubformula().isBoundedUntilFormula()){
-                    preprocessBoundedUntilFormula(formula.getSubformula().asBoundedUntilFormula(), data);
+                    preprocessBoundedUntilFormula(formula.getSubformula().asBoundedUntilFormula(), opInfo, data);
                 } else if (formula.getSubformula().isGloballyFormula()){
-                    preprocessGloballyFormula(formula.getSubformula().asGloballyFormula(), data);
+                    preprocessGloballyFormula(formula.getSubformula().asGloballyFormula(), opInfo, data);
                 } else if (formula.getSubformula().isEventuallyFormula()){
-                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), data);
+                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), opInfo, data);
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported.");
                 }
             }
 
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessRewardOperatorFormula(storm::logic::RewardOperatorFormula const& formula, PreprocessorData& data) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessRewardOperatorFormula(storm::logic::RewardOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
                 
                 STORM_LOG_THROW((formula.hasRewardModelName() && data.originalModel.hasRewardModel(formula.getRewardModelName()))
                                 || (!formula.hasRewardModelName() && data.originalModel.hasUniqueRewardModel()), storm::exceptions::InvalidPropertyException, "The reward model is not unique or the formula " << formula << " does not specify an existing reward model.");
@@ -205,43 +204,39 @@ namespace storm {
                 data.objectives.back()->lowerResultBound = storm::utility::zero<ValueType>();
                 
                 if (formula.getSubformula().isEventuallyFormula()){
-                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), data, rewardModelName);
+                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), opInfo, data, rewardModelName);
                 } else if (formula.getSubformula().isCumulativeRewardFormula()) {
-                    preprocessCumulativeRewardFormula(formula.getSubformula().asCumulativeRewardFormula(), data, rewardModelName);
+                    preprocessCumulativeRewardFormula(formula.getSubformula().asCumulativeRewardFormula(), opInfo, data, rewardModelName);
                 } else if (formula.getSubformula().isTotalRewardFormula()) {
-                    preprocessTotalRewardFormula(data, rewardModelName);
+                    preprocessTotalRewardFormula(opInfo, data, rewardModelName);
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported.");
                 }
             }
 
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessTimeOperatorFormula(storm::logic::TimeOperatorFormula const& formula, PreprocessorData& data) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessTimeOperatorFormula(storm::logic::TimeOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
                 // Time formulas are only supported for Markov automata
                 STORM_LOG_THROW(data.originalModel.isOfType(storm::models::ModelType::MarkovAutomaton), storm::exceptions::InvalidPropertyException, "Time operator formulas are only supported for Markov automata.");
                 
                 data.objectives.back()->lowerResultBound = storm::utility::zero<ValueType>();
                 
                 if (formula.getSubformula().isEventuallyFormula()){
-                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), data);
+                    preprocessEventuallyFormula(formula.getSubformula().asEventuallyFormula(), opInfo, data);
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The subformula of " << formula << " is not supported.");
                 }
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessUntilFormula(storm::logic::UntilFormula const& formula, PreprocessorData& data) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessUntilFormula(storm::logic::UntilFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, std::shared_ptr<storm::logic::Formula const> subformula) {
                 
                 storm::modelchecker::SparsePropositionalModelChecker<SparseModelType> mc(data.originalModel);
                 storm::storage::BitVector rightSubformulaResult = mc.check(formula.getRightSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector();
                 storm::storage::BitVector leftSubformulaResult = mc.check(formula.getLeftSubformula())->asExplicitQualitativeCheckResult().getTruthValuesVector();
                 // Check if the formula is already satisfied in the initial state because then the transformation to expected rewards will fail.
-                if (!data.objectives.back()->lowerTimeBound) {
-                    if (!(data.originalModel.getInitialStates() & rightSubformulaResult).empty()) {
-                        // TODO: Handle this case more properly
-                        STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "The Probability for the objective " << *data.objectives.back()->originalFormula << " is always one as the rhs of the until formula is true in the initial state. This (trivial) case is currently not implemented.");
-                    }
-                }
+                // TODO: Handle this case more properly
+                STORM_LOG_THROW((data.originalModel.getInitialStates() & rightSubformulaResult).empty(), storm::exceptions::NotImplementedException, "The Probability for the objective " << *data.objectives.back()->originalFormula << " is always one as the rhs of the until formula is true in the initial state. This (trivial) case is currently not implemented.");
                 
                 // Create a memory structure that stores whether a non-PhiState or a PsiState has already been reached
                 storm::storage::MemoryStructureBuilder<ValueType, RewardModelType> builder(2, data.originalModel);
@@ -257,7 +252,11 @@ namespace storm {
                 storm::storage::MemoryStructure objectiveMemory = builder.build();
                 data.memory = std::make_shared<storm::storage::MemoryStructure>(data.memory->product(objectiveMemory));
                 
-                data.objectives.back()->rewardModelName = data.rewardModelNamePrefix + std::to_string(data.objectives.size());
+                std::string rewardModelName = data.rewardModelNamePrefix + std::to_string(data.objectives.size());
+                if (subformula == nullptr) {
+                    subformula = std::make_shared<storm::logic::TotalRewardFormula>();
+                }
+                data.objectives.back()->formula = std::make_shared<storm::logic::RewardOperatorFormula>(subformula, rewardModelName, opInfo);
                 
                 auto relevantStatesFormula = std::make_shared<storm::logic::AtomicLabelFormula>(relevantStatesLabel);
                 data.tasks.push_back(std::make_shared<SparseMultiObjectivePreprocessorReachProbToTotalRewTask<SparseModelType>>(data.objectives.back(), relevantStatesFormula, formula.getRightSubformula().asSharedPointer()));
@@ -265,38 +264,37 @@ namespace storm {
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessBoundedUntilFormula(storm::logic::BoundedUntilFormula const& formula, PreprocessorData& data) {
-                if (formula.hasLowerBound()) {
-                    STORM_LOG_THROW(!formula.getLowerBound().containsVariables(), storm::exceptions::InvalidPropertyException, "The lower time bound for the formula " << formula << " still contains variables");
-                    if (!storm::utility::isZero(formula.getLowerBound<double>()) || formula.isLowerBoundStrict()) {
-                        data.objectives.back()->lowerTimeBound = storm::logic::TimeBound(formula.isLowerBoundStrict(), formula.getLowerBound());
-                    }
-                }
-                if (formula.hasUpperBound()) {
-                    STORM_LOG_THROW(!formula.getUpperBound().containsVariables(), storm::exceptions::InvalidPropertyException, "The Upper time bound for the formula " << formula << " still contains variables");
-                    if (!storm::utility::isInfinity(formula.getUpperBound<double>())) {
-                        data.objectives.back()->upperTimeBound = storm::logic::TimeBound(formula.isUpperBoundStrict(), formula.getUpperBound());
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessBoundedUntilFormula(storm::logic::BoundedUntilFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
+                
+                // Check how to handle this query
+                if (!formula.getTimeBoundReference().isRewardBound() && (!formula.hasLowerBound() || (!formula.isLowerBoundStrict() && storm::utility::isZero(formula.template getLowerBound<storm::RationalNumber>())))) {
+                    std::shared_ptr<storm::logic::Formula const> subformula;
+                    if (!formula.hasUpperBound()) {
+                        // The formula is actually unbounded
+                        subformula = std::make_shared<storm::logic::TotalRewardFormula>();
+                    } else {
+                        STORM_LOG_THROW(!data.originalModel.isOfType(storm::models::ModelType::MarkovAutomaton) || formula.getTimeBoundReference().isTimeBound(), storm::exceptions::InvalidPropertyException, "Bounded until formulas for Markov Automata are only allowed when time bounds are considered.");
+                        storm::logic::TimeBound bound(formula.isUpperBoundStrict(), formula.getUpperBound());
+                        subformula = std::make_shared<storm::logic::CumulativeRewardFormula>(bound, formula.getTimeBoundReference().getType());
                     }
+                    preprocessUntilFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), opInfo, data, subformula);
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "Property " << formula << "is not supported");
                 }
-                
-                data.objectives.back()->timeBoundReference = formula.getTimeBoundReference();
-                preprocessUntilFormula(storm::logic::UntilFormula(formula.getLeftSubformula().asSharedPointer(), formula.getRightSubformula().asSharedPointer()), data);
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessGloballyFormula(storm::logic::GloballyFormula const& formula, PreprocessorData& data) {
-                // The formula will be transformed to an until formula for the complementary event.
-                data.objectives.back()->considersComplementaryEvent = true;
-                
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessGloballyFormula(storm::logic::GloballyFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data) {
+                // The formula is transformed to an until formula for the complementary event.
                 auto negatedSubformula = std::make_shared<storm::logic::UnaryBooleanStateFormula>(storm::logic::UnaryBooleanStateFormula::OperatorType::Not, formula.getSubformula().asSharedPointer());
                 
-                preprocessUntilFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), data);
+                preprocessUntilFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), negatedSubformula), opInfo, data);
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessEventuallyFormula(storm::logic::EventuallyFormula const& formula, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessEventuallyFormula(storm::logic::EventuallyFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
                 if (formula.isReachabilityProbabilityFormula()){
-                    preprocessUntilFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), data);
+                    preprocessUntilFormula(storm::logic::UntilFormula(storm::logic::Formula::getTrueFormula(), formula.getSubformula().asSharedPointer()), opInfo, data);
                     return;
                 }
                 
@@ -320,7 +318,9 @@ namespace storm {
                 
                 auto relevantStatesFormula = std::make_shared<storm::logic::AtomicLabelFormula>(relevantStatesLabel);
                 
-                data.objectives.back()->rewardModelName = data.rewardModelNamePrefix + std::to_string(data.objectives.size());
+                std::string auxRewardModelName = data.rewardModelNamePrefix + std::to_string(data.objectives.size());
+                auto totalRewardFormula = std::make_shared<storm::logic::TotalRewardFormula>();
+                data.objectives.back()->formula = std::make_shared<storm::logic::RewardOperatorFormula>(totalRewardFormula, auxRewardModelName, opInfo);
                 data.finiteRewardCheckObjectives.set(data.objectives.size() - 1, true);
                 
                 if (formula.isReachabilityRewardFormula()) {
@@ -336,23 +336,19 @@ namespace storm {
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessCumulativeRewardFormula(storm::logic::CumulativeRewardFormula const& formula, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessCumulativeRewardFormula(storm::logic::CumulativeRewardFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
                 STORM_LOG_THROW(data.originalModel.isOfType(storm::models::ModelType::Mdp), storm::exceptions::InvalidPropertyException, "Cumulative reward formulas are not supported for the given model type.");
                 
-                STORM_LOG_THROW(!formula.getBound().containsVariables(), storm::exceptions::InvalidPropertyException, "The time bound for the formula " << formula << " still contains variables");
-                if (!storm::utility::isInfinity(formula.getBound<double>())) {
-                    data.objectives.back()->upperTimeBound = storm::logic::TimeBound(formula.isBoundStrict(), formula.getBound());
-                }
-
-                assert(optionalRewardModelName.is_initialized());
-                data.objectives.back()->rewardModelName = *optionalRewardModelName;
-
+                storm::logic::TimeBound bound(formula.isBoundStrict(), formula.getBound());
+                auto cumulativeRewardFormula = std::make_shared<storm::logic::CumulativeRewardFormula>(bound, storm::logic::TimeBoundType::Steps);
+                data.objectives.back()->formula = std::make_shared<storm::logic::RewardOperatorFormula>(cumulativeRewardFormula, *optionalRewardModelName, opInfo);
             }
             
             template<typename SparseModelType>
-            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessTotalRewardFormula(PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
-                assert(optionalRewardModelName.is_initialized());
-                data.objectives.back()->rewardModelName = *optionalRewardModelName;
+            void SparseMultiObjectivePreprocessor<SparseModelType>::preprocessTotalRewardFormula(storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName) {
+                
+                auto totalRewardFormula = std::make_shared<storm::logic::TotalRewardFormula>();
+                data.objectives.back()->formula = std::make_shared<storm::logic::RewardOperatorFormula>(totalRewardFormula, *optionalRewardModelName, opInfo);
                 data.finiteRewardCheckObjectives.set(data.objectives.size() - 1, true);
             }
             
@@ -374,9 +370,10 @@ namespace storm {
                 
                 std::set<std::string> relevantRewardModels;
                 for (auto const& obj : result.objectives) {
-                    relevantRewardModels.insert(*obj.rewardModelName);
-                    if (obj.timeBoundReference && obj.timeBoundReference->isRewardBound()) {
-                        relevantRewardModels.insert(obj.timeBoundReference->getRewardName());
+                    if (obj.formula->isRewardOperatorFormula()) {
+                        relevantRewardModels.insert(obj.formula->asRewardOperatorFormula().getRewardModelName());
+                    } else {
+                        STORM_LOG_ASSERT(false, "Unknown formula type.");
                     }
                 }
                 
@@ -409,7 +406,7 @@ namespace storm {
             typename SparseMultiObjectivePreprocessor<SparseModelType>::ReturnType::QueryType SparseMultiObjectivePreprocessor<SparseModelType>::getQueryType(std::vector<Objective<ValueType>> const& objectives) {
                 uint_fast64_t numOfObjectivesWithThreshold = 0;
                 for (auto& obj : objectives) {
-                    if (obj.bound) {
+                    if (obj.formula->hasBound()) {
                         ++numOfObjectivesWithThreshold;
                     }
                 }
@@ -436,8 +433,12 @@ namespace storm {
                 // Get the choices that yield non-zero reward
                 storm::storage::BitVector zeroRewardChoices(result.preprocessedModel->getNumberOfChoices(), true);
                 for (auto const& obj : result.objectives) {
-                    auto const& rewModel = result.preprocessedModel->getRewardModel(*obj.rewardModelName);
-                    zeroRewardChoices &= rewModel.getChoicesWithZeroReward(transitions);
+                    if (obj.formula->isRewardOperatorFormula()) {
+                        auto const& rewModel = result.preprocessedModel->getRewardModel(obj.formula->asRewardOperatorFormula().getRewardModelName());
+                        zeroRewardChoices &= rewModel.getChoicesWithZeroReward(transitions);
+                    } else {
+                        STORM_LOG_ASSERT(false, "Unknown formula type.");
+                    }
                 }
                 
                 // Get the states that have reward for at least one (or for all) choices assigned to it.
@@ -478,8 +479,9 @@ namespace storm {
                 storm::storage::BitVector maxRewardsToCheck(result.preprocessedModel->getNumberOfChoices(), true);
                 bool hasMinRewardToCheck = false;
                 for (auto const& objIndex : finiteRewardCheckObjectives) {
-                    auto const& rewModel = result.preprocessedModel->getRewardModel(result.objectives[objIndex].rewardModelName.get());
-                    if (storm::solver::minimize(result.objectives[objIndex].optimizationDirection)) {
+                    STORM_LOG_ASSERT(result.objectives[objIndex].formula->isRewardOperatorFormula(), "Objective needs to be checked for finite reward but has no reward operator.");
+                    auto const& rewModel = result.preprocessedModel->getRewardModel(result.objectives[objIndex].formula->asRewardOperatorFormula().getRewardModelName());
+                    if (storm::solver::minimize(result.objectives[objIndex].formula->getOptimalityType())) {
                         hasMinRewardToCheck = true;
                     } else {
                         maxRewardsToCheck &= rewModel.getChoicesWithZeroReward(transitions);
diff --git a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.h b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.h
index 318433575..505431e62 100644
--- a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.h
+++ b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessor.h
@@ -50,19 +50,20 @@ namespace storm {
                 /*!
                  * Apply the neccessary preprocessing for the given formula.
                  * @param formula the current (sub)formula
+                 * @param opInfo the information of the resulting operator formula
                  * @param data the current data. The currently processed objective is located at data.objectives.back()
                  * @param optionalRewardModelName the reward model name that is considered for the formula (if available)
                  */
                 static void preprocessOperatorFormula(storm::logic::OperatorFormula const& formula, PreprocessorData& data);
-                static void preprocessProbabilityOperatorFormula(storm::logic::ProbabilityOperatorFormula const& formula, PreprocessorData& data);
-                static void preprocessRewardOperatorFormula(storm::logic::RewardOperatorFormula const& formula, PreprocessorData& data);
-                static void preprocessTimeOperatorFormula(storm::logic::TimeOperatorFormula const& formula, PreprocessorData& data);
-                static void preprocessUntilFormula(storm::logic::UntilFormula const& formula, PreprocessorData& data);
-                static void preprocessBoundedUntilFormula(storm::logic::BoundedUntilFormula const& formula, PreprocessorData& data);
-                static void preprocessGloballyFormula(storm::logic::GloballyFormula const& formula, PreprocessorData& data);
-                static void preprocessEventuallyFormula(storm::logic::EventuallyFormula const& formula, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none);
-                static void preprocessCumulativeRewardFormula(storm::logic::CumulativeRewardFormula const& formula, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none);
-                static void preprocessTotalRewardFormula(PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none); // The total reward formula itself does not need to be provided as it is unique.
+                static void preprocessProbabilityOperatorFormula(storm::logic::ProbabilityOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data);
+                static void preprocessRewardOperatorFormula(storm::logic::RewardOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data);
+                static void preprocessTimeOperatorFormula(storm::logic::TimeOperatorFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data);
+                static void preprocessUntilFormula(storm::logic::UntilFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, std::shared_ptr<storm::logic::Formula const> subformula = nullptr);
+                static void preprocessBoundedUntilFormula(storm::logic::BoundedUntilFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data);
+                static void preprocessGloballyFormula(storm::logic::GloballyFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data);
+                static void preprocessEventuallyFormula(storm::logic::EventuallyFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none);
+                static void preprocessCumulativeRewardFormula(storm::logic::CumulativeRewardFormula const& formula, storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none);
+                static void preprocessTotalRewardFormula(storm::logic::OperatorInformation const& opInfo, PreprocessorData& data, boost::optional<std::string> const& optionalRewardModelName = boost::none); // The total reward formula itself does not need to be provided as it is unique.
                 
                 
                 /*!
diff --git a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessorTask.h b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessorTask.h
index e34e37584..f84b935b6 100644
--- a/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessorTask.h
+++ b/src/storm/modelchecker/multiobjective/SparseMultiObjectivePreprocessorTask.h
@@ -58,8 +58,9 @@ namespace storm {
                             objectiveRewards[row] = preprocessedModel.getTransitionMatrix().getConstrainedRowSum(row, goalStates);
                         }
                     }
-                    STORM_LOG_ASSERT(this->objective->rewardModelName.is_initialized(), "No reward model name has been specified");
-                    preprocessedModel.addRewardModel(this->objective->rewardModelName.get(), typename SparseModelType::RewardModelType(boost::none, std::move(objectiveRewards)));
+                    STORM_LOG_ASSERT(this->objective->formula->isRewardOperatorFormula(), "No reward operator formula.");
+                    STORM_LOG_ASSERT(this->objective->formula->asRewardOperatorFormula().hasRewardModelName(), "No reward model name has been specified");
+                    preprocessedModel.addRewardModel(this->objective->formula->asRewardOperatorFormula().getRewardModelName(), typename SparseModelType::RewardModelType(boost::none, std::move(objectiveRewards)));
                 }
                         
             private:
@@ -91,8 +92,9 @@ namespace storm {
                             std::fill_n(objectiveRewards.getStateActionRewardVector().begin() + preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state], preprocessedModel.getTransitionMatrix().getRowGroupSize(state), storm::utility::zero<typename SparseModelType::ValueType>());
                         }
                     }
-                    STORM_LOG_ASSERT(this->objective->rewardModelName.is_initialized(), "No reward model name has been specified");
-                    preprocessedModel.addRewardModel(this->objective->rewardModelName.get(), std::move(objectiveRewards));
+                    STORM_LOG_ASSERT(this->objective->formula->isRewardOperatorFormula(), "No reward operator formula.");
+                    STORM_LOG_ASSERT(this->objective->formula->asRewardOperatorFormula().hasRewardModelName(), "No reward model name has been specified");
+                    preprocessedModel.addRewardModel(this->objective->formula->asRewardOperatorFormula().getRewardModelName(), std::move(objectiveRewards));
                 }
                         
             private:
@@ -117,8 +119,10 @@ namespace storm {
                     
                     std::vector<typename SparseModelType::ValueType> timeRewards(preprocessedModel.getNumberOfStates(), storm::utility::zero<typename SparseModelType::ValueType>());
                     storm::utility::vector::setVectorValues(timeRewards, dynamic_cast<storm::models::sparse::MarkovAutomaton<typename SparseModelType::ValueType> const&>(preprocessedModel).getMarkovianStates() & relevantStates, storm::utility::one<typename SparseModelType::ValueType>());
-                    STORM_LOG_ASSERT(this->objective->rewardModelName.is_initialized(), "No reward model name has been specified");
-                    preprocessedModel.addRewardModel(this->objective->rewardModelName.get(), typename SparseModelType::RewardModelType(std::move(timeRewards)));
+
+                    STORM_LOG_ASSERT(this->objective->formula->isRewardOperatorFormula(), "No reward operator formula.");
+                    STORM_LOG_ASSERT(this->objective->formula->asRewardOperatorFormula().hasRewardModelName(), "No reward model name has been specified");
+                    preprocessedModel.addRewardModel(this->objective->formula->asRewardOperatorFormula().getRewardModelName(), typename SparseModelType::RewardModelType(std::move(timeRewards)));
                 }
                         
             private:
diff --git a/src/storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.cpp b/src/storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.cpp
index 86faf88a7..628fbb687 100644
--- a/src/storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.cpp
+++ b/src/storm/modelchecker/multiobjective/constraintbased/SparseCbAchievabilityQuery.cpp
@@ -26,7 +26,7 @@ namespace storm {
             
             template <class SparseModelType>
             SparseCbAchievabilityQuery<SparseModelType>::SparseCbAchievabilityQuery(SparseMultiObjectivePreprocessorReturnType<SparseModelType>& preprocessorResult) : SparseCbQuery<SparseModelType>(preprocessorResult) {
-                STORM_LOG_ASSERT(preprocessorResult.queryType==SparseMultiObjectivePreprocessorReturnType<SparseModelType>::QueryType::Achievability, "Invalid query Type");
+                STORM_LOG_ASSERT(preprocessorResult.queryType == SparseMultiObjectivePreprocessorReturnType<SparseModelType>::QueryType::Achievability, "Invalid query Type");
                 solver = storm::utility::solver::SmtSolverFactory().create(*this->expressionManager);
             }
 
@@ -134,36 +134,34 @@ namespace storm {
             void SparseCbAchievabilityQuery<SparseModelType>::addObjectiveConstraints() {
                 storm::expressions::Expression zero = this->expressionManager->rational(storm::utility::zero<ValueType>());
                 for (Objective<ValueType> const& obj : this->objectives) {
-                    if (obj.rewardModelName) {
-                        STORM_LOG_THROW(obj.bound, storm::exceptions::InvalidOperationException, "Invoked achievability query but no bound was specified for at least one objective.");
-                        STORM_LOG_THROW(!obj.lowerTimeBound && !obj.upperTimeBound, storm::exceptions::NotSupportedException, "Constraint based method currently does not support step bounds");
-                        std::vector<ValueType> rewards = getActionBasedExpectedRewards(*obj.rewardModelName);
-                        storm::expressions::Expression objValue = zero;
-                        for (uint_fast64_t choice = 0; choice < rewards.size(); ++choice) {
-                            if (!storm::utility::isZero(rewards[choice])) {
-                                objValue = objValue + (this->expressionManager->rational(rewards[choice]) * expectedChoiceVariables[choice].getExpression());
-                            }
-                        }
-                        // We need to actually evaluate the threshold as rational number. Otherwise a threshold like '<=16/9' might be considered as 1 due to integer division
-                        STORM_LOG_THROW(!obj.bound->threshold.containsVariables(), storm::exceptions::InvalidOperationException, "The threshold for one objective still contains undefined variables");
-                        storm::expressions::Expression threshold = this->expressionManager->rational(obj.bound->threshold.evaluateAsRational());
-                        switch (obj.bound->comparisonType) {
-                            case storm::logic::ComparisonType::Greater:
-                                solver->add( objValue > threshold);
-                                break;
-                            case storm::logic::ComparisonType::GreaterEqual:
-                                solver->add( objValue >= threshold);
-                                break;
-                            case storm::logic::ComparisonType::Less:
-                                solver->add( objValue < threshold);
-                                break;
-                            case storm::logic::ComparisonType::LessEqual:
-                                solver->add( objValue <= threshold);
-                                break;
-                            default:
-                                STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "One or more objectives have an invalid comparison type");
+                    STORM_LOG_THROW(obj.formula->isRewardOperatorFormula() && obj.formula->getSubformula().isTotalRewardFormula(), storm::exceptions::InvalidOperationException, "Constraint-based solver only supports total-reward objectives. Got " << *obj.formula << " instead.");
+                    STORM_LOG_THROW(obj.formula->hasBound(), storm::exceptions::InvalidOperationException, "Invoked achievability query but no bound was specified for at least one objective.");
+                    STORM_LOG_THROW(obj.formula->asRewardOperatorFormula().hasRewardModelName(), storm::exceptions::InvalidOperationException, "Expected reward operator with a reward model name. Got " << *obj.formula << " instead.");
+                    std::vector<ValueType> rewards = getActionBasedExpectedRewards(obj.formula->asRewardOperatorFormula().getRewardModelName());
+                    storm::expressions::Expression objValue = zero;
+                    for (uint_fast64_t choice = 0; choice < rewards.size(); ++choice) {
+                        if (!storm::utility::isZero(rewards[choice])) {
+                            objValue = objValue + (this->expressionManager->rational(rewards[choice]) * expectedChoiceVariables[choice].getExpression());
                         }
                     }
+                    // We need to actually evaluate the threshold as rational number. Otherwise a threshold like '<=16/9' might be considered as 1 due to integer division
+                    storm::expressions::Expression threshold = this->expressionManager->rational(obj.formula->getThreshold().evaluateAsRational());
+                    switch (obj.formula->getBound().comparisonType) {
+                        case storm::logic::ComparisonType::Greater:
+                            solver->add( objValue > threshold);
+                            break;
+                        case storm::logic::ComparisonType::GreaterEqual:
+                            solver->add( objValue >= threshold);
+                            break;
+                        case storm::logic::ComparisonType::Less:
+                            solver->add( objValue < threshold);
+                            break;
+                        case storm::logic::ComparisonType::LessEqual:
+                            solver->add( objValue <= threshold);
+                            break;
+                        default:
+                            STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "One or more objectives have an invalid comparison type");
+                    }
                 }
             }
             
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp
index e65aa6325..e695fa8db 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMaPcaaWeightVectorChecker.cpp
@@ -8,9 +8,11 @@
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
 #include "storm/solver/GmmxxLinearEquationSolver.h"
+#include "storm/logic/Formulas.h"
 
 #include "storm/exceptions/InvalidOperationException.h"
 #include "storm/exceptions/InvalidPropertyException.h"
+#include "storm/exceptions/UnexpectedException.h"
 
 namespace storm {
     namespace modelchecker {
@@ -25,13 +27,17 @@ namespace storm {
             SparsePcaaWeightVectorChecker<SparseMaModelType>(model, objectives, possibleECActions, possibleBottomStates) {
                 // Set the (discretized) state action rewards.
                 this->discreteActionRewards.resize(objectives.size());
-                for(auto objIndex : this->objectivesWithNoUpperTimeBound) {
-                    typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(*this->objectives[objIndex].rewardModelName);
+                for (auto objIndex : this->objectivesWithNoUpperTimeBound) {
+                    auto const& formula = *objectives[objIndex].formula;
+                    STORM_LOG_THROW(formula.isRewardOperatorFormula() && formula.asRewardOperatorFormula().hasRewardModelName(), storm::exceptions::UnexpectedException, "Unexpected type of operator formula: " << formula);
+                    STORM_LOG_THROW(formula.getSubformula().isTotalRewardFormula() || (formula.getSubformula().isCumulativeRewardFormula() && formula.getSubformula().asCumulativeRewardFormula().isTimeBounded()), storm::exceptions::UnexpectedException, "Unexpected type of sub-formula: " << formula.getSubformula());
+                    STORM_LOG_WARN_COND(!formula.getSubformula().isCumulativeRewardFormula() || (objectives[objIndex].originalFormula->isProbabilityOperatorFormula() && objectives[objIndex].originalFormula->asProbabilityOperatorFormula().getSubformula().isBoundedUntilFormula()), "Objective " << objectives[objIndex].originalFormula << " was simplified to a cumulative reward formula. Correctness of the algorithm is unknown for this type of property.");
+                    typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName());
                     STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected.");
                     this->discreteActionRewards[objIndex] = rewModel.hasStateActionRewards() ? rewModel.getStateActionRewardVector() : std::vector<ValueType>(this->model.getTransitionMatrix().getRowCount(), storm::utility::zero<ValueType>());
-                    if(rewModel.hasStateRewards()) {
+                    if (rewModel.hasStateRewards()) {
                         // Note that state rewards are earned over time and thus play no role for probabilistic states
-                        for(auto markovianState : this->model.getMarkovianStates()) {
+                        for (auto markovianState : this->model.getMarkovianStates()) {
                             this->discreteActionRewards[objIndex][this->model.getTransitionMatrix().getRowGroupIndices()[markovianState]] += rewModel.getStateReward(markovianState) / this->model.getExitRate(markovianState);
                         }
                     }
@@ -41,10 +47,6 @@ namespace storm {
             template <class SparseMaModelType>
             void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::boundedPhase(std::vector<ValueType> const& weightVector, std::vector<ValueType>& weightedRewardVector) {
                 
-                for (auto const& obj : this->objectives) {
-                    STORM_LOG_THROW(!obj.timeBoundReference || obj.timeBoundReference->isTimeBound(), storm::exceptions::InvalidPropertyException, "Multi-objective model checking of Markov automata is only supported for time-bounded formulass.");
-                }
-                
                 // Split the preprocessed model into transitions from/to probabilistic/Markovian states.
                 SubModel MS = createSubModel(true, weightedRewardVector);
                 SubModel PS = createSubModel(false, weightedRewardVector);
@@ -82,7 +84,7 @@ namespace storm {
                     
                     // Compute values that can be obtained at Markovian states after letting one (digitized) time unit pass.
                     // Only perform such a step if there is time left.
-                    if(currentEpoch>0) {
+                    if (currentEpoch>0) {
                         performMSStep(MS, PS, consideredObjectives, weightVector);
                         --currentEpoch;
                     } else {
@@ -93,7 +95,7 @@ namespace storm {
                 // compose the results from MS and PS
                 storm::utility::vector::setVectorValues(this->weightedResult, MS.states, MS.weightedSolutionVector);
                 storm::utility::vector::setVectorValues(this->weightedResult, PS.states, PS.weightedSolutionVector);
-                for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], MS.states, MS.objectiveSolutionVectors[objIndex]);
                     storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], PS.states, PS.objectiveSolutionVectors[objIndex]);
                 }
@@ -118,16 +120,16 @@ namespace storm {
                 result.weightedRewardVector.resize(result.getNumberOfChoices());
                 storm::utility::vector::selectVectorValues(result.weightedRewardVector, result.choices, weightedRewardVector);
                 result.objectiveRewardVectors.resize(this->objectives.size());
-                for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     std::vector<ValueType>& objVector = result.objectiveRewardVectors[objIndex];
                     objVector = std::vector<ValueType>(result.weightedRewardVector.size(), storm::utility::zero<ValueType>());
-                    if(this->objectivesWithNoUpperTimeBound.get(objIndex)) {
+                    if (this->objectivesWithNoUpperTimeBound.get(objIndex)) {
                         storm::utility::vector::selectVectorValues(objVector, result.choices, this->discreteActionRewards[objIndex]);
                     } else {
-                        typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(*this->objectives[objIndex].rewardModelName);
+                       typename SparseMaModelType::RewardModelType const& rewModel = this->model.getRewardModel(this->objectives[objIndex].formula->asRewardOperatorFormula().getRewardModelName());
                         STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected.");
                         STORM_LOG_ASSERT(!rewModel.hasStateRewards(), "State rewards for bounded objectives for MAs are not expected (bounded rewards are not supported).");
-                        if(rewModel.hasStateActionRewards()) {
+                        if (rewModel.hasStateActionRewards()) {
                             storm::utility::vector::selectVectorValues(objVector, result.choices, rewModel.getStateActionRewardVector());
                         }
                     }
@@ -136,7 +138,7 @@ namespace storm {
                 result.weightedSolutionVector.resize(result.getNumberOfStates());
                 storm::utility::vector::selectVectorValues(result.weightedSolutionVector, result.states, this->weightedResult);
                 result.objectiveSolutionVectors.resize(this->objectives.size());
-                for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     result.objectiveSolutionVectors[objIndex].resize(result.weightedSolutionVector.size());
                     storm::utility::vector::selectVectorValues(result.objectiveSolutionVectors[objIndex], result.states, this->objectiveResults[objIndex]);
                 }
@@ -163,11 +165,10 @@ namespace storm {
                 std::vector<VT> timeBounds;
                 std::vector<VT> eToPowerOfMinusMaxRateTimesBound;
                 VT smallestNonZeroBound = storm::utility::zero<VT>();
-                for(auto const& obj : this->objectives) {
-                    if(obj.upperTimeBound){
-                        STORM_LOG_THROW(!obj.upperTimeBound->getBound().containsVariables(), storm::exceptions::InvalidOperationException, "The time bound '" << obj.upperTimeBound->getBound() << " contains undefined variables");
-                        timeBounds.push_back(storm::utility::convertNumber<VT>(obj.upperTimeBound->getBound().evaluateAsRational()));
-                        STORM_LOG_ASSERT(!storm::utility::isZero(timeBounds.back()), "Got zero-valued upper time bound.");
+                for (auto const& obj : this->objectives) {
+                    if (obj.formula->getSubformula().isCumulativeRewardFormula()) {
+                        timeBounds.push_back(obj.formula->getSubformula().asCumulativeRewardFormula().template getBound<VT>());
+                        STORM_LOG_THROW(!storm::utility::isZero(timeBounds.back()), storm::exceptions::InvalidPropertyException, "Got zero-valued upper time bound. This is not suppoted.");
                         eToPowerOfMinusMaxRateTimesBound.push_back(std::exp(-maxRate * timeBounds.back()));
                         smallestNonZeroBound = storm::utility::isZero(smallestNonZeroBound) ? timeBounds.back() : std::min(smallestNonZeroBound, timeBounds.back());
                     } else {
@@ -175,7 +176,7 @@ namespace storm {
                         eToPowerOfMinusMaxRateTimesBound.push_back(storm::utility::zero<VT>());
                     }
                 }
-                if(storm::utility::isZero(smallestNonZeroBound)) {
+                if (storm::utility::isZero(smallestNonZeroBound)) {
                     // There are no time bounds. In this case, one is a valid digitization constant.
                     return storm::utility::one<VT>();
                 }
@@ -189,16 +190,16 @@ namespace storm {
                 VT delta = smallestNonZeroBound / smallestStepBound;
                 while(true) {
                     bool deltaValid = true;
-                    for(auto const& objIndex : objectivesWithTimeBound) {
+                    for (auto const& objIndex : objectivesWithTimeBound) {
                         auto const& timeBound = timeBounds[objIndex];
-                        if(timeBound/delta != std::floor(timeBound/delta)) {
+                        if (timeBound/delta != std::floor(timeBound/delta)) {
                             deltaValid = false;
                             break;
                         }
                     }
-                    if(deltaValid) {
+                    if (deltaValid) {
                         VT weightedPrecisionForCurrentDelta = storm::utility::zero<VT>();
-                        for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                        for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                             VT precisionOfObj = storm::utility::zero<VT>();
                             if (objectivesWithTimeBound.get(objIndex)) {
                                 precisionOfObj += storm::utility::one<VT>() - (eToPowerOfMinusMaxRateTimesBound[objIndex] * storm::utility::pow(storm::utility::one<VT>() + maxRate * delta, timeBounds[objIndex] / delta) );
@@ -207,7 +208,7 @@ namespace storm {
                         }
                         deltaValid &= weightedPrecisionForCurrentDelta <= goalPrecisionTimesNorm;
                     }
-                    if(deltaValid) {
+                    if (deltaValid) {
                         break;
                     }
                     ++smallestStepBound;
@@ -229,19 +230,19 @@ namespace storm {
             void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::digitize(SubModel& MS, VT const& digitizationConstant) const {
                 std::vector<VT> rateVector(MS.getNumberOfChoices());
                 storm::utility::vector::selectVectorValues(rateVector, MS.states, this->model.getExitRates());
-                for(uint_fast64_t row = 0; row < rateVector.size(); ++row) {
+                for (uint_fast64_t row = 0; row < rateVector.size(); ++row) {
                     VT const eToMinusRateTimesDelta = std::exp(-rateVector[row] * digitizationConstant);
-                    for(auto& entry : MS.toMS.getRow(row)) {
+                    for (auto& entry : MS.toMS.getRow(row)) {
                         entry.setValue((storm::utility::one<VT>() - eToMinusRateTimesDelta) * entry.getValue());
-                        if(entry.getColumn() == row) {
+                        if (entry.getColumn() == row) {
                             entry.setValue(entry.getValue() + eToMinusRateTimesDelta);
                         }
                     }
-                    for(auto& entry : MS.toPS.getRow(row)) {
+                    for (auto& entry : MS.toPS.getRow(row)) {
                         entry.setValue((storm::utility::one<VT>() - eToMinusRateTimesDelta) * entry.getValue());
                     }
                     MS.weightedRewardVector[row] *= storm::utility::one<VT>() - eToMinusRateTimesDelta;
-                    for(auto& objVector : MS.objectiveRewardVectors) {
+                    for (auto& objVector : MS.objectiveRewardVectors) {
                         objVector[row] *= storm::utility::one<VT>() - eToMinusRateTimesDelta;
                     }
                 }
@@ -258,13 +259,12 @@ namespace storm {
             void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::digitizeTimeBounds(TimeBoundMap& upperTimeBounds, VT const& digitizationConstant) {
                 
                 VT const maxRate = this->model.getMaximalExitRate();
-                for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     auto const& obj = this->objectives[objIndex];
-                    STORM_LOG_THROW(!obj.lowerTimeBound, storm::exceptions::InvalidPropertyException, "Lower time bounds are not supported by this model checker");
                     VT errorTowardsZero = storm::utility::zero<VT>();
                     VT errorAwayFromZero = storm::utility::zero<VT>();
-                    if(obj.upperTimeBound) {
-                        VT timeBound = storm::utility::convertNumber<VT>(obj.upperTimeBound->getBound().evaluateAsRational());
+                    if (obj.formula->getSubformula().isCumulativeRewardFormula()) {
+                        VT timeBound = obj.formula->getSubformula().asCumulativeRewardFormula().template getBound<VT>();
                         uint_fast64_t digitizedBound = storm::utility::convertNumber<uint_fast64_t>(timeBound/digitizationConstant);
                         auto timeBoundIt = upperTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->objectives.size(), false))).first;
                         timeBoundIt->second.set(objIndex);
@@ -272,7 +272,7 @@ namespace storm {
                         digitizationError -= std::exp(-maxRate * timeBound) * storm::utility::pow(storm::utility::one<VT>() + maxRate * digitizationConstant, digitizedBound);
                         errorAwayFromZero += digitizationError;
                     }
-                    if (storm::solver::maximize(obj.optimizationDirection)) {
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
                         this->offsetsToUnderApproximation[objIndex] = -errorTowardsZero;
                         this->offsetsToOverApproximation[objIndex] = errorAwayFromZero;
                     } else {
@@ -324,11 +324,11 @@ namespace storm {
             template <class SparseMaModelType>
             void SparseMaPcaaWeightVectorChecker<SparseMaModelType>::updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector<ValueType> const& weightVector, TimeBoundMap::iterator& upperTimeBoundIt, TimeBoundMap const& upperTimeBounds) {
                 
-                if(upperTimeBoundIt != upperTimeBounds.end() && currentEpoch == upperTimeBoundIt->first) {
+                if (upperTimeBoundIt != upperTimeBounds.end() && currentEpoch == upperTimeBoundIt->first) {
                     consideredObjectives |= upperTimeBoundIt->second;
-                    for(auto objIndex : upperTimeBoundIt->second) {
+                    for (auto objIndex : upperTimeBoundIt->second) {
                         // This objective now plays a role in the weighted sum
-                        ValueType factor = storm::solver::minimize(this->objectives[objIndex].optimizationDirection) ? -weightVector[objIndex] : weightVector[objIndex];
+                        ValueType factor = storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType()) ? -weightVector[objIndex] : weightVector[objIndex];
                         storm::utility::vector::addScaledVector(MS.weightedRewardVector, MS.objectiveRewardVectors[objIndex], factor);
                         storm::utility::vector::addScaledVector(PS.weightedRewardVector, PS.objectiveRewardVectors[objIndex], factor);
                     }
@@ -345,16 +345,16 @@ namespace storm {
                 // compute a choice vector for the probabilistic states that is optimal w.r.t. the weighted reward vector
                 minMax.solver->solveEquations(PS.weightedSolutionVector, minMax.b);
                 auto const& newChoices = minMax.solver->getSchedulerChoices();
-                if(consideredObjectives.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*consideredObjectives.begin()])) {
+                if (consideredObjectives.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*consideredObjectives.begin()])) {
                     // In this case there is no need to perform the computation on the individual objectives
                     optimalChoicesAtCurrentEpoch = newChoices;
                     PS.objectiveSolutionVectors[*consideredObjectives.begin()] = PS.weightedSolutionVector;
-                    if (storm::solver::minimize(this->objectives[*consideredObjectives.begin()].optimizationDirection)) {
+                    if (storm::solver::minimize(this->objectives[*consideredObjectives.begin()].formula->getOptimalityType())) {
                         storm::utility::vector::scaleVectorInPlace(PS.objectiveSolutionVectors[*consideredObjectives.begin()], -storm::utility::one<ValueType>());
                     }
                 } else {
                     // check whether the linEqSolver needs to be updated, i.e., whether the scheduler has changed
-                    if(linEq.solver == nullptr || newChoices != optimalChoicesAtCurrentEpoch) {
+                    if (linEq.solver == nullptr || newChoices != optimalChoicesAtCurrentEpoch) {
                         optimalChoicesAtCurrentEpoch = newChoices;
                         linEq.solver = nullptr;
                         storm::storage::SparseMatrix<ValueType> linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true);
@@ -365,17 +365,17 @@ namespace storm {
                     
                     // Get the results for the individual objectives.
                     // Note that we do not consider an estimate for each objective (as done in the unbounded phase) since the results from the previous epoch are already pretty close
-                    for(auto objIndex : consideredObjectives) {
+                    for (auto objIndex : consideredObjectives) {
                         auto const& objectiveRewardVectorPS = PS.objectiveRewardVectors[objIndex];
                         auto const& objectiveSolutionVectorMS = MS.objectiveSolutionVectors[objIndex];
                         // compute rhs of equation system, i.e., PS.toMS * x + Rewards
                         // To safe some time, only do this for the obtained optimal choices
                         auto itGroupIndex = PS.toPS.getRowGroupIndices().begin();
                         auto itChoiceOffset = optimalChoicesAtCurrentEpoch.begin();
-                        for(auto& bValue : linEq.b) {
+                        for (auto& bValue : linEq.b) {
                             uint_fast64_t row = (*itGroupIndex) + (*itChoiceOffset);
                             bValue = objectiveRewardVectorPS[row];
-                            for(auto const& entry : PS.toMS.getRow(row)){
+                            for (auto const& entry : PS.toMS.getRow(row)){
                                 bValue += entry.getValue() * objectiveSolutionVectorMS[entry.getColumn()];
                             }
                             ++itGroupIndex;
@@ -393,14 +393,14 @@ namespace storm {
                 storm::utility::vector::addVectors(MS.weightedRewardVector, MS.auxChoiceValues, MS.weightedSolutionVector);
                 MS.toPS.multiplyWithVector(PS.weightedSolutionVector, MS.auxChoiceValues);
                 storm::utility::vector::addVectors(MS.weightedSolutionVector, MS.auxChoiceValues, MS.weightedSolutionVector);
-                if(consideredObjectives.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*consideredObjectives.begin()])) {
+                if (consideredObjectives.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*consideredObjectives.begin()])) {
                     // In this case there is no need to perform the computation on the individual objectives
                     MS.objectiveSolutionVectors[*consideredObjectives.begin()] = MS.weightedSolutionVector;
-                    if (storm::solver::minimize(this->objectives[*consideredObjectives.begin()].optimizationDirection)) {
+                    if (storm::solver::minimize(this->objectives[*consideredObjectives.begin()].formula->getOptimalityType())) {
                         storm::utility::vector::scaleVectorInPlace(MS.objectiveSolutionVectors[*consideredObjectives.begin()], -storm::utility::one<ValueType>());
                     }
                 } else {
-                    for(auto objIndex : consideredObjectives) {
+                    for (auto objIndex : consideredObjectives) {
                         MS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues);
                         storm::utility::vector::addVectors(MS.objectiveRewardVectors[objIndex], MS.auxChoiceValues, MS.objectiveSolutionVectors[objIndex]);
                         MS.toPS.multiplyWithVector(PS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues);
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp
index 0916fd94c..0744f1d7d 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparseMdpPcaaWeightVectorChecker.cpp
@@ -5,6 +5,7 @@
 #include "storm/models/sparse/StandardRewardModel.h"
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
+#include "storm/logic/Formulas.h"
 #include "storm/exceptions/InvalidPropertyException.h"
 #include "storm/exceptions/IllegalArgumentException.h"
 #include "storm/exceptions/NotSupportedException.h"
@@ -21,24 +22,22 @@ namespace storm {
                                                                                                    storm::storage::BitVector const& possibleECActions,
                                                                                                    storm::storage::BitVector const& possibleBottomStates) :
                 SparsePcaaWeightVectorChecker<SparseMdpModelType>(model, objectives, possibleECActions, possibleBottomStates) {
-                // set the state action rewards
+                // set the state action rewards. Also do some sanity checks on the objectives.
                 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    typename SparseMdpModelType::RewardModelType const& rewModel = this->model.getRewardModel(*this->objectives[objIndex].rewardModelName);
-                    STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Reward model has transition rewards which is not expected.");
+                    auto const& formula = *objectives[objIndex].formula;
+                    STORM_LOG_THROW(formula.isRewardOperatorFormula() && formula.asRewardOperatorFormula().hasRewardModelName(), storm::exceptions::UnexpectedException, "Unexpected type of operator formula: " << formula);
+                    STORM_LOG_THROW(formula.getSubformula().isCumulativeRewardFormula() || formula.getSubformula().isTotalRewardFormula(), storm::exceptions::UnexpectedException, "Unexpected type of sub-formula: " << formula.getSubformula());
+                    typename SparseMdpModelType::RewardModelType const& rewModel = this->model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName());
+                    STORM_LOG_THROW(!rewModel.hasTransitionRewards(), storm::exceptions::NotSupportedException, "Reward model has transition rewards which is not expected.");
                     this->discreteActionRewards[objIndex] = rewModel.getTotalRewardVector(this->model.getTransitionMatrix());
                 }
             }
             
             template <class SparseMdpModelType>
             void SparseMdpPcaaWeightVectorChecker<SparseMdpModelType>::boundedPhase(std::vector<ValueType> const& weightVector, std::vector<ValueType>& weightedRewardVector) {
-                // Check whether reward bounded objectives occur.
+                // Currently, only step bounds are considered.
+                // TODO: Check whether reward bounded objectives occur.
                 bool containsRewardBoundedObjectives = false;
-                for (auto const& obj : this->objectives) {
-                    if (obj.timeBoundReference && obj.timeBoundReference->isRewardBound()) {
-                        containsRewardBoundedObjectives = true;
-                        break;
-                    }
-                }
                 
                 if (containsRewardBoundedObjectives) {
                     boundedPhaseWithRewardBounds(weightVector, weightedRewardVector);
@@ -56,12 +55,10 @@ namespace storm {
                 // Get for each occurring timeBound the indices of the objectives with that bound.
                 std::map<uint_fast64_t, storm::storage::BitVector, std::greater<uint_fast64_t>> stepBounds;
                 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    auto const& obj = this->objectives[objIndex];
-                    STORM_LOG_THROW(!obj.lowerTimeBound, storm::exceptions::InvalidPropertyException, "Lower step bounds are not supported by this model checker");
-                    if (obj.upperTimeBound) {
-                        STORM_LOG_THROW(!obj.upperTimeBound->getBound().containsVariables(), storm::exceptions::InvalidPropertyException, "The step bound '" << obj.upperTimeBound->getBound() << " contains undefined variables");
-                        uint_fast64_t stepBound = (uint_fast64_t) obj.upperTimeBound->getBound().evaluateAsInt();
-                        if (obj.upperTimeBound->isStrict()) {
+                    if (this->objectives[objIndex].formula->getSubformula().isCumulativeRewardFormula()) {
+                        auto const& subformula = this->objectives[objIndex].formula->getSubformula().asCumulativeRewardFormula();
+                        uint_fast64_t stepBound = subformula.template getBound<uint_fast64_t>();
+                        if (subformula.isBoundStrict()) {
                             --stepBound;
                         }
                         auto stepBoundIt = stepBounds.insert(std::make_pair(stepBound, storm::storage::BitVector(this->objectives.size(), false))).first;
@@ -85,7 +82,7 @@ namespace storm {
                         consideredObjectives |= stepBoundIt->second;
                         for(auto objIndex : stepBoundIt->second) {
                             // This objective now plays a role in the weighted sum
-                            ValueType factor = storm::solver::minimize(this->objectives[objIndex].optimizationDirection) ? -weightVector[objIndex] : weightVector[objIndex];
+                            ValueType factor = storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType()) ? -weightVector[objIndex] : weightVector[objIndex];
                             storm::utility::vector::addScaledVector(weightedRewardVector, this->discreteActionRewards[objIndex], factor);
                         }
                         ++stepBoundIt;
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp
index 74e646103..a1d457c11 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaAchievabilityQuery.cpp
@@ -31,18 +31,17 @@ namespace storm {
             void SparsePcaaAchievabilityQuery<SparseModelType, GeometryValueType>::initializeThresholdData() {
                 thresholds.reserve(this->objectives.size());
                 strictThresholds = storm::storage::BitVector(this->objectives.size(), false);
-                for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    auto const& obj = this->objectives[objIndex];
-                    STORM_LOG_ASSERT(obj.bound.is_initialized(), "Achievability query invoked but there is an objective without bound.");
-                    STORM_LOG_THROW(!obj.bound->threshold.containsVariables(), storm::exceptions::InvalidOperationException, "There is an objective whose bound contains undefined variables.");
-                    thresholds.push_back(storm::utility::convertNumber<GeometryValueType>(obj.bound->threshold.evaluateAsRational()));
-                    if (storm::solver::minimize(obj.optimizationDirection)) {
-                        STORM_LOG_ASSERT(!storm::logic::isLowerBound(obj.bound->comparisonType), "Minimizing objective should not specify an upper bound.");
+                for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
+                    auto const& formula = *this->objectives[objIndex].formula;
+                    STORM_LOG_ASSERT(formula.hasBound(), "Achievability query invoked but there is an objective without bound.");
+                    thresholds.push_back(formula.template getThresholdAs<GeometryValueType>());
+                    if (storm::solver::minimize(formula.getOptimalityType())) {
+                        STORM_LOG_ASSERT(!storm::logic::isLowerBound(formula.getBound().comparisonType), "Minimizing objective should not specify an upper bound.");
                         // Values for minimizing objectives will be negated in order to convert them to maximizing objectives.
                         // Hence, we also negate the threshold
                         thresholds.back() *= -storm::utility::one<GeometryValueType>();
                     }
-                    strictThresholds.set(objIndex, storm::logic::isStrict(obj.bound->comparisonType));
+                    strictThresholds.set(objIndex, storm::logic::isStrict(formula.getBound().comparisonType));
                 }
             }
             
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp
index c459ce7df..0405a6022 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuantitativeQuery.cpp
@@ -24,7 +24,7 @@ namespace storm {
                 STORM_LOG_ASSERT(preprocessorResult.queryType == SparseMultiObjectivePreprocessorReturnType<SparseModelType>::QueryType::Quantitative, "Invalid query Type");
                 
                 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    if (!this->objectives[objIndex].bound.is_initialized()) {
+                    if (!this->objectives[objIndex].formula->hasBound()) {
                         indexOfOptimizingObjective = objIndex;
                         break;
                     }
@@ -41,19 +41,19 @@ namespace storm {
                 thresholds.reserve(this->objectives.size());
                 strictThresholds = storm::storage::BitVector(this->objectives.size(), false);
                 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> thresholdConstraints;
-                thresholdConstraints.reserve(this->objectives.size()-1);
+                thresholdConstraints.reserve(this->objectives.size() - 1);
                 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    auto const& obj = this->objectives[objIndex];
-                    if (obj.bound) {
-                        STORM_LOG_THROW(!obj.bound->threshold.containsVariables(), storm::exceptions::InvalidOperationException, "There is an objective whose bound contains undefined variables.");
-                        thresholds.push_back(storm::utility::convertNumber<GeometryValueType>(obj.bound->threshold.evaluateAsRational()));
-                        if (storm::solver::minimize(obj.optimizationDirection)) {
-                            STORM_LOG_ASSERT(!storm::logic::isLowerBound(obj.bound->comparisonType), "Minimizing objective should not specify an upper bound.");
+                    auto const& formula = *this->objectives[objIndex].formula;
+
+                    if (formula.hasBound()) {
+                        thresholds.push_back(formula.template getThresholdAs<GeometryValueType>());
+                        if (storm::solver::minimize(formula.getOptimalityType())) {
+                            STORM_LOG_ASSERT(!storm::logic::isLowerBound(formula.getBound().comparisonType), "Minimizing objective should not specify an upper bound.");
                             // Values for minimizing objectives will be negated in order to convert them to maximizing objectives.
                             // Hence, we also negate the threshold
                             thresholds.back() *= -storm::utility::one<GeometryValueType>();
                         }
-                        strictThresholds.set(objIndex, storm::logic::isStrict(obj.bound->comparisonType));
+                        strictThresholds.set(objIndex, storm::logic::isStrict(formula.getBound().comparisonType));
                         WeightVector normalVector(this->objectives.size(), storm::utility::zero<GeometryValueType>());
                         normalVector[objIndex] = -storm::utility::one<GeometryValueType>();
                         thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back());
@@ -74,7 +74,7 @@ namespace storm {
                     
                     // transform the obtained result for the preprocessed model to a result w.r.t. the original model and return the checkresult
                     auto const& obj = this->objectives[indexOfOptimizingObjective];
-                    if (storm::solver::maximize(obj.optimizationDirection)) {
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
                         if (obj.considersComplementaryEvent) {
                             result = storm::utility::one<GeometryValueType>() - result;
                         }
@@ -95,7 +95,7 @@ namespace storm {
             
             template <class SparseModelType, typename GeometryValueType>
             bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::checkAchievability() {
-                 if (this->objectives.size()>1) {
+                 if (this->objectives.size() > 1) {
                     // We don't care for the optimizing objective at this point
                     this->diracWeightVectorsToBeChecked.set(indexOfOptimizingObjective, false);
                 
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp
index 9276613b1..a978dc53c 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaQuery.cpp
@@ -107,7 +107,7 @@ namespace storm {
                 step.upperBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getOverApproximationOfInitialStateResults());
                 // For the minimizing objectives, we need to scale the corresponding entries with -1 as we want to consider the downward closure
                 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
-                    if (storm::solver::minimize(this->objectives[objIndex].optimizationDirection)) {
+                    if (storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType())) {
                         step.lowerBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
                         step.upperBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
                     }
@@ -161,7 +161,7 @@ namespace storm {
                 result.reserve(point.size());
                 for(uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
                     auto const& obj = this->objectives[objIndex];
-                    if (storm::solver::maximize(obj.optimizationDirection)) {
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
                         if (obj.considersComplementaryEvent) {
                             result.push_back(storm::utility::one<GeometryValueType>() - point[objIndex]);
                         } else {
@@ -192,7 +192,7 @@ namespace storm {
                 transformationVector.reserve(numObjectives);
                 for(uint_fast64_t objIndex = 0; objIndex < numObjectives; ++objIndex) {
                     auto const& obj = this->objectives[objIndex];
-                    if (storm::solver::maximize(obj.optimizationDirection)) {
+                    if (storm::solver::maximize(obj.formula->getOptimalityType())) {
                         if (obj.considersComplementaryEvent) {
                             transformationMatrix[objIndex][objIndex] = -storm::utility::one<GeometryValueType>();
                             transformationVector.push_back(storm::utility::one<GeometryValueType>());
diff --git a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp
index be8727462..452254107 100644
--- a/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp
+++ b/src/storm/modelchecker/multiobjective/pcaa/SparsePcaaWeightVectorChecker.cpp
@@ -12,6 +12,7 @@
 #include "storm/utility/graph.h"
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
+#include "storm/logic/Formulas.h"
 
 #include "storm/exceptions/IllegalFunctionCallException.h"
 #include "storm/exceptions/UnexpectedException.h"
@@ -42,10 +43,11 @@ namespace storm {
                 
                 // set data for unbounded objectives
                 for(uint_fast64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) {
-                    auto const& obj = objectives[objIndex];
-                    if (!obj.upperTimeBound) {
+                    auto const& formula = *objectives[objIndex].formula;
+                    if (formula.getSubformula().isTotalRewardFormula()) {
                         objectivesWithNoUpperTimeBound.set(objIndex, true);
-                        actionsWithoutRewardInUnboundedPhase &= model.getRewardModel(*obj.rewardModelName).getChoicesWithZeroReward(model.getTransitionMatrix());
+                        STORM_LOG_ASSERT(formula.isRewardOperatorFormula() && formula.asRewardOperatorFormula().hasRewardModelName(), "Unexpected type of operator formula.");
+                        actionsWithoutRewardInUnboundedPhase &= model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName()).getChoicesWithZeroReward(model.getTransitionMatrix());
                     }
                 }
             }
@@ -59,26 +61,27 @@ namespace storm {
                 boost::optional<ValueType> weightedLowerResultBound = storm::utility::zero<ValueType>();
                 boost::optional<ValueType> weightedUpperResultBound = storm::utility::zero<ValueType>();
                 for (auto objIndex : objectivesWithNoUpperTimeBound) {
-                    if (storm::solver::minimize(objectives[objIndex].optimizationDirection)) {
-                        if (objectives[objIndex].lowerResultBound && weightedUpperResultBound) {
-                            weightedUpperResultBound.get() -= weightVector[objIndex] * objectives[objIndex].lowerResultBound.get();
+                    auto const& obj = objectives[objIndex];
+                    if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
+                        if (obj.lowerResultBound && weightedUpperResultBound) {
+                            weightedUpperResultBound.get() -= weightVector[objIndex] * obj.lowerResultBound.get();
                         } else {
                             weightedUpperResultBound = boost::none;
                         }
-                        if (objectives[objIndex].upperResultBound && weightedLowerResultBound) {
-                            weightedLowerResultBound.get() -= weightVector[objIndex] * objectives[objIndex].upperResultBound.get();
+                        if (obj.upperResultBound && weightedLowerResultBound) {
+                            weightedLowerResultBound.get() -= weightVector[objIndex] * obj.upperResultBound.get();
                         } else {
                             weightedLowerResultBound = boost::none;
                         }
                         storm::utility::vector::addScaledVector(weightedRewardVector, discreteActionRewards[objIndex], -weightVector[objIndex]);
                     } else {
-                        if (objectives[objIndex].lowerResultBound && weightedLowerResultBound) {
-                            weightedLowerResultBound.get() += weightVector[objIndex] * objectives[objIndex].lowerResultBound.get();
+                        if (obj.lowerResultBound && weightedLowerResultBound) {
+                            weightedLowerResultBound.get() += weightVector[objIndex] * obj.lowerResultBound.get();
                         } else {
                             weightedLowerResultBound = boost::none;
                         }
-                        if (objectives[objIndex].upperResultBound && weightedUpperResultBound) {
-                            weightedUpperResultBound.get() += weightVector[objIndex] * objectives[objIndex].upperResultBound.get();
+                        if (obj.upperResultBound && weightedUpperResultBound) {
+                            weightedUpperResultBound.get() += weightVector[objIndex] * obj.upperResultBound.get();
                         } else {
                             weightedUpperResultBound = boost::none;
                         }
@@ -90,8 +93,8 @@ namespace storm {
                 
                 unboundedIndividualPhase(weightVector);
                 // Only invoke boundedPhase if necessarry, i.e., if there is at least one objective with a time bound
-                for(auto const& obj : this->objectives) {
-                    if(obj.lowerTimeBound || obj.upperTimeBound) {
+                for (auto const& obj : this->objectives) {
+                    if (!obj.formula->getSubformula().isTotalRewardFormula()) {
                         boundedPhase(weightVector, weightedRewardVector);
                         break;
                     }
@@ -140,8 +143,8 @@ namespace storm {
             template <class SparseModelType>
             storm::storage::Scheduler<typename SparsePcaaWeightVectorChecker<SparseModelType>::ValueType> SparsePcaaWeightVectorChecker<SparseModelType>::computeScheduler() const {
                 STORM_LOG_THROW(this->checkHasBeenCalled, storm::exceptions::IllegalFunctionCallException, "Tried to retrieve results but check(..) has not been called before.");
-                for(auto const& obj : this->objectives) {
-                    STORM_LOG_THROW(!obj.lowerTimeBound && !obj.upperTimeBound, storm::exceptions::NotImplementedException, "Scheduler retrival is not implemented for timeBounded objectives.");
+                for (auto const& obj : this->objectives) {
+                    STORM_LOG_THROW(obj.formula->getSubformula().isTotalRewardFormula(), storm::exceptions::NotImplementedException, "Scheduler retrival is only implemented for objectives without time-bound.");
                 }
                 
                 storm::storage::Scheduler<ValueType> result(this->optimalChoices.size());
@@ -202,7 +205,7 @@ namespace storm {
                 if (objectivesWithNoUpperTimeBound.getNumberOfSetBits() == 1 && storm::utility::isOne(weightVector[*objectivesWithNoUpperTimeBound.begin()])) {
                    uint_fast64_t objIndex = *objectivesWithNoUpperTimeBound.begin();
                    objectiveResults[objIndex] = weightedResult;
-                   if (storm::solver::minimize(objectives[objIndex].optimizationDirection)) {
+                   if (storm::solver::minimize(objectives[objIndex].formula->getOptimalityType())) {
                        storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], -storm::utility::one<ValueType>());
                    }
                     for (uint_fast64_t objIndex2 = 0; objIndex2 < objectives.size(); ++objIndex2) {
@@ -235,7 +238,7 @@ namespace storm {
                            if (!storm::utility::isZero(weightVector[objIndex])) {
                                objectiveResults[objIndex] = weightedSumOfUncheckedObjectives;
                                ValueType scalingFactor = storm::utility::one<ValueType>() / sumOfWeightsOfUncheckedObjectives;
-                               if (storm::solver::minimize(obj.optimizationDirection)) {
+                               if (storm::solver::minimize(obj.formula->getOptimalityType())) {
                                    scalingFactor *= -storm::utility::one<ValueType>();
                                }
                                storm::utility::vector::scaleVectorInPlace(objectiveResults[objIndex], scalingFactor);
diff --git a/src/storm/models/sparse/NondeterministicModel.cpp b/src/storm/models/sparse/NondeterministicModel.cpp
index e823d455a..3382bfd95 100644
--- a/src/storm/models/sparse/NondeterministicModel.cpp
+++ b/src/storm/models/sparse/NondeterministicModel.cpp
@@ -2,6 +2,9 @@
 
 #include "storm/models/sparse/StandardRewardModel.h"
 #include "storm/models/sparse/MarkovAutomaton.h"
+#include "storm/storage/Scheduler.h"
+#include "storm/storage/memorystructure/MemoryStructureBuilder.h"
+#include "storm/storage/memorystructure/SparseModelMemoryProduct.h"
 
 #include "storm/adapters/RationalFunctionAdapter.h"
 
@@ -46,6 +49,26 @@ namespace storm {
                 }
             }
             
+            template<typename ValueType, typename RewardModelType>
+            std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> NondeterministicModel<ValueType, RewardModelType>::applyScheduler(storm::storage::Scheduler<ValueType> const& scheduler, bool dropUnreachableStates) {
+                if (scheduler.isMemorylessScheduler()) {
+                    auto memStruct = storm::storage::MemoryStructureBuilder<ValueType, RewardModelType>::buildTrivialMemoryStructure(*this);
+                    auto memoryProduct = memStruct.product(*this);
+                    if (!dropUnreachableStates) {
+                        memoryProduct.setBuildFullProduct();
+                    }
+                    return memoryProduct.build(scheduler);
+                } else {
+                    boost::optional<storm::storage::MemoryStructure> const& memStruct = scheduler.getMemoryStructure();
+                    STORM_LOG_ASSERT(memStruct, "Memoryless scheduler without memory structure.");
+                    auto memoryProduct = memStruct->product(*this);
+                    if (!dropUnreachableStates) {
+                        memoryProduct.setBuildFullProduct();
+                    }
+                    return memoryProduct.build(scheduler);
+                }
+            }
+            
             template<typename ValueType, typename RewardModelType>
             void NondeterministicModel<ValueType, RewardModelType>::printModelInformationToStream(std::ostream& out) const {
                 this->printModelInformationHeaderToStream(out);
diff --git a/src/storm/models/sparse/NondeterministicModel.h b/src/storm/models/sparse/NondeterministicModel.h
index 1f2bd157a..ab0a21a04 100644
--- a/src/storm/models/sparse/NondeterministicModel.h
+++ b/src/storm/models/sparse/NondeterministicModel.h
@@ -5,6 +5,13 @@
 #include "storm/utility/OsDetection.h"
 
 namespace storm {
+    
+    // Forward declare Scheduler class.
+    namespace storage {
+        template <typename ValueType>
+        class Scheduler;
+    }
+    
     namespace models {
         namespace sparse {
             
@@ -48,6 +55,13 @@ namespace storm {
                 
                 virtual void reduceToStateBasedRewards() override;
                 
+                /*!
+                 * Applies the given scheduler to this model.
+                 * @param scheduler the considered scheduler.
+                 * @param dropUnreachableStates if set, the resulting model only considers the states that are reachable from an initial state
+                 */
+                std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> applyScheduler(storm::storage::Scheduler<ValueType> const& scheduler, bool dropUnreachableStates = true);
+                
                 virtual void printModelInformationToStream(std::ostream& out) const override;
                 
                 virtual void writeDotToStream(std::ostream& outStream, bool includeLabeling = true, storm::storage::BitVector const* subsystem = nullptr, std::vector<ValueType> const* firstValue = nullptr, std::vector<ValueType> const* secondValue = nullptr, std::vector<uint_fast64_t> const* stateColoring = nullptr, std::vector<std::string> const* colors = nullptr, std::vector<uint_fast64_t>* scheduler = nullptr, bool finalizeOutput = true) const override;
diff --git a/src/storm/storage/Distribution.cpp b/src/storm/storage/Distribution.cpp
index 0c2ac73db..9d0d89d19 100644
--- a/src/storm/storage/Distribution.cpp
+++ b/src/storm/storage/Distribution.cpp
@@ -156,6 +156,18 @@ namespace storm {
             return false;
         }
         
+        template<typename ValueType, typename StateType>
+        ValueType Distribution<ValueType, StateType>::getProbability(StateType const& state) const {
+            auto it = this->distribution.find(state);
+            if (it == this->distribution.end()) {
+                return storm::utility::zero<ValueType>();
+            } else {
+                return it->second;
+            }
+        }
+        
+
+        
         template class Distribution<double>;
         template std::ostream& operator<<(std::ostream& out, Distribution<double> const& distribution);
         template class Distribution<double, uint_fast64_t>;
diff --git a/src/storm/storage/Distribution.h b/src/storm/storage/Distribution.h
index ff4b0689b..328e87f98 100644
--- a/src/storm/storage/Distribution.h
+++ b/src/storm/storage/Distribution.h
@@ -73,7 +73,7 @@ namespace storm {
              * entry is removed.
              */
             void shiftProbability(StateType const& fromState, StateType const& toState, ValueType const& probability, storm::utility::ConstantsComparator<ValueType> const& comparator = storm::utility::ConstantsComparator<ValueType>());
-                        
+
             /*!
              * Retrieves an iterator to the elements in this distribution.
              *
@@ -132,6 +132,14 @@ namespace storm {
             
             bool less(Distribution<ValueType, StateType> const& other, storm::utility::ConstantsComparator<ValueType> const& comparator) const;
             
+            
+            /*!
+             * Returns the probability of the given state
+             * @param state The state for which the probability is returned.
+             * @return The probability of the given state.
+             */
+            ValueType getProbability(StateType const& state) const;
+            
         private:
             // A list of states and the probabilities that are assigned to them.
             container_type distribution;
diff --git a/src/storm/storage/Scheduler.cpp b/src/storm/storage/Scheduler.cpp
index a1dc1ecbd..7ddc7e457 100644
--- a/src/storm/storage/Scheduler.cpp
+++ b/src/storm/storage/Scheduler.cpp
@@ -86,6 +86,11 @@ namespace storm {
             return memoryStructure ? memoryStructure->getNumberOfStates() : 1;
         }
 
+        template <typename ValueType>
+        boost::optional<storm::storage::MemoryStructure> const& Scheduler<ValueType>::getMemoryStructure() const {
+            return memoryStructure;
+        }
+
         template <typename ValueType>
         void Scheduler<ValueType>::printToStream(std::ostream& out, std::shared_ptr<storm::models::sparse::Model<ValueType>> model, bool skipUniqueChoices) const {
             STORM_LOG_THROW(model == nullptr || model->getNumberOfStates() == schedulerChoices.front().size(), storm::exceptions::InvalidOperationException, "The given model is not compatible with this scheduler.");
diff --git a/src/storm/storage/Scheduler.h b/src/storm/storage/Scheduler.h
index c1c3a5f88..b51e308b9 100644
--- a/src/storm/storage/Scheduler.h
+++ b/src/storm/storage/Scheduler.h
@@ -70,6 +70,11 @@ namespace storm {
              * Retrieves the number of memory states this scheduler considers.
              */
             uint_fast64_t getNumberOfMemoryStates() const;
+            
+            /*!
+             * Retrieves the memory structure associated with this scheduler
+             */
+            boost::optional<storm::storage::MemoryStructure> const& getMemoryStructure() const;
 
             /*!
              * Returns a copy of this scheduler with the new value type
diff --git a/src/storm/storage/memorystructure/MemoryStructure.cpp b/src/storm/storage/memorystructure/MemoryStructure.cpp
index 3d9348c44..499bc4db6 100644
--- a/src/storm/storage/memorystructure/MemoryStructure.cpp
+++ b/src/storm/storage/memorystructure/MemoryStructure.cpp
@@ -106,9 +106,9 @@ namespace storm {
             return MemoryStructure(std::move(resultTransitions), std::move(resultLabeling), std::move(resultInitialMemoryStates));
         }
             
-        template <typename ValueType>
-        SparseModelMemoryProduct<ValueType> MemoryStructure::product(storm::models::sparse::Model<ValueType> const& sparseModel) const {
-            return SparseModelMemoryProduct<ValueType>(sparseModel, *this);
+        template <typename ValueType, typename RewardModelType>
+        SparseModelMemoryProduct<ValueType, RewardModelType> MemoryStructure::product(storm::models::sparse::Model<ValueType, RewardModelType> const& sparseModel) const {
+            return SparseModelMemoryProduct<ValueType, RewardModelType>(sparseModel, *this);
         }
         
         std::string MemoryStructure::toString() const {
@@ -143,7 +143,9 @@ namespace storm {
         }
         
         template SparseModelMemoryProduct<double> MemoryStructure::product(storm::models::sparse::Model<double> const& sparseModel) const;
+        template SparseModelMemoryProduct<double, storm::models::sparse::StandardRewardModel<storm::Interval>> MemoryStructure::product(storm::models::sparse::Model<double, storm::models::sparse::StandardRewardModel<storm::Interval>> const& sparseModel) const;
         template SparseModelMemoryProduct<storm::RationalNumber> MemoryStructure::product(storm::models::sparse::Model<storm::RationalNumber> const& sparseModel) const;
+        template SparseModelMemoryProduct<storm::RationalFunction> MemoryStructure::product(storm::models::sparse::Model<storm::RationalFunction> const& sparseModel) const;
     }
 }
 
diff --git a/src/storm/storage/memorystructure/MemoryStructure.h b/src/storm/storage/memorystructure/MemoryStructure.h
index 9df320fd3..64bd7aaa9 100644
--- a/src/storm/storage/memorystructure/MemoryStructure.h
+++ b/src/storm/storage/memorystructure/MemoryStructure.h
@@ -6,11 +6,12 @@
 #include "storm/logic/Formula.h"
 #include "storm/models/sparse/StateLabeling.h"
 #include "storm/models/sparse/Model.h"
+#include "storm/models/sparse/StandardRewardModel.h"
 
 namespace storm {
     namespace storage {
         
-        template <typename ValueType>
+        template <typename ValueType, typename RewardModelType>
         class SparseModelMemoryProduct;
         
         /*!
@@ -53,8 +54,8 @@ namespace storm {
              * Builds the product of this memory structure and the given sparse model.
              * An exception is thrown if the state labelings of this memory structure and the given model are not disjoint.
              */
-            template <typename ValueType>
-            SparseModelMemoryProduct<ValueType> product(storm::models::sparse::Model<ValueType> const& sparseModel) const;
+            template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+            SparseModelMemoryProduct<ValueType, RewardModelType> product(storm::models::sparse::Model<ValueType, RewardModelType> const& sparseModel) const;
             
             std::string toString() const;
 
diff --git a/src/storm/storage/memorystructure/MemoryStructureBuilder.cpp b/src/storm/storage/memorystructure/MemoryStructureBuilder.cpp
index 9ee23eaa6..148cbf2ca 100644
--- a/src/storm/storage/memorystructure/MemoryStructureBuilder.cpp
+++ b/src/storm/storage/memorystructure/MemoryStructureBuilder.cpp
@@ -83,7 +83,15 @@ namespace storm {
             return MemoryStructure(std::move(transitions), std::move(stateLabeling), std::move(initialMemoryStates));
         }
         
+        template <typename ValueType, typename RewardModelType>
+        MemoryStructure MemoryStructureBuilder<ValueType, RewardModelType>::buildTrivialMemoryStructure(storm::models::sparse::Model<ValueType, RewardModelType> const& model) {
+            MemoryStructureBuilder<ValueType, RewardModelType> memoryBuilder(1, model);
+            memoryBuilder.setTransition(0,0, storm::storage::BitVector(model.getNumberOfStates(), true));
+            return memoryBuilder.build();
+        }
+
         template class MemoryStructureBuilder<double>;
+        template class MemoryStructureBuilder<double, storm::models::sparse::StandardRewardModel<storm::Interval>>;
         template class MemoryStructureBuilder<storm::RationalNumber>;
         template class MemoryStructureBuilder<storm::RationalFunction>;
         
diff --git a/src/storm/storage/memorystructure/MemoryStructureBuilder.h b/src/storm/storage/memorystructure/MemoryStructureBuilder.h
index 9fbe65c78..1cb51f384 100644
--- a/src/storm/storage/memorystructure/MemoryStructureBuilder.h
+++ b/src/storm/storage/memorystructure/MemoryStructureBuilder.h
@@ -57,6 +57,11 @@ namespace storm {
              */
             MemoryStructure build();
             
+            /*!
+             * Builds a trivial memory structure for the given model (consisting of a single memory state)
+             */
+            static MemoryStructure buildTrivialMemoryStructure(storm::models::sparse::Model<ValueType, RewardModelType> const& model);
+            
 
         private:
             storm::models::sparse::Model<ValueType, RewardModelType> const& model;
diff --git a/src/storm/storage/memorystructure/SparseModelMemoryProduct.cpp b/src/storm/storage/memorystructure/SparseModelMemoryProduct.cpp
index 56b3a6aa0..58f773cfe 100644
--- a/src/storm/storage/memorystructure/SparseModelMemoryProduct.cpp
+++ b/src/storm/storage/memorystructure/SparseModelMemoryProduct.cpp
@@ -17,17 +17,29 @@
 namespace storm {
     namespace storage {
 
-        template <typename ValueType>
-        SparseModelMemoryProduct<ValueType>::SparseModelMemoryProduct(storm::models::sparse::Model<ValueType> const& sparseModel, storm::storage::MemoryStructure const& memoryStructure) : model(sparseModel), memory(memoryStructure) {
-            // intentionally left empty
+        template <typename ValueType, typename RewardModelType>
+        SparseModelMemoryProduct<ValueType, RewardModelType>::SparseModelMemoryProduct(storm::models::sparse::Model<ValueType, RewardModelType> const& sparseModel, storm::storage::MemoryStructure const& memoryStructure) : model(sparseModel), memory(memoryStructure) {
+            reachableStates = storm::storage::BitVector(model.getNumberOfStates() * memory.getNumberOfStates(), false);
         }
+        
+        template <typename ValueType, typename RewardModelType>
+        void SparseModelMemoryProduct<ValueType, RewardModelType>::addReachableState(uint64_t const& modelState, uint64_t const& memoryState) {
+            reachableStates.set(modelState * memory.getNumberOfStates() + memoryState, true);
+        }
+
+       template <typename ValueType, typename RewardModelType>
+       void SparseModelMemoryProduct<ValueType, RewardModelType>::setBuildFullProduct() {
+            reachableStates.clear();
+            reachableStates.complement();
+       }
+        
+        template <typename ValueType, typename RewardModelType>
+        std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> SparseModelMemoryProduct<ValueType, RewardModelType>::build(boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) {
             
-        template <typename ValueType>
-        std::shared_ptr<storm::models::sparse::Model<ValueType>> SparseModelMemoryProduct<ValueType>::build() {
-            uint_fast64_t modelStateCount = model.getNumberOfStates();
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
+            uint64_t modelStateCount = model.getNumberOfStates();
+            uint64_t memoryStateCount = memory.getNumberOfStates();
             
-            std::vector<uint_fast64_t> memorySuccessors = computeMemorySuccessors();
+            std::vector<uint64_t> memorySuccessors = computeMemorySuccessors();
             
             // Get the initial states and reachable states. A stateIndex s corresponds to the model state (s / memoryStateCount) and memory state (s % memoryStateCount)
             storm::storage::BitVector initialStates(modelStateCount * memoryStateCount, false);
@@ -38,44 +50,49 @@ namespace storm {
             }
             STORM_LOG_ASSERT(memoryInitIt == memory.getInitialMemoryStates().end(), "Unexpected number of initial states.");
             
-            storm::storage::BitVector reachableStates = computeReachableStates(memorySuccessors, initialStates);
+            computeReachableStates(memorySuccessors, initialStates, scheduler);
             
             // Compute the mapping to the states of the result
-            uint_fast64_t reachableStateCount = 0;
-            toResultStateMapping = std::vector<uint_fast64_t> (model.getNumberOfStates() * memory.getNumberOfStates(), std::numeric_limits<uint_fast64_t>::max());
+            uint64_t reachableStateCount = 0;
+            toResultStateMapping = std::vector<uint64_t> (model.getNumberOfStates() * memory.getNumberOfStates(), std::numeric_limits<uint64_t>::max());
             for (auto const& reachableState : reachableStates) {
                 toResultStateMapping[reachableState] = reachableStateCount;
                 ++reachableStateCount;
             }
                 
             // Build the model components
-            storm::storage::SparseMatrix<ValueType> transitionMatrix = model.getTransitionMatrix().hasTrivialRowGrouping() ?
-                                                                       buildDeterministicTransitionMatrix(reachableStates, memorySuccessors) :
-                                                                       buildNondeterministicTransitionMatrix(reachableStates, memorySuccessors);
+            storm::storage::SparseMatrix<ValueType> transitionMatrix;
+            if (scheduler) {
+                transitionMatrix = buildTransitionMatrixForScheduler(memorySuccessors, scheduler.get());
+            } else if (model.getTransitionMatrix().hasTrivialRowGrouping()) {
+                transitionMatrix = buildDeterministicTransitionMatrix(memorySuccessors);
+            } else {
+                transitionMatrix = buildNondeterministicTransitionMatrix(memorySuccessors);
+            }
             storm::models::sparse::StateLabeling labeling = buildStateLabeling(transitionMatrix);
-            std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>> rewardModels = buildRewardModels(transitionMatrix, memorySuccessors);
+            std::unordered_map<std::string, RewardModelType> rewardModels = buildRewardModels(transitionMatrix, memorySuccessors, scheduler);
 
             // Add the label for the initial states. We need to translate the state indices w.r.t. the set of reachable states.
             labeling.addLabel("init", initialStates % reachableStates);
             
             
-            return buildResult(std::move(transitionMatrix), std::move(labeling), std::move(rewardModels));
+            return buildResult(std::move(transitionMatrix), std::move(labeling), std::move(rewardModels), scheduler);
 
         }
             
-        template <typename ValueType>
-        uint_fast64_t const& SparseModelMemoryProduct<ValueType>::getResultState(uint_fast64_t const& modelState, uint_fast64_t const& memoryState) const {
+        template <typename ValueType, typename RewardModelType>
+        uint64_t const& SparseModelMemoryProduct<ValueType, RewardModelType>::getResultState(uint64_t const& modelState, uint64_t const& memoryState) const {
                 return toResultStateMapping[modelState * memory.getNumberOfStates() + memoryState];
         }
             
-        template <typename ValueType>
-        std::vector<uint_fast64_t> SparseModelMemoryProduct<ValueType>::computeMemorySuccessors() const {
-            uint_fast64_t modelTransitionCount = model.getTransitionMatrix().getEntryCount();
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
-            std::vector<uint_fast64_t> result(modelTransitionCount * memoryStateCount, std::numeric_limits<uint_fast64_t>::max());
+        template <typename ValueType, typename RewardModelType>
+        std::vector<uint64_t> SparseModelMemoryProduct<ValueType, RewardModelType>::computeMemorySuccessors() const {
+            uint64_t modelTransitionCount = model.getTransitionMatrix().getEntryCount();
+            uint64_t memoryStateCount = memory.getNumberOfStates();
+            std::vector<uint64_t> result(modelTransitionCount * memoryStateCount, std::numeric_limits<uint64_t>::max());
             
-            for (uint_fast64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
-                for (uint_fast64_t transitionGoal = 0; transitionGoal < memoryStateCount; ++transitionGoal) {
+            for (uint64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
+                for (uint64_t transitionGoal = 0; transitionGoal < memoryStateCount; ++transitionGoal) {
                     auto const& memoryTransition = memory.getTransitionMatrix()[memoryState][transitionGoal];
                     if (memoryTransition) {
                         for (auto const& modelTransitionIndex : memoryTransition.get()) {
@@ -87,54 +104,76 @@ namespace storm {
             return result;
         }
             
-        template <typename ValueType>
-        storm::storage::BitVector SparseModelMemoryProduct<ValueType>::computeReachableStates(std::vector<uint_fast64_t> const& memorySuccessors, storm::storage::BitVector const& initialStates) const {
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
+        template <typename ValueType, typename RewardModelType>
+        void SparseModelMemoryProduct<ValueType, RewardModelType>::computeReachableStates(std::vector<uint64_t> const& memorySuccessors, storm::storage::BitVector const& initialStates,                     boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) {
+            uint64_t memoryStateCount = memory.getNumberOfStates();
             // Explore the reachable states via DFS.
             // A state s on the stack corresponds to the model state (s / memoryStateCount) and memory state (s % memoryStateCount)
-            storm::storage::BitVector reachableStates = initialStates;
-            std::vector<uint_fast64_t> stack(reachableStates.begin(), reachableStates.end());
-            while (!stack.empty()) {
-                uint_fast64_t stateIndex = stack.back();
-                stack.pop_back();
-                uint_fast64_t modelState = stateIndex / memoryStateCount;
-                uint_fast64_t memoryState = stateIndex % memoryStateCount;
-                
-                auto const& rowGroup = model.getTransitionMatrix().getRowGroup(modelState);
-                for (auto modelTransitionIt = rowGroup.begin(); modelTransitionIt != rowGroup.end(); ++modelTransitionIt) {
-                    if (!storm::utility::isZero(modelTransitionIt->getValue())) {
-                        uint_fast64_t successorModelState = modelTransitionIt->getColumn();
-                        uint_fast64_t modelTransitionId = modelTransitionIt - model.getTransitionMatrix().begin();
-                        uint_fast64_t successorMemoryState = memorySuccessors[modelTransitionId * memoryStateCount + memoryState];
-                        uint_fast64_t successorStateIndex = successorModelState * memoryStateCount + successorMemoryState;
-                        if (!reachableStates.get(successorStateIndex)) {
-                            reachableStates.set(successorStateIndex, true);
-                            stack.push_back(successorStateIndex);
+            reachableStates |= initialStates;
+            if (!reachableStates.full()) {
+                std::vector<uint64_t> stack(reachableStates.begin(), reachableStates.end());
+                while (!stack.empty()) {
+                    uint64_t stateIndex = stack.back();
+                    stack.pop_back();
+                    uint64_t modelState = stateIndex / memoryStateCount;
+                    uint64_t memoryState = stateIndex % memoryStateCount;
+                    
+                    if (scheduler) {
+                        auto choices = scheduler->getChoice(modelState, memoryState).getChoiceAsDistribution();
+                        uint64_t groupStart = model.getTransitionMatrix().getRowGroupIndices()[modelState];
+                        for (auto const& choice : choices) {
+                            STORM_LOG_ASSERT(groupStart + choice.first < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1], "Invalid choice " << choice.first << " at model state " << modelState << ".");
+                            auto const& row = model.getTransitionMatrix().getRow(groupStart + choice.first);
+                            for (auto modelTransitionIt = row.begin(); modelTransitionIt != row.end(); ++modelTransitionIt) {
+                                if (!storm::utility::isZero(modelTransitionIt->getValue())) {
+                                    uint64_t successorModelState = modelTransitionIt->getColumn();
+                                    uint64_t modelTransitionId = modelTransitionIt - model.getTransitionMatrix().begin();
+                                    uint64_t successorMemoryState = memorySuccessors[modelTransitionId * memoryStateCount + memoryState];
+                                    uint64_t successorStateIndex = successorModelState * memoryStateCount + successorMemoryState;
+                                    if (!reachableStates.get(successorStateIndex)) {
+                                        reachableStates.set(successorStateIndex, true);
+                                        stack.push_back(successorStateIndex);
+                                    }
+                                }
+                            }
+                        }
+                    } else {
+                        auto const& rowGroup = model.getTransitionMatrix().getRowGroup(modelState);
+                        for (auto modelTransitionIt = rowGroup.begin(); modelTransitionIt != rowGroup.end(); ++modelTransitionIt) {
+                            if (!storm::utility::isZero(modelTransitionIt->getValue())) {
+                                uint64_t successorModelState = modelTransitionIt->getColumn();
+                                uint64_t modelTransitionId = modelTransitionIt - model.getTransitionMatrix().begin();
+                                uint64_t successorMemoryState = memorySuccessors[modelTransitionId * memoryStateCount + memoryState];
+                                uint64_t successorStateIndex = successorModelState * memoryStateCount + successorMemoryState;
+                                if (!reachableStates.get(successorStateIndex)) {
+                                    reachableStates.set(successorStateIndex, true);
+                                    stack.push_back(successorStateIndex);
+                                }
+                            }
                         }
                     }
                 }
             }
-            return reachableStates;
         }
-            
-        template <typename ValueType>
-        storm::storage::SparseMatrix<ValueType> SparseModelMemoryProduct<ValueType>::buildDeterministicTransitionMatrix(storm::storage::BitVector const& reachableStates, std::vector<uint_fast64_t> const& memorySuccessors) const {
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
-            uint_fast64_t numResStates = reachableStates.getNumberOfSetBits();
-            uint_fast64_t numResTransitions = 0;
+        
+        template <typename ValueType, typename RewardModelType>
+        storm::storage::SparseMatrix<ValueType> SparseModelMemoryProduct<ValueType, RewardModelType>::buildDeterministicTransitionMatrix(std::vector<uint64_t> const& memorySuccessors) const {
+            uint64_t memoryStateCount = memory.getNumberOfStates();
+            uint64_t numResStates = reachableStates.getNumberOfSetBits();
+            uint64_t numResTransitions = 0;
             for (auto const& stateIndex : reachableStates) {
                 numResTransitions += model.getTransitionMatrix().getRow(stateIndex / memoryStateCount).getNumberOfEntries();
             }
             
             storm::storage::SparseMatrixBuilder<ValueType> builder(numResStates, numResStates, numResTransitions, true);
-            uint_fast64_t currentRow = 0;
+            uint64_t currentRow = 0;
             for (auto const& stateIndex : reachableStates) {
-                uint_fast64_t modelState = stateIndex / memoryStateCount;
-                uint_fast64_t memoryState = stateIndex % memoryStateCount;
+                uint64_t modelState = stateIndex / memoryStateCount;
+                uint64_t memoryState = stateIndex % memoryStateCount;
                 auto const& modelRow = model.getTransitionMatrix().getRow(modelState);
                 for (auto entryIt = modelRow.begin(); entryIt != modelRow.end(); ++entryIt) {
-                    uint_fast64_t transitionId = entryIt - model.getTransitionMatrix().begin();
-                    uint_fast64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                    uint64_t transitionId = entryIt - model.getTransitionMatrix().begin();
+                    uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
                     builder.addNextValue(currentRow, getResultState(entryIt->getColumn(), successorMemoryState), entryIt->getValue());
                 }
                 ++currentRow;
@@ -143,31 +182,31 @@ namespace storm {
             return builder.build();
         }
         
-        template <typename ValueType>
-        storm::storage::SparseMatrix<ValueType> SparseModelMemoryProduct<ValueType>::buildNondeterministicTransitionMatrix(storm::storage::BitVector const& reachableStates, std::vector<uint_fast64_t> const& memorySuccessors) const {
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
-            uint_fast64_t numResStates = reachableStates.getNumberOfSetBits();
-            uint_fast64_t numResChoices = 0;
-            uint_fast64_t numResTransitions = 0;
+        template <typename ValueType, typename RewardModelType>
+        storm::storage::SparseMatrix<ValueType> SparseModelMemoryProduct<ValueType, RewardModelType>::buildNondeterministicTransitionMatrix(std::vector<uint64_t> const& memorySuccessors) const {
+            uint64_t memoryStateCount = memory.getNumberOfStates();
+            uint64_t numResStates = reachableStates.getNumberOfSetBits();
+            uint64_t numResChoices = 0;
+            uint64_t numResTransitions = 0;
             for (auto const& stateIndex : reachableStates) {
-                uint_fast64_t modelState = stateIndex / memoryStateCount;
-                for (uint_fast64_t modelRow = model.getTransitionMatrix().getRowGroupIndices()[modelState]; modelRow < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]; ++modelRow) {
+                uint64_t modelState = stateIndex / memoryStateCount;
+                for (uint64_t modelRow = model.getTransitionMatrix().getRowGroupIndices()[modelState]; modelRow < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]; ++modelRow) {
                     ++numResChoices;
                     numResTransitions += model.getTransitionMatrix().getRow(modelRow).getNumberOfEntries();
                 }
             }
             
             storm::storage::SparseMatrixBuilder<ValueType> builder(numResChoices, numResStates, numResTransitions, true, true, numResStates);
-            uint_fast64_t currentRow = 0;
+            uint64_t currentRow = 0;
             for (auto const& stateIndex : reachableStates) {
-                uint_fast64_t modelState = stateIndex / memoryStateCount;
-                uint_fast64_t memoryState = stateIndex % memoryStateCount;
+                uint64_t modelState = stateIndex / memoryStateCount;
+                uint64_t memoryState = stateIndex % memoryStateCount;
                 builder.newRowGroup(currentRow);
-                for (uint_fast64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState]; modelRowIndex < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]; ++modelRowIndex) {
+                for (uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState]; modelRowIndex < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]; ++modelRowIndex) {
                     auto const& modelRow = model.getTransitionMatrix().getRow(modelRowIndex);
                     for (auto entryIt = modelRow.begin(); entryIt != modelRow.end(); ++entryIt) {
-                        uint_fast64_t transitionId = entryIt - model.getTransitionMatrix().begin();
-                        uint_fast64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                        uint64_t transitionId = entryIt - model.getTransitionMatrix().begin();
+                        uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
                         builder.addNextValue(currentRow, getResultState(entryIt->getColumn(), successorMemoryState), entryIt->getValue());
                     }
                     ++currentRow;
@@ -177,20 +216,115 @@ namespace storm {
             return builder.build();
         }
         
-        template <typename ValueType>
-        storm::models::sparse::StateLabeling SparseModelMemoryProduct<ValueType>::buildStateLabeling(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix) const {
-            uint_fast64_t modelStateCount = model.getNumberOfStates();
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
+        template <typename ValueType, typename RewardModelType>
+        storm::storage::SparseMatrix<ValueType> SparseModelMemoryProduct<ValueType, RewardModelType>::buildTransitionMatrixForScheduler(std::vector<uint64_t> const& memorySuccessors, storm::storage::Scheduler<ValueType> const& scheduler) const {
+            uint64_t memoryStateCount = memory.getNumberOfStates();
+            uint64_t numResStates = reachableStates.getNumberOfSetBits();
+            uint64_t numResChoices = 0;
+            uint64_t numResTransitions = 0;
+            bool hasTrivialNondeterminism = true;
+            for (auto const& stateIndex : reachableStates) {
+                uint64_t modelState = stateIndex / memoryStateCount;
+                uint64_t memoryState = stateIndex % memoryStateCount;
+                storm::storage::SchedulerChoice<ValueType> choice = scheduler.getChoice(modelState, memoryState);
+                if (choice.isDefined()) {
+                    ++numResChoices;
+                    if (choice.isDeterministic()) {
+                        uint64_t modelRow = model.getTransitionMatrix().getRowGroupIndices()[modelState] + choice.getDeterministicChoice();
+                        numResTransitions += model.getTransitionMatrix().getRow(modelRow).getNumberOfEntries();
+                    } else {
+                        std::set<uint64_t> successors;
+                        for (auto const& choiceIndex : choice.getChoiceAsDistribution()) {
+                            if (!storm::utility::isZero(choiceIndex.second)) {
+                                uint64_t modelRow = model.getTransitionMatrix().getRowGroupIndices()[modelState] + choiceIndex.first;
+                                for (auto const& entry : model.getTransitionMatrix().getRow(modelRow)) {
+                                    successors.insert(entry.getColumn());
+                                }
+                            }
+                        }
+                        numResTransitions += successors.size();
+                    }
+                } else {
+                    uint64_t modelRow = model.getTransitionMatrix().getRowGroupIndices()[modelState];
+                    uint64_t groupEnd = model.getTransitionMatrix().getRowGroupIndices()[modelState + 1];
+                    hasTrivialNondeterminism = hasTrivialNondeterminism && (groupEnd == modelRow + 1);
+                    for (; groupEnd; ++modelRow) {
+                        ++numResChoices;
+                        numResTransitions += model.getTransitionMatrix().getRow(modelRow).getNumberOfEntries();
+                    }
+                }
+            }
+            
+            storm::storage::SparseMatrixBuilder<ValueType> builder(numResChoices, numResStates, numResTransitions, true, !hasTrivialNondeterminism, hasTrivialNondeterminism ? 0 : numResStates);
+            uint64_t currentRow = 0;
+            for (auto const& stateIndex : reachableStates) {
+                uint64_t modelState = stateIndex / memoryStateCount;
+                uint64_t memoryState = stateIndex % memoryStateCount;
+                if (!hasTrivialNondeterminism) {
+                    builder.newRowGroup(currentRow);
+                }
+                storm::storage::SchedulerChoice<ValueType> choice = scheduler.getChoice(modelState, memoryState);
+                if (choice.isDefined()) {
+                    if (choice.isDeterministic()) {
+                        uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState] + choice.getDeterministicChoice();
+                        auto const& modelRow = model.getTransitionMatrix().getRow(modelRowIndex);
+                        for (auto entryIt = modelRow.begin(); entryIt != modelRow.end(); ++entryIt) {
+                            uint64_t transitionId = entryIt - model.getTransitionMatrix().begin();
+                            uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                            builder.addNextValue(currentRow, getResultState(entryIt->getColumn(), successorMemoryState), entryIt->getValue());
+                        }
+                    } else {
+                        std::map<uint64_t, ValueType> transitions;
+                        for (auto const& choiceIndex : choice.getChoiceAsDistribution()) {
+                            if (!storm::utility::isZero(choiceIndex.second)) {
+                                uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState] + choiceIndex.first;
+                                auto const& modelRow = model.getTransitionMatrix().getRow(modelRowIndex);
+                                for (auto entryIt = modelRow.begin(); entryIt != modelRow.end(); ++entryIt) {
+                                    uint64_t transitionId = entryIt - model.getTransitionMatrix().begin();
+                                    uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                                    ValueType transitionValue = choiceIndex.second * entryIt->getValue();
+                                    auto insertionRes = transitions.insert(std::make_pair(getResultState(entryIt->getColumn(), successorMemoryState), transitionValue));
+                                    if (!insertionRes.second) {
+                                        insertionRes.first->second += transitionValue;
+                                    }
+                                }
+                            }
+                        }
+                        for (auto const& transition : transitions) {
+                            builder.addNextValue(currentRow, transition.first, transition.second);
+                        }
+                    }
+                    ++currentRow;
+                } else {
+                    for (uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState]; modelRowIndex < model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]; ++modelRowIndex) {
+                        auto const& modelRow = model.getTransitionMatrix().getRow(modelRowIndex);
+                        for (auto entryIt = modelRow.begin(); entryIt != modelRow.end(); ++entryIt) {
+                            uint64_t transitionId = entryIt - model.getTransitionMatrix().begin();
+                            uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                            builder.addNextValue(currentRow, getResultState(entryIt->getColumn(), successorMemoryState), entryIt->getValue());
+                        }
+                        ++currentRow;
+                    }
+                }
+            }
+            
+            return builder.build();
+        }
+        
+        template <typename ValueType, typename RewardModelType>
+        storm::models::sparse::StateLabeling SparseModelMemoryProduct<ValueType, RewardModelType>::buildStateLabeling(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix) const {
+            uint64_t modelStateCount = model.getNumberOfStates();
+            uint64_t memoryStateCount = memory.getNumberOfStates();
             
-            uint_fast64_t numResStates = resultTransitionMatrix.getRowGroupCount();
+            uint64_t numResStates = resultTransitionMatrix.getRowGroupCount();
             storm::models::sparse::StateLabeling resultLabeling(numResStates);
             
             for (std::string modelLabel : model.getStateLabeling().getLabels()) {
                 if (modelLabel != "init") {
                     storm::storage::BitVector resLabeledStates(numResStates, false);
                     for (auto const& modelState : model.getStateLabeling().getStates(modelLabel)) {
-                        for (uint_fast64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
-                            uint_fast64_t const& resState = getResultState(modelState, memoryState);
+                        for (uint64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
+                            uint64_t const& resState = getResultState(modelState, memoryState);
                             // Check if the state exists in the result (i.e. if it is reachable)
                             if (resState < numResStates) {
                                 resLabeledStates.set(resState, true);
@@ -204,8 +338,8 @@ namespace storm {
                 STORM_LOG_THROW(!resultLabeling.containsLabel(memoryLabel), storm::exceptions::InvalidOperationException, "Failed to build the product of model and memory structure: State labelings are not disjoint as both structures contain the label " << memoryLabel << ".");
                 storm::storage::BitVector resLabeledStates(numResStates, false);
                 for (auto const& memoryState : memory.getStateLabeling().getStates(memoryLabel)) {
-                    for (uint_fast64_t modelState = 0; modelState < modelStateCount; ++modelState) {
-                        uint_fast64_t const& resState = getResultState(modelState, memoryState);
+                    for (uint64_t modelState = 0; modelState < modelStateCount; ++modelState) {
+                        uint64_t const& resState = getResultState(modelState, memoryState);
                         // Check if the state exists in the result (i.e. if it is reachable)
                         if (resState < numResStates) {
                             resLabeledStates.set(resState, true);
@@ -217,21 +351,24 @@ namespace storm {
             return resultLabeling;
         }
         
-        template <typename ValueType>
-        std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>> SparseModelMemoryProduct<ValueType>::buildRewardModels(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix, std::vector<uint_fast64_t> const& memorySuccessors) const {
-            std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>> result;
-            uint_fast64_t memoryStateCount = memory.getNumberOfStates();
-            uint_fast64_t numResStates = resultTransitionMatrix.getRowGroupCount();
+        template <typename ValueType, typename RewardModelType>
+        std::unordered_map<std::string, RewardModelType> SparseModelMemoryProduct<ValueType, RewardModelType>::buildRewardModels(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix, std::vector<uint64_t> const& memorySuccessors, boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) const {
+            
+            typedef typename RewardModelType::ValueType RewardValueType;
+            
+            std::unordered_map<std::string, RewardModelType> result;
+            uint64_t memoryStateCount = memory.getNumberOfStates();
+            uint64_t numResStates = resultTransitionMatrix.getRowGroupCount();
 
             for (auto const& rewardModel : model.getRewardModels()) {
-                boost::optional<std::vector<ValueType>> stateRewards;
+                boost::optional<std::vector<RewardValueType>> stateRewards;
                 if (rewardModel.second.hasStateRewards()) {
-                    stateRewards = std::vector<ValueType>(numResStates, storm::utility::zero<ValueType>());
-                    uint_fast64_t modelState = 0;
+                    stateRewards = std::vector<RewardValueType>(numResStates, storm::utility::zero<RewardValueType>());
+                    uint64_t modelState = 0;
                     for (auto const& modelStateReward : rewardModel.second.getStateRewardVector()) {
                         if (!storm::utility::isZero(modelStateReward)) {
-                            for (uint_fast64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
-                                uint_fast64_t const& resState = getResultState(modelState, memoryState);
+                            for (uint64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
+                                uint64_t const& resState = getResultState(modelState, memoryState);
                                 // Check if the state exists in the result (i.e. if it is reachable)
                                 if (resState < numResStates) {
                                     stateRewards.get()[resState] = modelStateReward;
@@ -241,45 +378,78 @@ namespace storm {
                         ++modelState;
                     }
                 }
-                boost::optional<std::vector<ValueType>> stateActionRewards;
+                boost::optional<std::vector<RewardValueType>> stateActionRewards;
                 if (rewardModel.second.hasStateActionRewards()) {
-                    stateActionRewards = std::vector<ValueType>(resultTransitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
-                    uint_fast64_t modelState = 0;
-                    uint_fast64_t modelRow = 0;
+                    stateActionRewards = std::vector<RewardValueType>(resultTransitionMatrix.getRowCount(), storm::utility::zero<RewardValueType>());
+                    uint64_t modelState = 0;
+                    uint64_t modelRow = 0;
                     for (auto const& modelStateActionReward : rewardModel.second.getStateActionRewardVector()) {
                         if (!storm::utility::isZero(modelStateActionReward)) {
                             while (modelRow >= model.getTransitionMatrix().getRowGroupIndices()[modelState + 1]) {
                                 ++modelState;
                             }
-                            uint_fast64_t rowOffset = modelRow - model.getTransitionMatrix().getRowGroupIndices()[modelState];
-                            for (uint_fast64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
-                                uint_fast64_t const& resState = getResultState(modelState, memoryState);
+                            uint64_t rowOffset = modelRow - model.getTransitionMatrix().getRowGroupIndices()[modelState];
+                            for (uint64_t memoryState = 0; memoryState < memoryStateCount; ++memoryState) {
+                                uint64_t const& resState = getResultState(modelState, memoryState);
                                 // Check if the state exists in the result (i.e. if it is reachable)
                                 if (resState < numResStates) {
-                                    stateActionRewards.get()[resultTransitionMatrix.getRowGroupIndices()[resState] + rowOffset] = modelStateActionReward;
+                                    if (scheduler && scheduler->getChoice(modelState, memoryState).isDefined()) {
+                                        ValueType factor = scheduler->getChoice(modelState, memoryState).getChoiceAsDistribution().getProbability(rowOffset);
+                                        stateActionRewards.get()[resultTransitionMatrix.getRowGroupIndices()[resState]] = factor * modelStateActionReward;
+                                    } else {
+                                        stateActionRewards.get()[resultTransitionMatrix.getRowGroupIndices()[resState] + rowOffset] = modelStateActionReward;
+                                    }
                                 }
                             }
                         }
                         ++modelRow;
                     }
                 }
-                boost::optional<storm::storage::SparseMatrix<ValueType>> transitionRewards;
+                boost::optional<storm::storage::SparseMatrix<RewardValueType>> transitionRewards;
                 if (rewardModel.second.hasTransitionRewards()) {
-                    storm::storage::SparseMatrixBuilder<ValueType> builder(resultTransitionMatrix.getRowCount(), resultTransitionMatrix.getColumnCount());
-                    uint_fast64_t stateIndex = 0;
+                    storm::storage::SparseMatrixBuilder<RewardValueType> builder(resultTransitionMatrix.getRowCount(), resultTransitionMatrix.getColumnCount());
+                    uint64_t stateIndex = 0;
                     for (auto const& resState : toResultStateMapping) {
                         if (resState < numResStates) {
-                            uint_fast64_t modelState = stateIndex / memoryStateCount;
-                            uint_fast64_t memoryState = stateIndex % memoryStateCount;
-                            uint_fast64_t rowGroupSize = resultTransitionMatrix.getRowGroupSize(resState);
-                            for (uint_fast64_t rowOffset = 0; rowOffset < rowGroupSize; ++rowOffset) {
-                                uint_fast64_t resRowIndex = resultTransitionMatrix.getRowGroupIndices()[resState] + rowOffset;
-                                uint_fast64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState] + rowOffset;
-                                auto const& rewardMatrixRow = rewardModel.second.getTransitionRewardMatrix().getRow(modelRowIndex);
-                                for (auto entryIt = rewardMatrixRow.begin(); entryIt != rewardMatrixRow.end(); ++entryIt) {
-                                    uint_fast64_t transitionId = entryIt - rewardModel.second.getTransitionRewardMatrix().begin();
-                                    uint_fast64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
-                                    builder.addNextValue(resRowIndex, getResultState(entryIt->getColumn(), successorMemoryState), entryIt->getValue());
+                            uint64_t modelState = stateIndex / memoryStateCount;
+                            uint64_t memoryState = stateIndex % memoryStateCount;
+                            uint64_t rowGroupSize = resultTransitionMatrix.getRowGroupSize(resState);
+                            if (scheduler && scheduler->getChoice(modelState, memoryState).isDefined()) {
+                                std::map<uint64_t, RewardValueType> rewards;
+                                for (uint64_t rowOffset = 0; rowOffset < rowGroupSize; ++rowOffset) {
+                                    uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState] + rowOffset;
+                                    auto transitionEntryIt = model.getTransitionMatrix().getRow(modelRowIndex).begin();
+                                    for (auto const& rewardEntry : rewardModel.second.getTransitionRewardMatrix().getRow(modelRowIndex)) {
+                                        while (transitionEntryIt->getColumn() != rewardEntry.getColumn()) {
+                                            STORM_LOG_ASSERT(transitionEntryIt != model.getTransitionMatrix().getRow(modelRowIndex).end(), "The reward transition matrix is not a submatrix of the model transition matrix.");
+                                            ++transitionEntryIt;
+                                        }
+                                        uint64_t transitionId = transitionEntryIt - model.getTransitionMatrix().begin();
+                                        uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                                        auto insertionRes = rewards.insert(std::make_pair(getResultState(rewardEntry.getColumn(), successorMemoryState), rewardEntry.getValue()));
+                                        if (!insertionRes.second) {
+                                            insertionRes.first->second += rewardEntry.getValue();
+                                        }
+                                    }
+                                }
+                                uint64_t resRowIndex = resultTransitionMatrix.getRowGroupIndices()[resState];
+                                for (auto& reward : rewards) {
+                                    builder.addNextValue(resRowIndex, reward.first, reward.second);
+                                }
+                            } else {
+                                for (uint64_t rowOffset = 0; rowOffset < rowGroupSize; ++rowOffset) {
+                                    uint64_t resRowIndex = resultTransitionMatrix.getRowGroupIndices()[resState] + rowOffset;
+                                    uint64_t modelRowIndex = model.getTransitionMatrix().getRowGroupIndices()[modelState] + rowOffset;
+                                    auto transitionEntryIt = model.getTransitionMatrix().getRow(modelRowIndex).begin();
+                                    for (auto const& rewardEntry : rewardModel.second.getTransitionRewardMatrix().getRow(modelRowIndex)) {
+                                        while (transitionEntryIt->getColumn() != rewardEntry.getColumn()) {
+                                            STORM_LOG_ASSERT(transitionEntryIt != model.getTransitionMatrix().getRow(modelRowIndex).end(), "The reward transition matrix is not a submatrix of the model transition matrix.");
+                                            ++transitionEntryIt;
+                                        }
+                                        uint64_t transitionId = transitionEntryIt - model.getTransitionMatrix().begin();
+                                        uint64_t const& successorMemoryState = memorySuccessors[transitionId * memoryStateCount + memoryState];
+                                        builder.addNextValue(resRowIndex, getResultState(rewardEntry.getColumn(), successorMemoryState), rewardEntry.getValue());
+                                    }
                                 }
                             }
                         }
@@ -287,32 +457,32 @@ namespace storm {
                     }
                     transitionRewards = builder.build();
                 }
-                result.insert(std::make_pair(rewardModel.first, storm::models::sparse::StandardRewardModel<ValueType>(std::move(stateRewards), std::move(stateActionRewards), std::move(transitionRewards))));
+                result.insert(std::make_pair(rewardModel.first, RewardModelType(std::move(stateRewards), std::move(stateActionRewards), std::move(transitionRewards))));
             }
             return result;
         }
             
-        template <typename ValueType>
-        std::shared_ptr<storm::models::sparse::Model<ValueType>> SparseModelMemoryProduct<ValueType>::buildResult(storm::storage::SparseMatrix<ValueType>&& matrix, storm::models::sparse::StateLabeling&& labeling, std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>>&& rewardModels) const {
-            storm::storage::sparse::ModelComponents<ValueType, storm::models::sparse::StandardRewardModel<ValueType>> components (std::move(matrix), std::move(labeling), std::move(rewardModels));
+        template <typename ValueType, typename RewardModelType>
+        std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> SparseModelMemoryProduct<ValueType, RewardModelType>::buildResult(storm::storage::SparseMatrix<ValueType>&& matrix, storm::models::sparse::StateLabeling&& labeling, std::unordered_map<std::string, RewardModelType>&& rewardModels, boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) const {
+            storm::storage::sparse::ModelComponents<ValueType, RewardModelType> components (std::move(matrix), std::move(labeling), std::move(rewardModels));
             
             if (model.isOfType(storm::models::ModelType::Ctmc)) {
                 components.rateTransitions = true;
             } else if (model.isOfType(storm::models::ModelType::MarkovAutomaton)) {
                 // We also need to translate the exit rates and the Markovian states
-                uint_fast64_t numResStates = components.transitionMatrix.getRowGroupCount();
-                uint_fast64_t memoryStateCount = memory.getNumberOfStates();
+                uint64_t numResStates = components.transitionMatrix.getRowGroupCount();
+                uint64_t memoryStateCount = memory.getNumberOfStates();
                 std::vector<ValueType> resultExitRates;
                 resultExitRates.reserve(components.transitionMatrix.getRowGroupCount());
                 storm::storage::BitVector resultMarkovianStates(numResStates, false);
-                auto const& modelExitRates = dynamic_cast<storm::models::sparse::MarkovAutomaton<ValueType> const&>(model).getExitRates();
-                auto const& modelMarkovianStates = dynamic_cast<storm::models::sparse::MarkovAutomaton<ValueType> const&>(model).getMarkovianStates();
+                auto const& modelExitRates = dynamic_cast<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType> const&>(model).getExitRates();
+                auto const& modelMarkovianStates = dynamic_cast<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType> const&>(model).getMarkovianStates();
                     
-                uint_fast64_t stateIndex = 0;
+                uint64_t stateIndex = 0;
                 for (auto const& resState : toResultStateMapping) {
                     if (resState < numResStates) {
                         assert(resState == resultExitRates.size());
-                        uint_fast64_t modelState = stateIndex / memoryStateCount;
+                        uint64_t modelState = stateIndex / memoryStateCount;
                         resultExitRates.push_back(modelExitRates[modelState]);
                         if (modelMarkovianStates.get(modelState)) {
                             resultMarkovianStates.set(resState, true);
@@ -324,12 +494,21 @@ namespace storm {
                 components.exitRates = std::move(resultExitRates);
             }
             
-            return storm::utility::builder::buildModelFromComponents(model.getType(), std::move(components));
+            storm::models::ModelType resultType = model.getType();
+            if (scheduler && !scheduler->isPartialScheduler()) {
+                if (model.isOfType(storm::models::ModelType::Mdp)) {
+                    resultType = storm::models::ModelType::Dtmc;
+                }
+                // Note that converting deterministic MAs to CTMCs via state elimination does not preserve all properties (e.g. step bounded)
+            }
+            
+            return storm::utility::builder::buildModelFromComponents(resultType, std::move(components));
         }
         
-        template  class SparseModelMemoryProduct<double>;
-        template  class SparseModelMemoryProduct<storm::RationalNumber>;
-        template  class SparseModelMemoryProduct<storm::RationalFunction>;
+        template class SparseModelMemoryProduct<double>;
+        template class SparseModelMemoryProduct<double, storm::models::sparse::StandardRewardModel<storm::Interval>>;
+        template class SparseModelMemoryProduct<storm::RationalNumber>;
+        template class SparseModelMemoryProduct<storm::RationalFunction>;
             
     }
 }
diff --git a/src/storm/storage/memorystructure/SparseModelMemoryProduct.h b/src/storm/storage/memorystructure/SparseModelMemoryProduct.h
index 0c60923e7..fb959daa5 100644
--- a/src/storm/storage/memorystructure/SparseModelMemoryProduct.h
+++ b/src/storm/storage/memorystructure/SparseModelMemoryProduct.h
@@ -8,61 +8,71 @@
 #include "storm/storage/memorystructure/MemoryStructure.h"
 #include "storm/models/sparse/Model.h"
 #include "storm/models/sparse/StandardRewardModel.h"
+#include "storm/storage/Scheduler.h"
 
 namespace storm {
     namespace storage {
         /*!
          * This class builds the product of the given sparse model and the given memory structure.
          * This is similar to the well-known product of a model with a deterministic rabin automaton.
-         * Note: we already do one memory-structure-step for the initial state, i.e., if s is an initial state of
-         * the given model and s satisfies memoryStructure.getTransitionMatrix[0][n], then (s,n) will be the corresponding
-         * initial state of the resulting model.
+         * The product contains only the reachable states of the product
          *
          * The states of the resulting sparse model will have the original state labels plus the labels of this
          * memory structure.
          * An exception is thrown if the state labelings are not disjoint.
          */
-        template <typename ValueType>
+        template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
         class SparseModelMemoryProduct {
         public:
             
-            SparseModelMemoryProduct(storm::models::sparse::Model<ValueType> const& sparseModel, storm::storage::MemoryStructure const& memoryStructure);
+            SparseModelMemoryProduct(storm::models::sparse::Model<ValueType, RewardModelType> const& sparseModel, storm::storage::MemoryStructure const& memoryStructure);
             
-            // Invokes the building of the product
-            std::shared_ptr<storm::models::sparse::Model<ValueType>> build();
+            // Enforces that the given model and memory state as well as the successor(s) are considered reachable -- even if they are not reachable from an initial state.
+            void addReachableState(uint64_t const& modelState, uint64_t const& memoryState);
+    
+            // Enforces that every state is considered reachable. If this is set, the result has size #modelStates * #memoryStates
+            void setBuildFullProduct();
+            
+            // Invokes the building of the product under the specified scheduler (if given).
+            std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> build(boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler = boost::none);
             
             // Retrieves the state of the resulting model that represents the given memory and model state.
-            // An invalid index is returned, if the specifyied state does not exist (i.e., if it is not reachable).
-            uint_fast64_t const& getResultState(uint_fast64_t const& modelState, uint_fast64_t const& memoryState) const;
+            // Should only be called AFTER calling build();
+            // An invalid index is returned, if the specifyied state does not exist (i.e., if it is not part of product).
+            uint64_t const& getResultState(uint64_t const& modelState, uint64_t const& memoryState) const;
             
         private:
             
             // Computes for each pair of memory state and model transition index the successor memory state
             // The resulting vector maps (modelTransition * memoryStateCount) + memoryState to the corresponding successor state of the memory structure
-            std::vector<uint_fast64_t> computeMemorySuccessors() const;
+            std::vector<uint64_t> computeMemorySuccessors() const;
             
             // Computes the reachable states of the resulting model
-            storm::storage::BitVector computeReachableStates(std::vector<uint_fast64_t> const& memorySuccessors, storm::storage::BitVector const& initialStates) const;
+            void computeReachableStates(std::vector<uint64_t> const& memorySuccessors, storm::storage::BitVector const& initialStates, boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler);
             
             // Methods that build the model components
             // Matrix for deterministic models
-            storm::storage::SparseMatrix<ValueType> buildDeterministicTransitionMatrix(storm::storage::BitVector const& reachableStates, std::vector<uint_fast64_t> const& memorySuccessors) const;
+            storm::storage::SparseMatrix<ValueType> buildDeterministicTransitionMatrix(std::vector<uint64_t> const& memorySuccessors) const;
             // Matrix for nondeterministic models
-            storm::storage::SparseMatrix<ValueType> buildNondeterministicTransitionMatrix(storm::storage::BitVector const& reachableStates, std::vector<uint_fast64_t> const& memorySuccessors) const;
+            storm::storage::SparseMatrix<ValueType> buildNondeterministicTransitionMatrix(std::vector<uint64_t> const& memorySuccessors) const;
+            // Matrix for models that consider a scheduler
+            storm::storage::SparseMatrix<ValueType> buildTransitionMatrixForScheduler(std::vector<uint64_t> const& memorySuccessors, storm::storage::Scheduler<ValueType> const& scheduler) const;
             // State labeling. Note: DOES NOT ADD A LABEL FOR THE INITIAL STATES
             storm::models::sparse::StateLabeling buildStateLabeling(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix) const;
             // Reward models
-            std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>> buildRewardModels(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix, std::vector<uint_fast64_t> const& memorySuccessors) const;
+            std::unordered_map<std::string, RewardModelType> buildRewardModels(storm::storage::SparseMatrix<ValueType> const& resultTransitionMatrix, std::vector<uint64_t> const& memorySuccessors, boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) const;
             
             // Builds the resulting model
-            std::shared_ptr<storm::models::sparse::Model<ValueType>> buildResult(storm::storage::SparseMatrix<ValueType>&& matrix, storm::models::sparse::StateLabeling&& labeling, std::unordered_map<std::string, storm::models::sparse::StandardRewardModel<ValueType>>&& rewardModels) const;
+            std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> buildResult(storm::storage::SparseMatrix<ValueType>&& matrix, storm::models::sparse::StateLabeling&& labeling, std::unordered_map<std::string, RewardModelType>&& rewardModels, boost::optional<storm::storage::Scheduler<ValueType>> const& scheduler) const;
             
             
             // Maps (modelState * memoryStateCount) + memoryState to the state in the result that represents (memoryState,modelState)
-            std::vector<uint_fast64_t> toResultStateMapping;
+            std::vector<uint64_t> toResultStateMapping;
             
-
-            storm::models::sparse::Model<ValueType> const& model;
+            // Indicates which states are considered reachable. (s, m) is reachable if this BitVector is true at (s * memoryStateCount) + m
+            storm::storage::BitVector reachableStates;
+            
+            storm::models::sparse::Model<ValueType, RewardModelType> const& model;
             storm::storage::MemoryStructure const& memory;
             
         };