diff --git a/examples/ctmc/tandem/tandem.sm b/examples/ctmc/tandem/tandem.sm
index 439d33813..c35246e1a 100644
--- a/examples/ctmc/tandem/tandem.sm
+++ b/examples/ctmc/tandem/tandem.sm
@@ -35,4 +35,8 @@ endmodule
 // reward - number of customers in network
 rewards "customers"
 	true : sc + sm;
-endrewards
\ No newline at end of file
+endrewards
+
+label "network_full" = sc=c&sm=c&ph=2;
+label "first_queue_full" = sc=c;
+label "second_queue_full" = sm=c;
diff --git a/src/builder/DdPrismModelBuilder.cpp b/src/builder/DdPrismModelBuilder.cpp
index 6378a3c56..eae18d06e 100644
--- a/src/builder/DdPrismModelBuilder.cpp
+++ b/src/builder/DdPrismModelBuilder.cpp
@@ -562,10 +562,12 @@ namespace storm {
                 }
                 
                 storm::dd::Add<Type> transitionRewardDd = synchronization * states * rewards;
-                if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::MDP) {
-                    transitionRewardDd = transitions.notZero().toAdd() * transitionRewardDd;
-                } else {
+                if (generationInfo.program.getModelType() == storm::prism::Program::ModelType::DTMC) {
+                    // For DTMCs we need to keep the weighting for the scaling that follows.
                     transitionRewardDd = transitions * transitionRewardDd;
+                } else {
+                    // For all other model types, we do not scale the rewards.
+                    transitionRewardDd = transitions.notZero().toAdd() * transitionRewardDd;
                 }
                 
                 // Perform some sanity checks.
diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp
index 05cf9c52e..4315bfdde 100644
--- a/src/builder/ExplicitPrismModelBuilder.cpp
+++ b/src/builder/ExplicitPrismModelBuilder.cpp
@@ -127,21 +127,6 @@ namespace storm {
 #endif
             }
             
-            storm::prism::RewardModel rewardModel = storm::prism::RewardModel();
-            
-            // Select the appropriate reward model.
-            if (options.buildRewards) {
-                // If a specific reward model was selected or one with the empty name exists, select it.
-                if (options.rewardModelName) {
-                    rewardModel = preparedProgram.getRewardModel(options.rewardModelName.get());
-                } else if (preparedProgram.hasRewardModel("")) {
-                    rewardModel = preparedProgram.getRewardModel("");
-                } else if (preparedProgram.hasRewardModel()) {
-                    // Otherwise, we select the first one.
-                    rewardModel = preparedProgram.getRewardModel(0);
-                }
-            }
-            
             // If the set of labels we are supposed to built is restricted, we need to remove the other labels from the program.
             if (options.labelsToBuild) {
                 preparedProgram.filterLabels(options.labelsToBuild.get());
@@ -162,6 +147,20 @@ namespace storm {
             
             // Now that the program is fixed, we we need to substitute all constants with their concrete value.
             preparedProgram = preparedProgram.substituteConstants();
+            
+            // Select the appropriate reward model (after the constants have been substituted).
+            storm::prism::RewardModel rewardModel = storm::prism::RewardModel();
+            if (options.buildRewards) {
+                // If a specific reward model was selected or one with the empty name exists, select it.
+                if (options.rewardModelName) {
+                    rewardModel = preparedProgram.getRewardModel(options.rewardModelName.get());
+                } else if (preparedProgram.hasRewardModel("")) {
+                    rewardModel = preparedProgram.getRewardModel("");
+                } else if (preparedProgram.hasRewardModel()) {
+                    // Otherwise, we select the first one.
+                    rewardModel = preparedProgram.getRewardModel(0);
+                }
+            }
                 
             ModelComponents modelComponents = buildModelComponents(preparedProgram, rewardModel, options);
             
diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
new file mode 100644
index 000000000..e371ad8ce
--- /dev/null
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
@@ -0,0 +1,364 @@
+#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
+#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
+#include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
+
+#include "src/storage/dd/CuddOdd.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/graph.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, class ValueType>
+        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::LinearEquationSolverFactory<ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        bool HybridCtmcCslModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const {
+            return formula.isCslStateFormula() || formula.isCslPathFormula() || formula.isRewardPathFormula();
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        storm::dd::Add<DdType> HybridCtmcCslModelChecker<DdType, ValueType>::computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector) {
+            return rateMatrix / exitRateVector;
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            
+            return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), HybridDtmcPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector())));
+        }
+
+        template<storm::dd::DdType DdType, class ValueType>
+        storm::models::symbolic::Ctmc<DdType> const& HybridCtmcCslModelChecker<DdType, ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::symbolic::Ctmc<DdType>>();
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        storm::dd::Add<DdType> HybridCtmcCslModelChecker<DdType, ValueType>::computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate) {
+            STORM_LOG_DEBUG("Computing uniformized matrix using uniformization rate " << uniformizationRate << ".");
+            STORM_LOG_DEBUG("Keeping " << maybeStates.getNonZeroCount() << " rows.");
+            
+            // Cut all non-maybe rows/columns from the transition matrix.
+            storm::dd::Add<DdType> uniformizedMatrix = transitionMatrix * maybeStates.toAdd() * maybeStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+            
+            // Now perform the uniformization.
+            uniformizedMatrix = uniformizedMatrix / model.getManager().getConstant(uniformizationRate);
+            storm::dd::Add<DdType> diagonal = model.getRowColumnIdentity() * maybeStates.toAdd();
+            storm::dd::Add<DdType> diagonalOffset = diagonal;
+            diagonalOffset -= diagonal * (exitRateVector / model.getManager().getConstant(uniformizationRate));
+            uniformizedMatrix += diagonalOffset;
+            
+            return uniformizedMatrix;
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+
+            boost::optional<storm::dd::Add<DdType>> modifiedStateRewardVector;
+            if (this->getModel().hasStateRewards()) {
+                modifiedStateRewardVector = this->getModel().getStateRewardVector() / this->getModel().getExitRateVector();
+            }
+            
+            return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative);
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            double lowerBound = 0;
+            double upperBound = 0;
+            if (!pathFormula.hasDiscreteTimeBound()) {
+                std::pair<double, double> const& intervalBounds =  pathFormula.getIntervalBounds();
+                lowerBound = intervalBounds.first;
+                upperBound = intervalBounds.second;
+            } else {
+                upperBound = pathFormula.getDiscreteTimeBound();
+            }
+            
+            return this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound);
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, storm::dd::Add<DdType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const {
+            // If the time bounds are [0, inf], we rather call untimed reachability.
+            storm::utility::ConstantsComparator<ValueType> comparator;
+            if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
+                return HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
+            }
+            
+            // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0
+            // or t' != inf.
+            
+            // If we identify the states that have probability 0 of reaching the target states, we can exclude them from the
+            // further computations.
+            storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), this->getModel().getTransitionMatrix().notZero(), phiStates, psiStates);
+            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0.getNonZeroCount() << " states with probability greater 0.");
+            storm::dd::Bdd<DdType> statesWithProbabilityGreater0NonPsi = statesWithProbabilityGreater0 && !psiStates;
+            STORM_LOG_INFO("Found " << statesWithProbabilityGreater0NonPsi.getNonZeroCount() << " 'maybe' states.");
+            
+            if (!statesWithProbabilityGreater0NonPsi.isZero()) {
+                if (comparator.isZero(upperBound)) {
+                    // In this case, the interval is of the form [0, 0].
+                    return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), psiStates.toAdd()));
+                } else {
+                    if (comparator.isZero(lowerBound)) {
+                        // In this case, the interval is of the form [0, t].
+                        // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier.
+                        
+                        // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
+                        ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRates).getMax();
+                        STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                        
+                        // Compute the uniformized matrix.
+                        storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0NonPsi, uniformizationRate);
+                        
+                        // Compute the vector that is to be added as a compensation for removing the absorbing states.
+                        storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * this->getModel().getTransitionMatrix() * psiStates.swapVariables(this->getModel().getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(this->getModel().getColumnVariables()) / this->getModel().getManager().getConstant(uniformizationRate);
+                        
+                        // Create an ODD for the translation to an explicit representation.
+                        storm::dd::Odd<DdType> odd(statesWithProbabilityGreater0NonPsi);
+                        
+                        // Convert the symbolic parts to their explicit representation.
+                        storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                        std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
+                        
+                        // Finally compute the transient probabilities.
+                        std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
+                        std::vector<ValueType> subresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound, uniformizationRate, values, *this->linearEquationSolverFactory);
+
+                        return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(),
+                                                                                                      (psiStates || !statesWithProbabilityGreater0) && this->getModel().getReachableStates(),
+                                                                                                      psiStates.toAdd(),
+                                                                                                      statesWithProbabilityGreater0NonPsi,
+                                                                                                      odd, subresult));
+                    } else if (comparator.isInfinity(upperBound)) {
+                        // In this case, the interval is of the form [t, inf] with t != 0.
+                        
+                        // Start by computing the (unbounded) reachability probabilities of reaching psi states while
+                        // staying in phi states.
+                        std::unique_ptr<CheckResult> unboundedResult = HybridDtmcPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(this->getModel(), this->computeProbabilityMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
+                        
+                        // Compute the set of relevant states.
+                        storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
+
+                        // Filter the unbounded result such that it only contains values for the relevant states.
+                        unboundedResult->filter(SymbolicQualitativeCheckResult<DdType>(this->getModel().getReachableStates(), relevantStates));
+
+                        // Build an ODD for the relevant states.
+                        storm::dd::Odd<DdType> odd(relevantStates);
+                        
+                        std::vector<ValueType> result;
+                        if (unboundedResult->isHybridQuantitativeCheckResult()) {
+                            std::unique_ptr<CheckResult> explicitUnboundedResult = unboundedResult->asHybridQuantitativeCheckResult<DdType>().toExplicitQuantitativeCheckResult();
+                            result = std::move(explicitUnboundedResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
+                        } else {
+                            STORM_LOG_THROW(unboundedResult->isSymbolicQuantitativeCheckResult(), storm::exceptions::InvalidStateException, "Expected check result of different type.");
+                            result = unboundedResult->asSymbolicQuantitativeCheckResult<DdType>().getValueVector().template toVector<ValueType>(odd);
+                        }
+                        
+                        // Determine the uniformization rate for the transient probability computation.
+                        ValueType uniformizationRate = 1.02 * (relevantStates.toAdd() * exitRates).getMax();
+                        
+                        // Compute the uniformized matrix.
+                        storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  relevantStates, uniformizationRate);
+                        storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                        
+                        // Compute the transient probabilities.
+                        result = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, result, *this->linearEquationSolverFactory);
+                        
+                        return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, result));
+                    } else {
+                        // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
+                        
+                        if (lowerBound != upperBound) {
+                            // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'.
+
+                            // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
+                            ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0NonPsi.toAdd() * exitRates).getMax();
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                        
+                            // Compute the (first) uniformized matrix.
+                            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0NonPsi, uniformizationRate);
+                            
+                            // Create the one-step vector.
+                            storm::dd::Add<DdType> b = (statesWithProbabilityGreater0NonPsi.toAdd() * this->getModel().getTransitionMatrix() * psiStates.swapVariables(this->getModel().getRowColumnMetaVariablePairs()).toAdd()).sumAbstract(this->getModel().getColumnVariables()) / this->getModel().getManager().getConstant(uniformizationRate);
+                            
+                            // Build an ODD for the relevant states and translate the symbolic parts to their explicit representation.
+                            storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0NonPsi);
+                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                            std::vector<ValueType> explicitB = b.template toVector<ValueType>(odd);
+                            
+                            // Compute the transient probabilities.
+                            std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNonZeroCount(), storm::utility::zero<ValueType>());
+                            std::vector<ValueType> subResult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, &explicitB, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
+
+                            // Transform the explicit result to a hybrid check result, so we can easily convert it to
+                            // a symbolic qualitative format.
+                            HybridQuantitativeCheckResult<DdType> hybridResult(this->getModel().getReachableStates(), psiStates || (!statesWithProbabilityGreater0 && this->getModel().getReachableStates()),
+                                                                               psiStates.toAdd(), statesWithProbabilityGreater0NonPsi, odd, subResult);
+                            
+                            // Compute the set of relevant states.
+                            storm::dd::Bdd<DdType> relevantStates = statesWithProbabilityGreater0 && phiStates;
+                            
+                            // Filter the unbounded result such that it only contains values for the relevant states.
+                            hybridResult.filter(SymbolicQualitativeCheckResult<DdType>(this->getModel().getReachableStates(), relevantStates));
+                            
+                            // Build an ODD for the relevant states.
+                            odd = storm::dd::Odd<DdType>(relevantStates);
+                            
+                            std::unique_ptr<CheckResult> explicitResult = hybridResult.toExplicitQuantitativeCheckResult();
+                            std::vector<ValueType> newSubresult = std::move(explicitResult->asExplicitQuantitativeCheckResult<ValueType>().getValueVector());
+                            
+                            // Then compute the transient probabilities of being in such a state after t time units. For this,
+                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                            uniformizationRate =  1.02 * (relevantStates.toAdd() * exitRates).getMax();
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                            
+                            // If the lower and upper bounds coincide, we have only determined the relevant states at this
+                            // point, but we still need to construct the starting vector.
+                            if (lowerBound == upperBound) {
+                                odd = storm::dd::Odd<DdType>(relevantStates);
+                                newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
+                            }
+                            
+                            // Finally, we compute the second set of transient probabilities.
+                            uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  relevantStates, uniformizationRate);
+                            explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                            
+                            newSubresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            
+                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !relevantStates && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), relevantStates, odd, newSubresult));
+                        } else {
+                            // In this case, the interval is of the form [t, t] with t != 0, t != inf.
+
+                            // Build an ODD for the relevant states.
+                            storm::dd::Odd<DdType> odd = storm::dd::Odd<DdType>(statesWithProbabilityGreater0);
+                            
+                            std::vector<ValueType> newSubresult = psiStates.toAdd().template toVector<ValueType>(odd);
+                            
+                            // Then compute the transient probabilities of being in such a state after t time units. For this,
+                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                            ValueType uniformizationRate =  1.02 * (statesWithProbabilityGreater0.toAdd() * exitRates).getMax();
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                            
+                            // Finally, we compute the second set of transient probabilities.
+                            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), exitRates,  statesWithProbabilityGreater0, uniformizationRate);
+                            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                            
+                            newSubresult = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            
+                            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), !statesWithProbabilityGreater0 && this->getModel().getReachableStates(), this->getModel().getManager().getAddZero(), statesWithProbabilityGreater0, odd, newSubresult));
+                        }
+                    }
+                }
+            } else {
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), psiStates.toAdd()));
+            }
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            return this->computeInstantaneousRewardsHelper(rewardPathFormula.getContinuousTimeBound());
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(double timeBound) const {
+            // Only compute the result if the model has a state-based reward this->getModel().
+            STORM_LOG_THROW(this->getModel().hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Create ODD for the translation.
+            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
+            
+            // Initialize result to state rewards of the this->getModel().
+            std::vector<ValueType> result = this->getModel().getStateRewardVector().template toVector<ValueType>(odd);
+            
+            // If the time-bound is not zero, we need to perform a transient analysis.
+            if (timeBound > 0) {
+                ValueType uniformizationRate = 1.02 * this->getModel().getExitRateVector().getMax();
+                STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                
+                storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(),  this->getModel().getReachableStates(), uniformizationRate);
+                
+                storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+                result = SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory);
+            }
+            
+            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), odd, result));
+        }
+
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            return this->computeCumulativeRewardsHelper(rewardPathFormula.getContinuousTimeBound());
+        }
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        std::unique_ptr<CheckResult> HybridCtmcCslModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(double timeBound) const {
+            // Only compute the result if the model has a state-based reward this->getModel().
+            STORM_LOG_THROW(this->getModel().hasStateRewards() || this->getModel().hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // If the time bound is zero, the result is the constant zero vector.
+            if (timeBound == 0) {
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getAddZero()));
+            }
+            
+            // Otherwise, we need to perform some computations.
+            
+            // Start with the uniformization.
+            ValueType uniformizationRate = 1.02 * this->getModel().getExitRateVector().getMax();
+            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+            
+            // Create ODD for the translation.
+            storm::dd::Odd<DdType> odd(this->getModel().getReachableStates());
+            
+            // Compute the uniformized matrix.
+            storm::dd::Add<DdType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector(),  this->getModel().getReachableStates(), uniformizationRate);
+            storm::storage::SparseMatrix<ValueType> explicitUniformizedMatrix = uniformizedMatrix.toMatrix(odd, odd);
+
+            // Then compute the state reward vector to use in the computation.
+            storm::dd::Add<DdType> totalRewardVector = this->getModel().hasStateRewards() ? this->getModel().getStateRewardVector() : this->getModel().getManager().getAddZero();
+            if (this->getModel().hasTransitionRewards()) {
+                totalRewardVector += (this->getModel().getTransitionMatrix() * this->getModel().getTransitionRewardMatrix()).sumAbstract(this->getModel().getColumnVariables());
+            }
+            std::vector<ValueType> explicitTotalRewardVector = totalRewardVector.template toVector<ValueType>(odd);
+            
+            // Finally, compute the transient probabilities.
+            std::vector<ValueType> result = SparseCtmcCslModelChecker<ValueType>::template computeTransientProbabilities<true>(explicitUniformizedMatrix, nullptr, timeBound, uniformizationRate, explicitTotalRewardVector, *this->linearEquationSolverFactory);
+            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->getModel().getManager().getBddZero(), this->getModel().getManager().getAddZero(), this->getModel().getReachableStates(), odd, result));
+        }
+        
+        // Explicitly instantiate the model checker.
+        template class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double>;
+        
+    } // namespace modelchecker
+} // namespace storm
\ No newline at end of file
diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.h b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
new file mode 100644
index 000000000..f0e619533
--- /dev/null
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
@@ -0,0 +1,65 @@
+#ifndef STORM_MODELCHECKER_HYBRIDCTMCCSLMODELCHECKER_H_
+#define STORM_MODELCHECKER_HYBRIDCTMCCSLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+#include "src/models/symbolic/Ctmc.h"
+#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
+
+namespace storm {
+    namespace modelchecker {
+        
+        template<storm::dd::DdType DdType, class ValueType>
+        class HybridCtmcCslModelChecker : public SymbolicPropositionalModelChecker<DdType> {
+        public:
+            explicit HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType> const& model);
+            explicit HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            
+            // The implemented methods of the AbstractModelChecker interface.
+            virtual bool canHandle(storm::logic::Formula const& formula) const override;
+            virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            
+        protected:
+            storm::models::symbolic::Ctmc<DdType> const& getModel() const override;
+            
+        private:
+            /*!
+             * Converts the given rate-matrix into a time-abstract probability matrix.
+             *
+             * @param model The symbolic model.
+             * @param rateMatrix The rate matrix.
+             * @param exitRateVector The exit rate vector of the model.
+             * @return The probability matrix.
+             */
+            static storm::dd::Add<DdType> computeProbabilityMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& rateMatrix, storm::dd::Add<DdType> const& exitRateVector);
+
+            /*!
+             * Computes the matrix representing the transitions of the uniformized CTMC.
+             *
+             * @param model The symbolic model.
+             * @param transitionMatrix The matrix to uniformize.
+             * @param exitRateVector The exit rate vector.
+             * @param maybeStates The states that need to be considered.
+             * @param uniformizationRate The rate to be used for uniformization.
+             * @return The uniformized matrix.
+             */
+            static storm::dd::Add<DdType> computeUniformizedMatrix(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Add<DdType> const& exitRateVector, storm::dd::Bdd<DdType> const& maybeStates, ValueType uniformizationRate);
+            
+            // The methods that perform the actual checking.
+            std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, storm::dd::Add<DdType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const;
+            std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(double timeBound) const;
+            std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(double timeBound) const;
+            
+            // An object that is used for solving linear equations and performing matrix-vector multiplication.
+            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_HYBRIDCTMCCSLMODELCHECKER_H_ */
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index 74540327b..ecd550f4c 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -49,9 +49,7 @@ namespace storm {
                 upperBound = pathFormula.getDiscreteTimeBound();
             }
             
-            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound)));
-            
-            return result;
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound)));
         }
         
         template<class ValueType>
@@ -140,8 +138,7 @@ namespace storm {
                         result = this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), backwardTransitions, phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
                         
                         // Determine the set of states that must be considered further.
-                        storm::storage::BitVector relevantStates = storm::utility::vector::filterGreaterZero(result);
-                        relevantStates = storm::utility::graph::performProbGreater0(backwardTransitions, phiStates, relevantStates & phiStates);
+                        storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates;
                         std::vector<ValueType> subResult(relevantStates.getNumberOfSetBits());
                         storm::utility::vector::selectVectorValues(subResult, relevantStates, result);
                         
@@ -164,17 +161,11 @@ namespace storm {
                     } else {
                         // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
                         
-                        // Prepare some variables that are used by the two following blocks.
-                        storm::storage::BitVector relevantStates;
-                        ValueType uniformizationRate = 0;
-                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix;
-                        std::vector<ValueType> newSubresult;
-                        
-                        if (lowerBound == upperBound) {
-                            relevantStates = statesWithProbabilityGreater0;
-                        } else {
+                        if (lowerBound != upperBound) {
+                            // In this case, the interval is of the form [t, t'] with t != 0, t' != inf and t != t'.
+                            
                             // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
-                            uniformizationRate = 0;
+                            ValueType uniformizationRate = storm::utility::zero<ValueType>();
                             for (auto const& state : statesWithProbabilityGreater0NonPsi) {
                                 uniformizationRate = std::max(uniformizationRate, exitRates[state]);
                             }
@@ -182,7 +173,7 @@ namespace storm {
                             STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
                             
                             // Compute the (first) uniformized matrix.
-                            uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates);
+                            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0NonPsi, uniformizationRate, exitRates);
                             
                             // Compute the vector that is to be added as a compensation for removing the absorbing states.
                             std::vector<ValueType> b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(statesWithProbabilityGreater0NonPsi, psiStates);
@@ -192,40 +183,54 @@ namespace storm {
                             
                             // Start by computing the transient probabilities of reaching a psi state in time t' - t.
                             std::vector<ValueType> values(statesWithProbabilityGreater0NonPsi.getNumberOfSetBits(), storm::utility::zero<ValueType>());
-                            std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
+                            std::vector<ValueType> subresult = computeTransientProbabilities(uniformizedMatrix, &b, upperBound - lowerBound, uniformizationRate, values, *this->linearEquationSolverFactory);
                             
-                            // Determine the set of states that must be considered further.
-                            relevantStates = storm::utility::vector::filterGreaterZero(subresult);
-                            relevantStates = storm::utility::graph::performProbGreater0(uniformizedMatrix.transpose(), phiStates % statesWithProbabilityGreater0NonPsi, relevantStates & (phiStates % statesWithProbabilityGreater0NonPsi));
+                            storm::storage::BitVector relevantStates = statesWithProbabilityGreater0 & phiStates;
+                            std::vector<ValueType> newSubresult = std::vector<ValueType>(relevantStates.getNumberOfSetBits());
+                            storm::utility::vector::setVectorValues(newSubresult, statesWithProbabilityGreater0NonPsi % relevantStates, subresult);
+                            storm::utility::vector::setVectorValues(newSubresult, psiStates % relevantStates, storm::utility::one<ValueType>());
+
+                            // Then compute the transient probabilities of being in such a state after t time units. For this,
+                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                            uniformizationRate = storm::utility::zero<ValueType>();
+                            for (auto const& state : relevantStates) {
+                                uniformizationRate = std::max(uniformizationRate, exitRates[state]);
+                            }
+                            uniformizationRate *= 1.02;
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
                             
-                            newSubresult = std::vector<ValueType>(relevantStates.getNumberOfSetBits());
-                            storm::utility::vector::selectVectorValues(newSubresult, relevantStates, subresult);
-                        }
-                        
-                        // Then compute the transient probabilities of being in such a state after t time units. For this,
-                        // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
-                        uniformizationRate = 0;
-                        for (auto const& state : relevantStates) {
-                            uniformizationRate = std::max(uniformizationRate, exitRates[state]);
-                        }
-                        uniformizationRate *= 1.02;
-                        STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
-                        
-                        // If the lower and upper bounds coincide, we have only determined the relevant states at this
-                        // point, but we still need to construct the starting vector.
-                        if (lowerBound == upperBound) {
-                            newSubresult = std::vector<ValueType>(relevantStates.getNumberOfSetBits());
+                            // Finally, we compute the second set of transient probabilities.
+                            uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), relevantStates, uniformizationRate, exitRates);
+                            newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            
+                            // Fill in the correct values.
+                            result = std::vector<ValueType>(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
+                            storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero<ValueType>());
+                            storm::utility::vector::setVectorValues(result, relevantStates, newSubresult);
+                        } else {
+                            // In this case, the interval is of the form [t, t] with t != 0, t != inf.
+                            
+                            std::vector<ValueType> newSubresult = std::vector<ValueType>(statesWithProbabilityGreater0.getNumberOfSetBits());
                             storm::utility::vector::setVectorValues(newSubresult, psiStates % statesWithProbabilityGreater0, storm::utility::one<ValueType>());
+                            
+                            // Then compute the transient probabilities of being in such a state after t time units. For this,
+                            // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
+                            ValueType uniformizationRate = storm::utility::zero<ValueType>();
+                            for (auto const& state : statesWithProbabilityGreater0) {
+                                uniformizationRate = std::max(uniformizationRate, exitRates[state]);
+                            }
+                            uniformizationRate *= 1.02;
+                            STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                            
+                            // Finally, we compute the second set of transient probabilities.
+                            storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, uniformizationRate, exitRates);
+                            newSubresult = computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
+                            
+                            // Fill in the correct values.
+                            result = std::vector<ValueType>(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
+                            storm::utility::vector::setVectorValues(result, ~statesWithProbabilityGreater0, storm::utility::zero<ValueType>());
+                            storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, newSubresult);
                         }
-                        
-                        // Finally, we compute the second set of transient probabilities.
-                        uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), relevantStates, uniformizationRate, exitRates);
-                        newSubresult = this->computeTransientProbabilities(uniformizedMatrix, nullptr, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
-                        
-                        // Fill in the correct values.
-                        result = std::vector<ValueType>(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
-                        storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero<ValueType>());
-                        storm::utility::vector::setVectorValues(result, relevantStates, newSubresult);
                     }
                 }
             }
@@ -249,7 +254,7 @@ namespace storm {
             for (auto const& state : maybeStates) {
                 for (auto& element : uniformizedMatrix.getRow(currentRow)) {
                     if (element.getColumn() == currentRow) {
-                        element.setValue(-(exitRates[state] - element.getValue()) / uniformizationRate + storm::utility::one<ValueType>());
+                        element.setValue((element.getValue() - exitRates[state]) / uniformizationRate + storm::utility::one<ValueType>());
                     } else {
                         element.setValue(element.getValue() / uniformizationRate);
                     }
@@ -262,7 +267,7 @@ namespace storm {
         
         template<class ValueType>
         template<bool computeCumulativeReward>
-        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) const {
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
             
             ValueType lambda = timeBound * uniformizationRate;
             
@@ -307,7 +312,7 @@ namespace storm {
                 if (computeCumulativeReward) {
                     result = std::vector<ValueType>(values.size());
                     std::function<ValueType (ValueType const&)> scaleWithUniformizationRate = [&uniformizationRate] (ValueType const& a) -> ValueType { return a / uniformizationRate; };
-                    storm::utility::vector::applyPointwise(result, result, scaleWithUniformizationRate);
+                    storm::utility::vector::applyPointwise(values, result, scaleWithUniformizationRate);
                 } else {
                     result = std::vector<ValueType>(values.size());
                 }
@@ -339,7 +344,7 @@ namespace storm {
                 weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
                 storm::utility::vector::applyPointwise(result, values, result, addAndScale);
             }
-            
+
             return result;
         }
         
@@ -440,6 +445,7 @@ namespace storm {
             boost::optional<std::vector<ValueType>> modifiedStateRewardVector;
             if (this->getModel().hasStateRewards()) {
                 modifiedStateRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
+                
                 typename std::vector<ValueType>::const_iterator it2 = this->getModel().getExitRateVector().begin();
                 for (typename std::vector<ValueType>::iterator it1 = modifiedStateRewardVector.get().begin(), ite1 = modifiedStateRewardVector.get().end(); it1 != ite1; ++it1, ++it2) {
                     *it1 /= *it2;
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
index b74967ac5..5108a183f 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -3,15 +3,21 @@
 
 #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
 #include "src/models/sparse/Ctmc.h"
+#include "src/storage/dd/DdType.h"
 #include "src/utility/solver.h"
 #include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
         
+        template<storm::dd::DdType DdType, typename ValueType>
+        class HybridCtmcCslModelChecker;
+        
         template<class ValueType>
         class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
+            friend class HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, ValueType>;
+            
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model);
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
@@ -33,7 +39,6 @@ namespace storm {
             static std::vector<ValueType> computeUntilProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             std::vector<ValueType> computeInstantaneousRewardsHelper(double timeBound) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const;
-            std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const;
 
             /*!
              * Computes the matrix representing the transitions of the uniformized CTMC.
@@ -61,7 +66,7 @@ namespace storm {
              * @return The vector of transient probabilities.
              */
             template<bool useMixedPoissonProbabilities = false>
-            std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) const;
+            static std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             
             /*!
              * Converts the given rate-matrix into a time-abstract probability matrix.
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
index 4ddda3339..25a02adff 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
@@ -231,13 +231,14 @@ namespace storm {
         std::unique_ptr<CheckResult> HybridDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
-            return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
+            return this->computeReachabilityRewardsHelper(this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<CheckResult> HybridDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
-            // Only compute the result if the model has at least one reward model.
-            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+        std::unique_ptr<CheckResult> HybridDtmcPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
+            
+            // Only compute the result if there is at least one reward model.
+            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
 
             // Determine which states have a reward of infinity by definition.
             storm::dd::Bdd<DdType> infinityStates = storm::utility::graph::performProb1(model, transitionMatrix.notZero(), model.getReachableStates(), targetStates);
@@ -246,7 +247,7 @@ namespace storm {
             STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
             STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
             STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
-
+            
             // Check whether we need to compute exact rewards for some states.
             if (qualitative) {
                 // Set the values for all maybe-states to 1 to indicate that their reward values
@@ -266,9 +267,9 @@ namespace storm {
                     storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
                     
                     // Then compute the state reward vector to use in the computation.
-                    storm::dd::Add<DdType> subvector = model.hasStateRewards() ? maybeStatesAdd * model.getStateRewardVector() : model.getManager().getAddZero();
-                    if (model.hasTransitionRewards()) {
-                        subvector += (submatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                    if (transitionRewardMatrix) {
+                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
                     }
                     
                     // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
index bb7a263a9..b3ab51279 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
@@ -7,9 +7,14 @@
 
 namespace storm {
     namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class HybridCtmcCslModelChecker;
+        
         template<storm::dd::DdType DdType, typename ValueType>
         class HybridDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
         public:
+            friend class HybridCtmcCslModelChecker<DdType, ValueType>;
+            
             explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model);
             explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
@@ -32,7 +37,7 @@ namespace storm {
             static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
-            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
+            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(storm::models::symbolic::Model<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
 
             // An object that is used for retrieving linear equation solvers.
             std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp
new file mode 100644
index 000000000..f54bed203
--- /dev/null
+++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.cpp
@@ -0,0 +1,330 @@
+#include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h"
+
+#include "src/storage/dd/CuddOdd.h"
+
+#include "src/utility/macros.h"
+#include "src/utility/graph.h"
+
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        HybridMdpPrctlModelChecker<DdType, ValueType>::HybridMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        HybridMdpPrctlModelChecker<DdType, ValueType>::HybridMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model) : SymbolicPropositionalModelChecker<DdType>(model), linearEquationSolverFactory(new storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>()) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        bool HybridMdpPrctlModelChecker<DdType, ValueType>::canHandle(storm::logic::Formula const& formula) const {
+            return formula.isPctlStateFormula() || formula.isPctlPathFormula() || formula.isRewardPathFormula();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 and 1 of satisfying the until-formula.
+            std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01;
+            if (minimize) {
+                statesWithProbability01 = storm::utility::graph::performProb01Min(model, phiStates, psiStates);
+            } else {
+                statesWithProbability01 = storm::utility::graph::performProb01Max(model, phiStates, psiStates);
+            }
+            storm::dd::Bdd<DdType> maybeStates = !statesWithProbability01.first && !statesWithProbability01.second && model.getReachableStates();
+            
+            // Perform some logging.
+            STORM_LOG_INFO("Found " << statesWithProbability01.first.getNonZeroCount() << " 'no' states.");
+            STORM_LOG_INFO("Found " << statesWithProbability01.second.getNonZeroCount() << " 'yes' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact probabilities for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 0.5 to indicate that their probability values are neither 0 nor 1.
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd() + maybeStates.toAdd() * model.getManager().getConstant(0.5)));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                    // maybe states.
+                    storm::dd::Add<DdType> prob1StatesAsColumn = statesWithProbability01.second.toAdd();
+                    prob1StatesAsColumn = prob1StatesAsColumn.swapVariables(model.getRowColumnMetaVariablePairs());
+                    storm::dd::Add<DdType> subvector = submatrix * prob1StatesAsColumn;
+                    subvector = subvector.sumAbstract(model.getColumnVariables());
+                    
+                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                    
+                    // Finally cut away all columns targeting non-maybe states.
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    
+                    // Create the solution vector.
+                    std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
+                    
+                    // Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
+                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                    
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                    solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
+                    
+                    // Return a hybrid check result that stores the numerical values explicitly.
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, statesWithProbability01.second.toAdd(), maybeStates, odd, x));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), statesWithProbability01.second.toAdd()));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(this->getModel().getReachableStates(), this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector())));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType> HybridMdpPrctlModelChecker<DdType, ValueType>::computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates) {
+            storm::dd::Add<DdType> result = transitionMatrix * nextStates.swapVariables(model.getRowColumnMetaVariablePairs()).toAdd();
+            return result.sumAbstract(model.getColumnVariables());
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            STORM_LOG_THROW(pathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
+            std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& leftResult = leftResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            SymbolicQualitativeCheckResult<DdType> const& rightResult = rightResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeBoundedUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
+            // probability 0 or 1 of satisfying the until-formula.
+            storm::dd::Bdd<DdType> statesWithProbabilityGreater0;
+            if (minimize) {
+                statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(model, transitionMatrix.notZero(), phiStates, psiStates);
+            } else {
+                statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(model, transitionMatrix.notZero(), phiStates, psiStates);
+            }
+            storm::dd::Bdd<DdType> maybeStates = statesWithProbabilityGreater0 && !psiStates && model.getReachableStates();
+                        
+            // If there are maybe states, we need to perform matrix-vector multiplications.
+            if (!maybeStates.isZero()) {
+                // Create the ODD for the translation between symbolic and explicit storage.
+                storm::dd::Odd<DdType> odd(maybeStates);
+                
+                // Create the matrix and the vector for the equation system.
+                storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                
+                // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                // non-maybe states in the matrix.
+                storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                
+                // Then compute the vector that contains the one-step probabilities to a state with probability 1 for all
+                // maybe states.
+                storm::dd::Add<DdType> prob1StatesAsColumn = psiStates.toAdd().swapVariables(model.getRowColumnMetaVariablePairs());
+                storm::dd::Add<DdType> subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables());
+                
+                // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                
+                // Finally cut away all columns targeting non-maybe states.
+                submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                
+                // Create the solution vector.
+                std::vector<ValueType> x(maybeStates.getNonZeroCount(), storm::utility::zero<ValueType>());
+                
+                // Translate the symbolic matrix/vector to their explicit representations.
+                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                
+                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                solver->performMatrixVectorMultiplication(minimize, x, &explicitRepresentation.second, stepBound);
+                
+                // Return a hybrid check result that stores the numerical values explicitly.
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.toAdd(), maybeStates, odd, x));
+            } else {
+                return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), psiStates.toAdd()));
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards() || model.hasTransitionRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Compute the reward vector to add in each step based on the available reward models.
+            storm::dd::Add<DdType> totalRewardVector = model.hasStateRewards() ? model.getStateRewardVector() : model.getManager().getAddZero();
+            if (model.hasTransitionRewards()) {
+                totalRewardVector += (transitionMatrix * model.getTransitionRewardMatrix()).sumAbstract(model.getColumnVariables());
+            }
+            
+            // Create the ODD for the translation between symbolic and explicit storage.
+            storm::dd::Odd<DdType> odd(model.getReachableStates());
+            
+            // Create the solution vector.
+            std::vector<ValueType> x(model.getNumberOfStates(), storm::utility::zero<ValueType>());
+            
+            // Translate the symbolic matrix/vector to their explicit representations.
+            storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
+            std::vector<ValueType> b = totalRewardVector.template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
+            solver->performMatrixVectorMultiplication(minimize, x, &b, stepBound);
+            
+            // Return a hybrid check result that stores the numerical values explicitly.
+            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
+            return this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), rewardPathFormula.getDiscreteTimeBound(), *this->linearEquationSolverFactory);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            // Only compute the result if the model has at least one reward this->getModel().
+            STORM_LOG_THROW(model.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Create the ODD for the translation between symbolic and explicit storage.
+            storm::dd::Odd<DdType> odd(model.getReachableStates());
+
+            // Translate the symbolic matrix to its explicit representations.
+            storm::storage::SparseMatrix<ValueType> explicitMatrix = transitionMatrix.toMatrix(model.getNondeterminismVariables(), odd, odd);
+
+            // Create the solution vector (and initialize it to the state rewards of the model).
+            std::vector<ValueType> x = model.getStateRewardVector().template toVector<ValueType>(model.getNondeterminismVariables(), odd, explicitMatrix.getRowGroupIndices());
+            
+            // Perform the matrix-vector multiplication.
+            std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitMatrix);
+            solver->performMatrixVectorMultiplication(minimize, x, nullptr, stepBound);
+            
+            // Return a hybrid check result that stores the numerical values explicitly.
+            return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().getAddZero(), model.getReachableStates(), odd, x));
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
+            std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
+            SymbolicQualitativeCheckResult<DdType> const& subResult = subResultPointer->asSymbolicQualitativeCheckResult<DdType>();
+            return this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel(), this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<CheckResult> HybridMdpPrctlModelChecker<DdType, ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
+            
+            // Only compute the result if there is at least one reward model.
+            STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
+            
+            // Determine which states have a reward of infinity by definition.
+            storm::dd::Bdd<DdType> infinityStates;
+            storm::dd::Bdd<DdType> transitionMatrixBdd = transitionMatrix.notZero();
+            if (minimize) {
+                infinityStates = storm::utility::graph::performProb1A(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0A(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
+            } else {
+                infinityStates = storm::utility::graph::performProb1E(model, transitionMatrixBdd, model.getReachableStates(), targetStates, storm::utility::graph::performProbGreater0E(model, transitionMatrixBdd, model.getReachableStates(), targetStates));
+            }
+            infinityStates = !infinityStates && model.getReachableStates();
+            storm::dd::Bdd<DdType> maybeStates = (!targetStates && !infinityStates) && model.getReachableStates();
+            STORM_LOG_INFO("Found " << infinityStates.getNonZeroCount() << " 'infinity' states.");
+            STORM_LOG_INFO("Found " << targetStates.getNonZeroCount() << " 'target' states.");
+            STORM_LOG_INFO("Found " << maybeStates.getNonZeroCount() << " 'maybe' states.");
+            
+            // Check whether we need to compute exact rewards for some states.
+            if (qualitative) {
+                // Set the values for all maybe-states to 1 to indicate that their reward values
+                // are neither 0 nor infinity.
+                return std::unique_ptr<CheckResult>(new SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()) + maybeStates.toAdd() * model.getManager().getConstant(storm::utility::one<ValueType>())));
+            } else {
+                // If there are maybe states, we need to solve an equation system.
+                if (!maybeStates.isZero()) {
+                    // Create the ODD for the translation between symbolic and explicit storage.
+                    storm::dd::Odd<DdType> odd(maybeStates);
+                    
+                    // Create the matrix and the vector for the equation system.
+                    storm::dd::Add<DdType> maybeStatesAdd = maybeStates.toAdd();
+                    
+                    // Start by cutting away all rows that do not belong to maybe states. Note that this leaves columns targeting
+                    // non-maybe states in the matrix.
+                    storm::dd::Add<DdType> submatrix = transitionMatrix * maybeStatesAdd;
+                    
+                    // Then compute the state reward vector to use in the computation.
+                    storm::dd::Add<DdType> subvector = stateRewardVector ? maybeStatesAdd * stateRewardVector.get() : model.getManager().getAddZero();
+                    if (transitionRewardMatrix) {
+                        subvector += (submatrix * transitionRewardMatrix.get()).sumAbstract(model.getColumnVariables());
+                    }
+                    
+                    // Before cutting the non-maybe columns, we need to compute the sizes of the row groups.
+                    std::vector<uint_fast64_t> rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).toAdd().sumAbstract(model.getNondeterminismVariables()).template toVector<uint_fast64_t>(odd);
+                    
+                    // Finally cut away all columns targeting non-maybe states.
+                    submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
+                    
+                    // Create the solution vector.
+                    std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
+                    
+                    // Translate the symbolic matrix/vector to their explicit representations.
+                    std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
+                    
+                    // Now solve the resulting equation system.
+                    std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitRepresentation.first);
+                    solver->solveEquationSystem(minimize, x, explicitRepresentation.second);
+                                        
+                    // Return a hybrid check result that stores the numerical values explicitly.
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>()), maybeStates, odd, x));
+                } else {
+                    return std::unique_ptr<CheckResult>(new storm::modelchecker::SymbolicQuantitativeCheckResult<DdType>(model.getReachableStates(), infinityStates.toAdd() * model.getManager().getConstant(storm::utility::infinity<ValueType>())));
+                }
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::models::symbolic::Mdp<DdType> const& HybridMdpPrctlModelChecker<DdType, ValueType>::getModel() const {
+            return this->template getModelAs<storm::models::symbolic::Mdp<DdType>>();
+        }
+        
+        template class HybridMdpPrctlModelChecker<storm::dd::DdType::CUDD, double>;
+    }
+}
\ No newline at end of file
diff --git a/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h
new file mode 100644
index 000000000..24a33ee79
--- /dev/null
+++ b/src/modelchecker/prctl/HybridMdpPrctlModelChecker.h
@@ -0,0 +1,44 @@
+#ifndef STORM_MODELCHECKER_HYBRIDMDPPRCTLMODELCHECKER_H_
+#define STORM_MODELCHECKER_HYBRIDMDPPRCTLMODELCHECKER_H_
+
+#include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
+#include "src/models/symbolic/Mdp.h"
+#include "src/utility/solver.h"
+
+namespace storm {
+    namespace modelchecker {
+        template<storm::dd::DdType DdType, typename ValueType>
+        class HybridMdpPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType> {
+        public:
+            explicit HybridMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model);
+            explicit HybridMdpPrctlModelChecker(storm::models::symbolic::Mdp<DdType> const& model, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            
+            // The implemented methods of the AbstractModelChecker interface.
+            virtual bool canHandle(storm::logic::Formula const& formula) const override;
+            virtual std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            virtual std::unique_ptr<CheckResult> computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative = false, boost::optional<storm::logic::OptimalityType> const& optimalityType = boost::optional<storm::logic::OptimalityType>()) override;
+            
+        protected:
+            storm::models::symbolic::Mdp<DdType> const& getModel() const override;
+            
+        private:
+            // The methods that perform the actual checking.
+            static std::unique_ptr<CheckResult> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static storm::dd::Add<DdType> computeNextProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
+            static std::unique_ptr<CheckResult> computeUntilProbabilitiesHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeCumulativeRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeInstantaneousRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, uint_fast64_t stepBound, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static std::unique_ptr<CheckResult> computeReachabilityRewardsHelper(bool minimize, storm::models::symbolic::NondeterministicModel<DdType> const& model, storm::dd::Add<DdType> const& transitionMatrix, boost::optional<storm::dd::Add<DdType>> const& stateRewardVector, boost::optional<storm::dd::Add<DdType>> const& transitionRewardMatrix, storm::dd::Bdd<DdType> const& targetStates, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
+            
+            // An object that is used for retrieving linear equation solvers.
+            std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+        };
+        
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_HYBRIDMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index 8f63110c8..f16951d12 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -251,7 +251,7 @@ namespace storm {
                 // Converting the matrix from the fixpoint notation to the form needed for the equation
                 // system. That is, we go from x = A*x + b to (I-A)x = b.
                 submatrix.convertToEquationSystem();
-                
+                                
                 // Initialize the x vector with 1 for each element. This is the initial guess for
                 // the iterative solvers.
                 std::vector<ValueType> x(submatrix.getColumnCount(), storm::utility::one<ValueType>());
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 0a9007c22..518d21972 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -72,7 +72,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
@@ -96,7 +96,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, subResult.getTruthValuesVector())));
@@ -166,7 +166,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilities(storm::logic::UntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
@@ -206,7 +206,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeCumulativeRewards(storm::logic::CumulativeRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeCumulativeRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound())));
         }
@@ -228,7 +228,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeInstantaneousRewards(storm::logic::InstantaneousRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             STORM_LOG_THROW(rewardPathFormula.hasDiscreteTimeBound(), storm::exceptions::InvalidArgumentException, "Formula needs to have a discrete time bound.");
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeInstantaneousRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, rewardPathFormula.getDiscreteTimeBound())));
         }
@@ -303,6 +303,7 @@ namespace storm {
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = MinMaxLinearEquationSolverFactory.create(submatrix);
                 solver->solveEquationSystem(minimize, x, b);
                 
+                
                 // Set values of resulting vector according to result.
                 storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
             }
@@ -316,7 +317,7 @@ namespace storm {
         
         template<typename ValueType>
         std::unique_ptr<CheckResult> SparseMdpPrctlModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
+            STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->MinMaxLinearEquationSolverFactory, qualitative)));
diff --git a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp
index 1552235d9..4fa0dbafe 100644
--- a/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp
+++ b/src/modelchecker/propositional/SymbolicPropositionalModelChecker.cpp
@@ -1,6 +1,7 @@
 #include "src/modelchecker/propositional/SymbolicPropositionalModelChecker.h"
 
 #include "src/models/symbolic/Dtmc.h"
+#include "src/models/symbolic/Ctmc.h"
 #include "src/models/symbolic/Mdp.h"
 
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
@@ -48,6 +49,7 @@ namespace storm {
         
         // Explicitly instantiate the template class.
         template storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD> const& SymbolicPropositionalModelChecker<storm::dd::DdType::CUDD>::getModelAs() const;
+        template storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD> const& SymbolicPropositionalModelChecker<storm::dd::DdType::CUDD>::getModelAs() const;
         template storm::models::symbolic::Mdp<storm::dd::DdType::CUDD> const& SymbolicPropositionalModelChecker<storm::dd::DdType::CUDD>::getModelAs() const;
         template class SymbolicPropositionalModelChecker<storm::dd::DdType::CUDD>;
     }
diff --git a/src/modelchecker/results/HybridQuantitativeCheckResult.cpp b/src/modelchecker/results/HybridQuantitativeCheckResult.cpp
index 8132b75b6..abf56e163 100644
--- a/src/modelchecker/results/HybridQuantitativeCheckResult.cpp
+++ b/src/modelchecker/results/HybridQuantitativeCheckResult.cpp
@@ -1,5 +1,6 @@
 #include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/storage/dd/CuddDdManager.h"
 
 #include "src/exceptions/InvalidOperationException.h"
@@ -32,6 +33,16 @@ namespace storm {
             return std::unique_ptr<SymbolicQualitativeCheckResult<Type>>(new SymbolicQualitativeCheckResult<Type>(reachableStates, symbolicResult));
         }
         
+        template<storm::dd::DdType Type>
+        std::unique_ptr<CheckResult> HybridQuantitativeCheckResult<Type>::toExplicitQuantitativeCheckResult() const {
+            storm::dd::Bdd<Type> allStates = symbolicStates || explicitStates;
+            storm::dd::Odd<Type> allStatesOdd(allStates);
+            
+            std::vector<double> fullExplicitValues = symbolicValues.template toVector<double>(allStatesOdd);
+            this->odd.expandExplicitVector(allStatesOdd, this->explicitValues, fullExplicitValues);
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<double>(std::move(fullExplicitValues)));
+        }
+        
         template<storm::dd::DdType Type>
         bool HybridQuantitativeCheckResult<Type>::isHybrid() const {
             return true;
@@ -139,7 +150,7 @@ namespace storm {
         
         template<storm::dd::DdType Type>
         double HybridQuantitativeCheckResult<Type>::getMax() const {
-            double max = this->symbolicValues.getMin();
+            double max = this->symbolicValues.getMax();
             if (!explicitStates.isZero()) {
                 for (auto const& element : explicitValues) {
                     max = std::max(max, element);
diff --git a/src/modelchecker/results/HybridQuantitativeCheckResult.h b/src/modelchecker/results/HybridQuantitativeCheckResult.h
index 160bbbb6d..05033f78f 100644
--- a/src/modelchecker/results/HybridQuantitativeCheckResult.h
+++ b/src/modelchecker/results/HybridQuantitativeCheckResult.h
@@ -25,6 +25,8 @@ namespace storm {
             
             virtual std::unique_ptr<CheckResult> compareAgainstBound(storm::logic::ComparisonType comparisonType, double bound) const override;
             
+            std::unique_ptr<CheckResult> toExplicitQuantitativeCheckResult() const;
+            
             virtual bool isHybrid() const override;
             virtual bool isResultForAllStates() const override;
             
diff --git a/src/models/symbolic/Ctmc.cpp b/src/models/symbolic/Ctmc.cpp
index b90aa5afd..01f4cb4d7 100644
--- a/src/models/symbolic/Ctmc.cpp
+++ b/src/models/symbolic/Ctmc.cpp
@@ -18,7 +18,12 @@ namespace storm {
                              boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector,
                              boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix)
             : DeterministicModel<Type>(storm::models::ModelType::Ctmc, manager, reachableStates, initialStates, transitionMatrix, rowVariables, rowExpressionAdapter, columnVariables, columnExpressionAdapter, rowColumnMetaVariablePairs, labelToExpressionMap, optionalStateRewardVector, optionalTransitionRewardMatrix) {
-                // Intentionally left empty.
+                exitRates = this->getTransitionMatrix().sumAbstract(this->getColumnVariables());
+            }
+            
+            template<storm::dd::DdType Type>
+            storm::dd::Add<Type> const& Ctmc<Type>::getExitRateVector() const {
+                return exitRates;
             }
             
             // Explicitly instantiate the template class.
diff --git a/src/models/symbolic/Ctmc.h b/src/models/symbolic/Ctmc.h
index 7d46dcd73..e2f975f76 100644
--- a/src/models/symbolic/Ctmc.h
+++ b/src/models/symbolic/Ctmc.h
@@ -52,6 +52,16 @@ namespace storm {
                      std::map<std::string, storm::expressions::Expression> labelToExpressionMap = std::map<std::string, storm::expressions::Expression>(),
                      boost::optional<storm::dd::Add<Type>> const& optionalStateRewardVector = boost::optional<storm::dd::Dd<Type>>(),
                      boost::optional<storm::dd::Add<Type>> const& optionalTransitionRewardMatrix = boost::optional<storm::dd::Dd<Type>>());
+                
+                /*!
+                 * Retrieves the exit rate vector of the CTMC.
+                 *
+                 * @return The exit rate vector.
+                 */
+                storm::dd::Add<Type> const& getExitRateVector() const;
+                
+            private:
+                storm::dd::Add<Type> exitRates;
             };
             
         } // namespace symbolic
diff --git a/src/models/symbolic/Model.cpp b/src/models/symbolic/Model.cpp
index 0fe4f0d13..8b72f3687 100644
--- a/src/models/symbolic/Model.cpp
+++ b/src/models/symbolic/Model.cpp
@@ -87,11 +87,21 @@ namespace storm {
                 return transitionRewardMatrix.get();
             }
             
+            template<storm::dd::DdType Type>
+            boost::optional<storm::dd::Add<Type>> const& Model<Type>::getOptionalTransitionRewardMatrix() const {
+                return transitionRewardMatrix;
+            }
+            
             template<storm::dd::DdType Type>
             storm::dd::Add<Type> const& Model<Type>::getStateRewardVector() const {
                 return stateRewardVector.get();
             }
             
+            template<storm::dd::DdType Type>
+            boost::optional<storm::dd::Add<Type>> const& Model<Type>::getOptionalStateRewardVector() const {
+                return stateRewardVector;
+            }
+            
             template<storm::dd::DdType Type>
             bool Model<Type>::hasStateRewards() const {
                 return static_cast<bool>(stateRewardVector);
diff --git a/src/models/symbolic/Model.h b/src/models/symbolic/Model.h
index 211a4cc66..23d2030f3 100644
--- a/src/models/symbolic/Model.h
+++ b/src/models/symbolic/Model.h
@@ -134,6 +134,13 @@ namespace storm {
                  */
                 storm::dd::Add<Type>& getTransitionMatrix();
                 
+                /*!
+                 * Retrieves the (optional) matrix representing the transition rewards of the model.
+                 *
+                 * @return The matrix representing the transition rewards of the model.
+                 */
+                boost::optional<storm::dd::Add<Type>> const& getOptionalTransitionRewardMatrix() const;
+                
                 /*!
                  * Retrieves the matrix representing the transition rewards of the model. Note that calling this method
                  * is only valid if the model has transition rewards.
@@ -157,6 +164,13 @@ namespace storm {
                  * @return A vector representing the state rewards of the model.
                  */
                 storm::dd::Add<Type> const& getStateRewardVector() const;
+
+                /*!
+                 * Retrieves an (optional) vector representing the state rewards of the model.
+                 *
+                 * @return A vector representing the state rewards of the model.
+                 */
+                boost::optional<storm::dd::Add<Type>> const& getOptionalStateRewardVector() const;
                 
                 /*!
                  * Retrieves whether this model has state rewards.
diff --git a/src/settings/Argument.h b/src/settings/Argument.h
index c697fdb45..7aad2c048 100644
--- a/src/settings/Argument.h
+++ b/src/settings/Argument.h
@@ -73,12 +73,12 @@ namespace storm {
 				return this->setFromTypeValue(newValue);
 			}
             
-			bool setFromTypeValue(T const& newValue) {
+			bool setFromTypeValue(T const& newValue, bool hasBeenSet = true) {
 				if (!this->validate(newValue)) {
                     return false;
 				}
 				this->argumentValue = newValue;
-				this->hasBeenSet = true;
+				this->hasBeenSet = hasBeenSet;
 				return true;
 			}
             
@@ -120,7 +120,7 @@ namespace storm {
             
 			void setFromDefaultValue() override {
                 STORM_LOG_THROW(this->hasDefaultValue, storm::exceptions::IllegalFunctionCallException, "Unable to set value from default value, because the argument has none.");
-				bool result = this->setFromTypeValue(this->defaultValue);
+				bool result = this->setFromTypeValue(this->defaultValue, false);
                 STORM_LOG_THROW(result, storm::exceptions::IllegalArgumentValueException, "Unable to assign default value to argument, because it was rejected.");
 			}
             
diff --git a/src/storage/dd/CuddAdd.cpp b/src/storage/dd/CuddAdd.cpp
index 8152083cd..6b7f2067f 100644
--- a/src/storage/dd/CuddAdd.cpp
+++ b/src/storage/dd/CuddAdd.cpp
@@ -424,7 +424,47 @@ namespace storm {
         std::vector<ValueType> Add<DdType::CUDD>::toVector(Odd<DdType::CUDD> const& rowOdd) const {
             std::vector<ValueType> result(rowOdd.getTotalOffset());
             std::vector<uint_fast64_t> ddVariableIndices = this->getSortedVariableIndices();
-            addToVectorRec(this->getCuddDdNode(), 0, ddVariableIndices.size(), 0, rowOdd, ddVariableIndices, result);
+            addToVector(rowOdd, ddVariableIndices, result);
+            return result;
+        }
+        
+        template<typename ValueType>
+        std::vector<ValueType> Add<DdType::CUDD>::toVector(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, std::vector<uint_fast64_t> groupOffsets) const {
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            
+            // Prepare the proper sets of meta variables.
+            for (auto const& variable : this->getContainedMetaVariables()) {
+                if (groupMetaVariables.find(variable) != groupMetaVariables.end()) {
+                    continue;
+                }
+                
+                rowMetaVariables.insert(variable);
+            }
+            std::vector<uint_fast64_t> ddGroupVariableIndices;
+            for (auto const& variable : groupMetaVariables) {
+                DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddGroupVariableIndices.push_back(ddVariable.getIndex());
+                }
+            }
+            std::vector<uint_fast64_t> ddRowVariableIndices;
+            for (auto const& variable : rowMetaVariables) {
+                DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddRowVariableIndices.push_back(ddVariable.getIndex());
+                }
+            }
+            
+            // Start by splitting the symbolic vector into groups.
+            std::vector<Add<DdType::CUDD>> groups;
+            splitGroupsRec(this->getCuddDdNode(), groups, ddGroupVariableIndices, 0, ddGroupVariableIndices.size(), rowMetaVariables);
+            
+            // Now iterate over the groups and add them to the resulting vector.
+            std::vector<ValueType> result(groupOffsets.back(), storm::utility::zero<ValueType>());
+            for (uint_fast64_t i = 0; i < groups.size(); ++i) {
+                toVectorRec(groups[i].getCuddDdNode(), result, groupOffsets, rowOdd, 0, ddRowVariableIndices.size(), 0, ddRowVariableIndices);
+            }
+            
             return result;
         }
         
@@ -520,11 +560,33 @@ namespace storm {
             return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(trivialRowGroupIndices));
         }
         
+        storm::storage::SparseMatrix<double> Add<DdType::CUDD>::toMatrix(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            for (auto const& variable : this->getContainedMetaVariables()) {
+                // If the meta variable is a group meta variable, we do not insert it into the set of row/column meta variables.
+                if (groupMetaVariables.find(variable) != groupMetaVariables.end()) {
+                    continue;
+                }
+                
+                if (variable.getName().size() > 0 && variable.getName().back() == '\'') {
+                    columnMetaVariables.insert(variable);
+                } else {
+                    rowMetaVariables.insert(variable);
+                }
+            }
+            
+            // Create the canonical row group sizes and build the matrix.
+            return toMatrix(rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd);
+        }
+        
         storm::storage::SparseMatrix<double> Add<DdType::CUDD>::toMatrix(std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
             std::vector<uint_fast64_t> ddRowVariableIndices;
             std::vector<uint_fast64_t> ddColumnVariableIndices;
             std::vector<uint_fast64_t> ddGroupVariableIndices;
             std::set<storm::expressions::Variable> rowAndColumnMetaVariables;
+            boost::optional<std::vector<double>> optionalExplicitVector;
             
             for (auto const& variable : rowMetaVariables) {
                 DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
@@ -550,8 +612,6 @@ namespace storm {
             }
             std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end());
             
-            // TODO: assert that the group variables are at the very top of the variable ordering?
-            
             // Start by computing the offsets (in terms of rows) for each row group.
             Add<DdType::CUDD> stateToNumberOfChoices = this->notZero().existsAbstract(columnMetaVariables).toAdd().sumAbstract(groupMetaVariables);
             std::vector<uint_fast64_t> rowGroupIndices = stateToNumberOfChoices.toVector<uint_fast64_t>(rowOdd);
@@ -575,21 +635,22 @@ namespace storm {
             
             // Now compute the indices at which the individual rows start.
             std::vector<uint_fast64_t> rowIndications(rowGroupIndices.back() + 1);
+            
             std::vector<storm::dd::Add<DdType::CUDD>> statesWithGroupEnabled(groups.size());
+            storm::dd::Add<storm::dd::DdType::CUDD> stateToRowGroupCount = this->getDdManager()->getAddZero();
             for (uint_fast64_t i = 0; i < groups.size(); ++i) {
                 auto const& dd = groups[i];
                 
                 toMatrixRec(dd.getCuddDdNode(), rowIndications, columnsAndValues, rowGroupIndices, rowOdd, columnOdd, 0, 0, ddRowVariableIndices.size() + ddColumnVariableIndices.size(), 0, 0, ddRowVariableIndices, ddColumnVariableIndices, false);
                 
                 statesWithGroupEnabled[i] = dd.notZero().existsAbstract(columnMetaVariables).toAdd();
-                addToVectorRec(statesWithGroupEnabled[i].getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices);
+                stateToRowGroupCount += statesWithGroupEnabled[i];
+                statesWithGroupEnabled[i].addToVector(rowOdd, ddRowVariableIndices, rowGroupIndices);
             }
             
             // Since we modified the rowGroupIndices, we need to restore the correct values.
-            for (uint_fast64_t i = rowGroupIndices.size() - 1; i > 0; --i) {
-                rowGroupIndices[i] = rowGroupIndices[i - 1];
-            }
-            rowGroupIndices[0] = 0;
+            std::function<uint_fast64_t (uint_fast64_t const&, double const&)> fct = [] (uint_fast64_t const& a, double const& b) -> uint_fast64_t { return a - static_cast<uint_fast64_t>(b); };
+            modifyVectorRec(stateToRowGroupCount.getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices, fct);
             
             // Now that we computed the number of entries in each row, compute the corresponding offsets in the entry vector.
             tmp = 0;
@@ -607,22 +668,162 @@ namespace storm {
                 
                 toMatrixRec(dd.getCuddDdNode(), rowIndications, columnsAndValues, rowGroupIndices, rowOdd, columnOdd, 0, 0, ddRowVariableIndices.size() + ddColumnVariableIndices.size(), 0, 0, ddRowVariableIndices, ddColumnVariableIndices, true);
                 
-                addToVectorRec(statesWithGroupEnabled[i].getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices);
+                statesWithGroupEnabled[i].addToVector(rowOdd, ddRowVariableIndices, rowGroupIndices);
             }
             
             // Since we modified the rowGroupIndices, we need to restore the correct values.
-            for (uint_fast64_t i = rowGroupIndices.size() - 1; i > 0; --i) {
-                rowGroupIndices[i] = rowGroupIndices[i - 1];
+            modifyVectorRec(stateToRowGroupCount.getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices, fct);
+            
+            // Since the last call to toMatrixRec modified the rowIndications, we need to restore the correct values.
+            for (uint_fast64_t i = rowIndications.size() - 1; i > 0; --i) {
+                rowIndications[i] = rowIndications[i - 1];
+            }
+            rowIndications[0] = 0;
+            
+            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices));
+        }
+        
+        std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> Add<DdType::CUDD>::toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
+            std::set<storm::expressions::Variable> rowMetaVariables;
+            std::set<storm::expressions::Variable> columnMetaVariables;
+            
+            for (auto const& variable : this->getContainedMetaVariables()) {
+                // If the meta variable is a group meta variable, we do not insert it into the set of row/column meta variables.
+                if (groupMetaVariables.find(variable) != groupMetaVariables.end()) {
+                    continue;
+                }
+                
+                if (variable.getName().size() > 0 && variable.getName().back() == '\'') {
+                    columnMetaVariables.insert(variable);
+                } else {
+                    rowMetaVariables.insert(variable);
+                }
+            }
+            
+            // Create the canonical row group sizes and build the matrix.
+            return toMatrixVector(vector, std::move(rowGroupSizes), rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd);
+        }
+        
+        std::pair<storm::storage::SparseMatrix<double>,std::vector<double>> Add<DdType::CUDD>::toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupIndices, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const {
+            std::vector<uint_fast64_t> ddRowVariableIndices;
+            std::vector<uint_fast64_t> ddColumnVariableIndices;
+            std::vector<uint_fast64_t> ddGroupVariableIndices;
+            std::set<storm::expressions::Variable> rowAndColumnMetaVariables;
+
+            for (auto const& variable : rowMetaVariables) {
+                DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddRowVariableIndices.push_back(ddVariable.getIndex());
+                }
+                rowAndColumnMetaVariables.insert(variable);
+            }
+            std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
+            for (auto const& variable : columnMetaVariables) {
+                DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddColumnVariableIndices.push_back(ddVariable.getIndex());
+                }
+                rowAndColumnMetaVariables.insert(variable);
+            }
+            std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
+            for (auto const& variable : groupMetaVariables) {
+                DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
+                for (auto const& ddVariable : metaVariable.getDdVariables()) {
+                    ddGroupVariableIndices.push_back(ddVariable.getIndex());
+                }
+            }
+            std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end());
+            
+            // Transform the row group sizes to the actual row group indices.
+            rowGroupIndices.resize(rowGroupIndices.size() + 1);
+            uint_fast64_t tmp = 0;
+            uint_fast64_t tmp2 = 0;
+            for (uint_fast64_t i = 1; i < rowGroupIndices.size(); ++i) {
+                tmp2 = rowGroupIndices[i];
+                rowGroupIndices[i] = rowGroupIndices[i - 1] + tmp;
+                std::swap(tmp, tmp2);
             }
             rowGroupIndices[0] = 0;
             
+            // Create the explicit vector we need to fill later.
+            std::vector<double> explicitVector(rowGroupIndices.back());
+            
+            // Next, we split the matrix into one for each group. This only works if the group variables are at the very
+            // top.
+            std::vector<std::pair<Add<DdType::CUDD>, Add<DdType::CUDD>>> groups;
+            splitGroupsRec(this->getCuddDdNode(), vector.getCuddDdNode(), groups, ddGroupVariableIndices, 0, ddGroupVariableIndices.size(), rowAndColumnMetaVariables, rowMetaVariables);
+            
+            // Create the actual storage for the non-zero entries.
+            std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>> columnsAndValues(this->getNonZeroCount());
+            
+            // Now compute the indices at which the individual rows start.
+            std::vector<uint_fast64_t> rowIndications(rowGroupIndices.back() + 1);
+
+            std::vector<storm::dd::Add<DdType::CUDD>> statesWithGroupEnabled(groups.size());
+            storm::dd::Add<storm::dd::DdType::CUDD> stateToRowGroupCount = this->getDdManager()->getAddZero();
+            for (uint_fast64_t i = 0; i < groups.size(); ++i) {
+                std::pair<storm::dd::Add<DdType::CUDD>, storm::dd::Add<DdType::CUDD>> ddPair = groups[i];
+                
+                toMatrixRec(ddPair.first.getCuddDdNode(), rowIndications, columnsAndValues, rowGroupIndices, rowOdd, columnOdd, 0, 0, ddRowVariableIndices.size() + ddColumnVariableIndices.size(), 0, 0, ddRowVariableIndices, ddColumnVariableIndices, false);
+                toVectorRec(ddPair.second.getCuddDdNode(), explicitVector, rowGroupIndices, rowOdd, 0, ddRowVariableIndices.size(), 0, ddRowVariableIndices);
+                
+                statesWithGroupEnabled[i] = (ddPair.first.notZero().existsAbstract(columnMetaVariables) || ddPair.second.notZero()).toAdd();
+                stateToRowGroupCount += statesWithGroupEnabled[i];
+                statesWithGroupEnabled[i].addToVector(rowOdd, ddRowVariableIndices, rowGroupIndices);
+            }
+            
+            // Since we modified the rowGroupIndices, we need to restore the correct values.
+            std::function<uint_fast64_t (uint_fast64_t const&, double const&)> fct = [] (uint_fast64_t const& a, double const& b) -> uint_fast64_t { return a - static_cast<uint_fast64_t>(b); };
+            modifyVectorRec(stateToRowGroupCount.getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices, fct);
+            
+            // Now that we computed the number of entries in each row, compute the corresponding offsets in the entry vector.
+            tmp = 0;
+            tmp2 = 0;
+            for (uint_fast64_t i = 1; i < rowIndications.size(); ++i) {
+                tmp2 = rowIndications[i];
+                rowIndications[i] = rowIndications[i - 1] + tmp;
+                std::swap(tmp, tmp2);
+            }
+            rowIndications[0] = 0;
+            
+            // Now actually fill the entry vector.
+            for (uint_fast64_t i = 0; i < groups.size(); ++i) {
+                auto const& dd = groups[i].first;
+                
+                toMatrixRec(dd.getCuddDdNode(), rowIndications, columnsAndValues, rowGroupIndices, rowOdd, columnOdd, 0, 0, ddRowVariableIndices.size() + ddColumnVariableIndices.size(), 0, 0, ddRowVariableIndices, ddColumnVariableIndices, true);
+                
+                statesWithGroupEnabled[i].addToVector(rowOdd, ddRowVariableIndices, rowGroupIndices);
+            }
+            
+            // Since we modified the rowGroupIndices, we need to restore the correct values.
+            modifyVectorRec(stateToRowGroupCount.getCuddDdNode(), 0, ddRowVariableIndices.size(), 0, rowOdd, ddRowVariableIndices, rowGroupIndices, fct);
+            
             // Since the last call to toMatrixRec modified the rowIndications, we need to restore the correct values.
             for (uint_fast64_t i = rowIndications.size() - 1; i > 0; --i) {
                 rowIndications[i] = rowIndications[i - 1];
             }
             rowIndications[0] = 0;
             
-            return storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices));
+            return std::make_pair(storm::storage::SparseMatrix<double>(columnOdd.getTotalOffset(), std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices)), std::move(explicitVector));
+        }
+        
+        template<typename ValueType>
+        void Add<DdType::CUDD>::toVectorRec(DdNode const* dd, std::vector<ValueType>& result, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const {
+            // For the empty DD, we do not need to add any entries.
+            if (dd == Cudd_ReadZero(this->getDdManager()->getCuddManager().getManager())) {
+                return;
+            }
+
+            // If we are at the maximal level, the value to be set is stored as a constant in the DD.
+            if (currentRowLevel == maxLevel) {
+                result[rowGroupOffsets[currentRowOffset]] = Cudd_V(dd);
+            } else if (ddRowVariableIndices[currentRowLevel] < dd->index) {
+                toVectorRec(dd, result, rowGroupOffsets, rowOdd.getElseSuccessor(), currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices);
+                toVectorRec(dd, result, rowGroupOffsets, rowOdd.getThenSuccessor(), currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices);
+            } else {
+                toVectorRec(Cudd_E(dd), result, rowGroupOffsets, rowOdd.getElseSuccessor(), currentRowLevel + 1, maxLevel, currentRowOffset, ddRowVariableIndices);
+                toVectorRec(Cudd_T(dd), result, rowGroupOffsets, rowOdd.getThenSuccessor(), currentRowLevel + 1, maxLevel, currentRowOffset + rowOdd.getElseOffset(), ddRowVariableIndices);
+            }
         }
         
         void Add<DdType::CUDD>::toMatrixRec(DdNode const* dd, std::vector<uint_fast64_t>& rowIndications, std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>>& columnsAndValues, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, Odd<DdType::CUDD> const& columnOdd, uint_fast64_t currentRowLevel, uint_fast64_t currentColumnLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, uint_fast64_t currentColumnOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices, std::vector<uint_fast64_t> const& ddColumnVariableIndices, bool generateValues) const {
@@ -694,8 +895,39 @@ namespace storm {
             }
         }
         
+        void Add<DdType::CUDD>::splitGroupsRec(DdNode* dd1, DdNode* dd2, std::vector<std::pair<Add<DdType::CUDD>, Add<DdType::CUDD>>>& groups, std::vector<uint_fast64_t> const& ddGroupVariableIndices, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::set<storm::expressions::Variable> const& remainingMetaVariables1, std::set<storm::expressions::Variable> const& remainingMetaVariables2) const {
+            // For the empty DD, we do not need to create a group.
+            if (dd1 == Cudd_ReadZero(this->getDdManager()->getCuddManager().getManager()) && dd2 == Cudd_ReadZero(this->getDdManager()->getCuddManager().getManager())) {
+                return;
+            }
+            
+            if (currentLevel == maxLevel) {
+                groups.push_back(std::make_pair(Add<DdType::CUDD>(this->getDdManager(), ADD(this->getDdManager()->getCuddManager(), dd1), remainingMetaVariables1), Add<DdType::CUDD>(this->getDdManager(), ADD(this->getDdManager()->getCuddManager(), dd2), remainingMetaVariables2)));
+            } else if (ddGroupVariableIndices[currentLevel] < dd1->index) {
+                if (ddGroupVariableIndices[currentLevel] < dd2->index) {
+                    splitGroupsRec(dd1, dd2, groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                    splitGroupsRec(dd1, dd2, groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                } else {
+                    splitGroupsRec(dd1, Cudd_T(dd2), groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                    splitGroupsRec(dd1, Cudd_E(dd2), groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                }
+            } else if (ddGroupVariableIndices[currentLevel] < dd2->index) {
+                splitGroupsRec(Cudd_T(dd1), dd2, groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                splitGroupsRec(Cudd_E(dd1), dd2, groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+            } else {
+                splitGroupsRec(Cudd_T(dd1), Cudd_T(dd2), groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+                splitGroupsRec(Cudd_E(dd1), Cudd_E(dd2), groups, ddGroupVariableIndices, currentLevel + 1, maxLevel, remainingMetaVariables1, remainingMetaVariables2);
+            }
+        }
+        
         template<typename ValueType>
-        void Add<DdType::CUDD>::addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const {
+        void Add<DdType::CUDD>::addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const {
+            std::function<ValueType (ValueType const&, double const&)> fct = [] (ValueType const& a, double const& b) -> ValueType { return a + b; };
+            modifyVectorRec(this->getCuddDdNode(), 0, ddVariableIndices.size(), 0, odd, ddVariableIndices, targetVector, fct);
+        }
+        
+        template<typename ValueType>
+        void Add<DdType::CUDD>::modifyVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector, std::function<ValueType (ValueType const&, double const&)> const& function) const {
             // For the empty DD, we do not need to add any entries.
             if (dd == Cudd_ReadZero(this->getDdManager()->getCuddManager().getManager())) {
                 return;
@@ -703,16 +935,16 @@ namespace storm {
             
             // If we are at the maximal level, the value to be set is stored as a constant in the DD.
             if (currentLevel == maxLevel) {
-                targetVector[currentOffset] += static_cast<ValueType>(Cudd_V(dd));
+                targetVector[currentOffset] = function(targetVector[currentOffset], Cudd_V(dd));
             } else if (ddVariableIndices[currentLevel] < dd->index) {
                 // If we skipped a level, we need to enumerate the explicit entries for the case in which the bit is set
                 // and for the one in which it is not set.
-                addToVectorRec(dd, currentLevel + 1, maxLevel, currentOffset, odd.getElseSuccessor(), ddVariableIndices, targetVector);
-                addToVectorRec(dd, currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector);
+                modifyVectorRec(dd, currentLevel + 1, maxLevel, currentOffset, odd.getElseSuccessor(), ddVariableIndices, targetVector, function);
+                modifyVectorRec(dd, currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector, function);
             } else {
                 // Otherwise, we simply recursively call the function for both (different) cases.
-                addToVectorRec(Cudd_E(dd), currentLevel + 1, maxLevel, currentOffset, odd.getElseSuccessor(), ddVariableIndices, targetVector);
-                addToVectorRec(Cudd_T(dd), currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector);
+                modifyVectorRec(Cudd_E(dd), currentLevel + 1, maxLevel, currentOffset, odd.getElseSuccessor(), ddVariableIndices, targetVector, function);
+                modifyVectorRec(Cudd_T(dd), currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector, function);
             }
         }
         
@@ -830,9 +1062,14 @@ namespace storm {
         // Explicitly instantiate some templated functions.
         template std::vector<double> Add<DdType::CUDD>::toVector() const;
         template std::vector<double> Add<DdType::CUDD>::toVector(Odd<DdType::CUDD> const& rowOdd) const;
-        template void Add<DdType::CUDD>::addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<double>& targetVector) const;
+        template void Add<DdType::CUDD>::addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<double>& targetVector) const;
+        template void Add<DdType::CUDD>::modifyVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<double>& targetVector, std::function<double (double const&, double const&)> const& function) const;
+        template std::vector<double> Add<DdType::CUDD>::toVector(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, std::vector<uint_fast64_t> groupOffsets) const;
+        template void Add<DdType::CUDD>::toVectorRec(DdNode const* dd, std::vector<double>& result, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const;
         template std::vector<uint_fast64_t> Add<DdType::CUDD>::toVector() const;
         template std::vector<uint_fast64_t> Add<DdType::CUDD>::toVector(Odd<DdType::CUDD> const& rowOdd) const;
-        template void Add<DdType::CUDD>::addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<uint_fast64_t>& targetVector) const;
+        template void Add<DdType::CUDD>::addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<uint_fast64_t>& targetVector) const;
+        template void Add<DdType::CUDD>::modifyVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<uint_fast64_t>& targetVector, std::function<uint_fast64_t (uint_fast64_t const&, double const&)> const& function) const;
+
     }
 }
\ No newline at end of file
diff --git a/src/storage/dd/CuddAdd.h b/src/storage/dd/CuddAdd.h
index 8d01839d6..01a57b629 100644
--- a/src/storage/dd/CuddAdd.h
+++ b/src/storage/dd/CuddAdd.h
@@ -1,6 +1,8 @@
 #ifndef STORM_STORAGE_DD_CUDDADD_H_
 #define STORM_STORAGE_DD_CUDDADD_H_
 
+#include <boost/optional.hpp>
+
 #include "src/storage/dd/Add.h"
 #include "src/storage/dd/CuddDd.h"
 #include "src/storage/dd/CuddDdForwardIterator.h"
@@ -512,7 +514,7 @@ namespace storm {
             /*!
              * Converts the ADD to a vector.
              *
-             * @return The double vector that is represented by this ADD.
+             * @return The vector that is represented by this ADD.
              */
             template<typename ValueType>
             std::vector<ValueType> toVector() const;
@@ -522,11 +524,23 @@ namespace storm {
              * each entry.
              *
              * @param rowOdd The ODD used for determining the correct row.
-             * @return The double vector that is represented by this ADD.
+             * @return The vector that is represented by this ADD.
              */
             template<typename ValueType>
             std::vector<ValueType> toVector(storm::dd::Odd<DdType::CUDD> const& rowOdd) const;
             
+            /*!
+             * Converts the ADD to a row-grouped vector. The given offset-labeled DD is used to determine the correct
+             * row group of each entry. Note that the group meta variables are assumed to be at the very top in the
+             * variable ordering.
+             *
+             * @param groupMetaVariables The meta variables responsible for the row-grouping.
+             * @param rowOdd The ODD used for determining the correct row.
+             * @return The vector that is represented by this ADD.
+             */
+            template<typename ValueType>
+            std::vector<ValueType> toVector(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, std::vector<uint_fast64_t> groupOffsets) const;
+            
             /*!
              * Converts the ADD to a (sparse) double matrix. All contained non-primed variables are assumed to encode the
              * row, whereas all primed variables are assumed to encode the column.
@@ -563,14 +577,27 @@ namespace storm {
              * determine the correct row and column, respectively, for each entry. Note: this function assumes that
              * the meta variables used to distinguish different row groups are at the very top of the ADD.
              *
-             * @param rowMetaVariables The meta variables that encode the rows of the matrix.
-             * @param columnMetaVariables The meta variables that encode the columns of the matrix.
              * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
              * @param rowOdd The ODD used for determining the correct row.
              * @param columnOdd The ODD used for determining the correct column.
              * @return The matrix that is represented by this ADD.
              */
-            storm::storage::SparseMatrix<double> toMatrix(std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const;
+            storm::storage::SparseMatrix<double> toMatrix(std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const;
+                        
+            /*!
+             * Converts the ADD to a row-grouped (sparse) double matrix and the given vector to a row-grouped vector.
+             * The given offset-labeled DDs are used to determine the correct row and column, respectively, for each
+             * entry. Note: this function assumes that the meta variables used to distinguish different row groups are
+             * at the very top of the ADD.
+             *
+             * @param vector The symbolic vector to convert.
+             * @param rowGroupSizes A vector specifying the sizes of the row groups.
+             * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
+             * @param rowOdd The ODD used for determining the correct row.
+             * @param columnOdd The ODD used for determining the correct column.
+             * @return The matrix that is represented by this ADD.
+             */
+            std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const;
             
             /*!
              * Exports the DD to the given file in the dot format.
@@ -632,6 +659,40 @@ namespace storm {
              */
             Add(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, ADD cuddAdd, std::set<storm::expressions::Variable> const& containedMetaVariables = std::set<storm::expressions::Variable>());
             
+            /*!
+             * Converts the ADD to a row-grouped (sparse) double matrix. If the optional vector is given, it is also
+             * translated to an explicit row-grouped vector with the same row-grouping. The given offset-labeled DDs
+             * are used to determine the correct row and column, respectively, for each entry. Note: this function
+             * assumes that the meta variables used to distinguish different row groups are at the very top of the ADD.
+             *
+             * @param rowMetaVariables The meta variables that encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables that encode the columns of the matrix.
+             * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
+             * @param rowOdd The ODD used for determining the correct row.
+             * @param columnOdd The ODD used for determining the correct column.
+             * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector
+             * (if it was given).
+             */
+            storm::storage::SparseMatrix<double> toMatrix(std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const;
+            
+            /*!
+             * Converts the ADD to a row-grouped (sparse) double matrix and the given vector to an equally row-grouped
+             * explicit vector. The given offset-labeled DDs are used to determine the correct row and column,
+             * respectively, for each entry. Note: this function assumes that the meta variables used to distinguish
+             * different row groups are at the very top of the ADD.
+             *
+             * @param vector The vector that is to be transformed to an equally grouped explicit vector.
+             * @param rowGroupSizes A vector specifying the sizes of the row groups.
+             * @param rowMetaVariables The meta variables that encode the rows of the matrix.
+             * @param columnMetaVariables The meta variables that encode the columns of the matrix.
+             * @param groupMetaVariables The meta variables that are used to distinguish different row groups.
+             * @param rowOdd The ODD used for determining the correct row.
+             * @param columnOdd The ODD used for determining the correct column.
+             * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector
+             * (if it was given).
+             */
+            std::pair<storm::storage::SparseMatrix<double>,std::vector<double>> toMatrixVector(storm::dd::Add<storm::dd::DdType::CUDD> const& vector, std::vector<uint_fast64_t>&& rowGroupSizes, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& groupMetaVariables, storm::dd::Odd<DdType::CUDD> const& rowOdd, storm::dd::Odd<DdType::CUDD> const& columnOdd) const;
+            
             /*!
              * Helper function to convert the DD into a (sparse) matrix.
              *
@@ -659,6 +720,21 @@ namespace storm {
              */
             void toMatrixRec(DdNode const* dd, std::vector<uint_fast64_t>& rowIndications, std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>>& columnsAndValues, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, Odd<DdType::CUDD> const& columnOdd, uint_fast64_t currentRowLevel, uint_fast64_t currentColumnLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, uint_fast64_t currentColumnOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices, std::vector<uint_fast64_t> const& ddColumnVariableIndices, bool generateValues = true) const;
             
+            /*!
+             * Helper function to convert the DD into a (sparse) vector.
+             *
+             * @param dd The DD to convert.
+             * @param result The vector that will hold the values upon successful completion.
+             * @param rowGroupOffsets The row offsets at which a given row group starts.
+             * @param rowOdd The ODD used for the row translation.
+             * @param currentRowLevel The currently considered row level in the DD.
+             * @param maxLevel The number of levels that need to be considered.
+             * @param currentRowOffset The current row offset.
+             * @param ddRowVariableIndices The (sorted) indices of all DD row variables that need to be considered.
+             */
+            template<typename ValueType>
+            void toVectorRec(DdNode const* dd, std::vector<ValueType>& result, std::vector<uint_fast64_t> const& rowGroupOffsets, Odd<DdType::CUDD> const& rowOdd, uint_fast64_t currentRowLevel, uint_fast64_t maxLevel, uint_fast64_t currentRowOffset, std::vector<uint_fast64_t> const& ddRowVariableIndices) const;
+            
             /*!
              * Splits the given matrix DD into the groups using the given group variables.
              *
@@ -672,7 +748,32 @@ namespace storm {
             void splitGroupsRec(DdNode* dd, std::vector<Add<DdType::CUDD>>& groups, std::vector<uint_fast64_t> const& ddGroupVariableIndices, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::set<storm::expressions::Variable> const& remainingMetaVariables) const;
             
             /*!
-             * Performs a recursive step to add the given DD-based vector to the given explicit vector.
+             * Splits the given matrix and vector DDs into the groups using the given group variables.
+             *
+             * @param dd1 The matrix DD to split.
+             * @param dd2 The vector DD to split.
+             * @param groups A vector that is to be filled with the pairs of matrix/vector DDs for the individual groups.
+             * @param ddGroupVariableIndices The (sorted) indices of all DD group variables that need to be considered.
+             * @param currentLevel The currently considered level in the DD.
+             * @param maxLevel The number of levels that need to be considered.
+             * @param remainingMetaVariables1 The meta variables that remain in the matrix DD after the groups have been split.
+             * @param remainingMetaVariables2 The meta variables that remain in the vector DD after the groups have been split.
+             */
+            void splitGroupsRec(DdNode* dd1, DdNode* dd2, std::vector<std::pair<Add<DdType::CUDD>, Add<DdType::CUDD>>>& groups, std::vector<uint_fast64_t> const& ddGroupVariableIndices, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::set<storm::expressions::Variable> const& remainingMetaVariables1, std::set<storm::expressions::Variable> const& remainingMetaVariables2) const;
+            
+            /*!
+             * Adds the current (DD-based) vector to the given explicit one.
+             *
+             * @param odd The ODD used for the translation.
+             * @param ddVariableIndices The (sorted) indices of all DD variables that need to be considered.
+             * @param targetVector The vector to which the translated DD-based vector is to be added.
+             */
+            template<typename ValueType>
+            void addToVector(Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const;
+            
+            /*!
+             * Performs a recursive step to perform the given function between the given DD-based vector to the given
+             * explicit vector.
              *
              * @param dd The DD to add to the explicit vector.
              * @param currentLevel The currently considered level in the DD.
@@ -683,7 +784,7 @@ namespace storm {
              * @param targetVector The vector to which the translated DD-based vector is to be added.
              */
             template<typename ValueType>
-            void addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const;
+            void modifyVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector, std::function<ValueType (ValueType const&, double const&)> const& function) const;
             
             /*!
              * Builds an ADD representing the given vector.
diff --git a/src/storage/dd/CuddOdd.cpp b/src/storage/dd/CuddOdd.cpp
index 46ed679b6..29301957e 100644
--- a/src/storage/dd/CuddOdd.cpp
+++ b/src/storage/dd/CuddOdd.cpp
@@ -106,6 +106,22 @@ namespace storm {
             return result;
         }
         
+        void Odd<DdType::CUDD>::expandExplicitVector(storm::dd::Odd<DdType::CUDD> const& newOdd, std::vector<double> const& oldValues, std::vector<double>& newValues) const {
+            expandValuesToVectorRec(0, *this, oldValues, 0, newOdd, newValues);
+        }
+        
+        void Odd<DdType::CUDD>::expandValuesToVectorRec(uint_fast64_t oldOffset, storm::dd::Odd<DdType::CUDD> const& oldOdd, std::vector<double> const& oldValues, uint_fast64_t newOffset, storm::dd::Odd<DdType::CUDD> const& newOdd, std::vector<double>& newValues) {
+            if (oldOdd.isTerminalNode()) {
+                STORM_LOG_THROW(newOdd.isTerminalNode(), storm::exceptions::InvalidArgumentException, "The ODDs for the translation must have the same height.");
+                if (oldOdd.getThenOffset() != 0) {
+                    newValues[newOffset] += oldValues[oldOffset];
+                }
+            } else {
+                expandValuesToVectorRec(oldOffset, oldOdd.getElseSuccessor(), oldValues, newOffset, newOdd.getElseSuccessor(), newValues);
+                expandValuesToVectorRec(oldOffset + oldOdd.getElseOffset(), oldOdd.getThenSuccessor(), oldValues, newOffset + newOdd.getElseOffset(), newOdd.getThenSuccessor(), newValues);
+            }
+        }
+        
         void Odd<DdType::CUDD>::addSelectedValuesToVectorRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, bool complement, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, uint_fast64_t currentOffset, storm::dd::Odd<DdType::CUDD> const& odd, std::vector<double>& result, uint_fast64_t& currentIndex, std::vector<double> const& values) {
             // If there are no more values to select, we can directly return.
             if (dd == Cudd_ReadLogicZero(manager.getManager()) && !complement) {
diff --git a/src/storage/dd/CuddOdd.h b/src/storage/dd/CuddOdd.h
index 4a03a3c26..e00bdd2c7 100644
--- a/src/storage/dd/CuddOdd.h
+++ b/src/storage/dd/CuddOdd.h
@@ -112,6 +112,16 @@ namespace storm {
              */
             std::vector<double> filterExplicitVector(storm::dd::Bdd<DdType::CUDD> const& selectedValues, std::vector<double> const& values) const;
             
+            /*!
+             * Adds the old values to the new values. It does so by writing the old values at their correct positions
+             * wrt. to the new ODD.
+             *
+             * @param newOdd The new ODD to use.
+             * @param oldValues The old vector of values (which is being read).
+             * @param newValues The new vector of values (which is being written).
+             */
+            void expandExplicitVector(storm::dd::Odd<DdType::CUDD> const& newOdd, std::vector<double> const& oldValues, std::vector<double>& newValues) const;
+            
         private:
             // Declare a hash functor that is used for the unique tables in the construction process.
             class HashFunctor {
@@ -159,7 +169,6 @@ namespace storm {
              */
             static std::shared_ptr<Odd<DdType::CUDD>> buildOddFromBddRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, bool complement, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<std::unordered_map<std::pair<DdNode*, bool>, std::shared_ptr<Odd<DdType::CUDD>>, HashFunctor>>& uniqueTableForLevels);
             
-            
             /*!
              * Adds the selected values the target vector.
              *
@@ -176,6 +185,19 @@ namespace storm {
              */
             static void addSelectedValuesToVectorRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, bool complement, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, uint_fast64_t currentOffset, storm::dd::Odd<DdType::CUDD> const& odd, std::vector<double>& result, uint_fast64_t& currentIndex, std::vector<double> const& values);
             
+            /*!
+             * Adds the values of the old explicit values to the new explicit values where the positions in the old vector
+             * are given by the current old ODD and the positions in the new vector are given by the new ODD.
+             *
+             * @param oldOffset The offset in the old explicit values.
+             * @param oldOdd The ODD to use for the old explicit values.
+             * @param oldValues The vector of old values.
+             * @param newOffset The offset in the new explicit values.
+             * @param newOdd The ODD to use for the new explicit values.
+             * @param newValues The vector of new values.
+             */
+            static void expandValuesToVectorRec(uint_fast64_t oldOffset, storm::dd::Odd<DdType::CUDD> const& oldOdd, std::vector<double> const& oldValues, uint_fast64_t newOffset, storm::dd::Odd<DdType::CUDD> const& newOdd, std::vector<double>& newValues);
+            
             // The then- and else-nodes.
             std::shared_ptr<Odd<DdType::CUDD>> elseNode;
             std::shared_ptr<Odd<DdType::CUDD>> thenNode;
diff --git a/src/storage/prism/Program.cpp b/src/storage/prism/Program.cpp
index 2dbc33eb9..258d4e14d 100644
--- a/src/storage/prism/Program.cpp
+++ b/src/storage/prism/Program.cpp
@@ -595,6 +595,7 @@ namespace storm {
             // Check the commands of the modules.
             bool hasProbabilisticCommand = false;
             bool hasMarkovianCommand = false;
+            bool hasLabeledMarkovianCommand = false;
             for (auto const& module : this->getModules()) {
                 std::set<storm::expressions::Variable> legalVariables = globalVariables;
                 for (auto const& variable : module.getBooleanVariables()) {
@@ -621,11 +622,7 @@ namespace storm {
                     // If the command is Markovian and labeled, we throw an error or raise a warning, depending on
                     // whether or not the PRISM compatibility mode was enabled.
                     if (command.isMarkovian() && command.isLabeled()) {
-                        if (storm::settings::generalSettings().isPrismCompatibilityEnabled()) {
-                            STORM_LOG_WARN_COND(false, "The model uses synchronizing Markovian commands. This may lead to unexpected verification results, because of unclear semantics.");
-                        } else {
-                            STORM_LOG_THROW(false, storm::exceptions::WrongFormatException, "The model uses synchronizing Markovian commands. This may lead to unexpected verification results, because of unclear semantics.");
-                        }
+                        hasLabeledMarkovianCommand = true;
                     }
                     
                     // Check all updates.
@@ -665,6 +662,14 @@ namespace storm {
                 }
             }
             
+            if (hasLabeledMarkovianCommand) {
+                if (storm::settings::generalSettings().isPrismCompatibilityEnabled()) {
+                    STORM_LOG_WARN_COND(false, "The model uses synchronizing Markovian commands. This may lead to unexpected verification results, because of unclear semantics.");
+                } else {
+                    STORM_LOG_THROW(false, storm::exceptions::WrongFormatException, "The model uses synchronizing Markovian commands. This may lead to unexpected verification results, because of unclear semantics.");
+                }
+            }
+            
             if (this->getModelType() == Program::ModelType::DTMC || this->getModelType() == Program::ModelType::MDP) {
                 STORM_LOG_THROW(!hasMarkovianCommand, storm::exceptions::WrongFormatException, "Discrete-time model must not have Markovian commands.");
             } else if (this->getModelType() == Program::ModelType::CTMC) {
diff --git a/src/utility/cli.h b/src/utility/cli.h
index 9d21378c4..351584e78 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -68,6 +68,8 @@ log4cplus::Logger printer;
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
+#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
+#include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
 
@@ -317,8 +319,12 @@ namespace storm {
                 std::string constants = settings.getConstantDefinitionString();
                 
                 bool buildRewards = false;
-                if (formula) {
+                boost::optional<std::string> rewardModelName;
+                if (formula || settings.isSymbolicRewardModelNameSet()) {
                     buildRewards = formula.get()->isRewardOperatorFormula() || formula.get()->isRewardPathFormula();
+                    if (settings.isSymbolicRewardModelNameSet()) {
+                        rewardModelName = settings.getSymbolicRewardModelName();
+                    }
                 }
                 
                 // Customize and perform model-building.
@@ -328,6 +334,8 @@ namespace storm {
                         options = typename storm::builder::ExplicitPrismModelBuilder<ValueType>::Options(*formula.get());
                     }
                     options.addConstantDefinitionsFromString(program, settings.getConstantDefinitionString());
+                    options.buildRewards = buildRewards;
+                    options.rewardModelName = rewardModelName;
                     
                     // Generate command labels if we are going to build a counterexample later.
                     if (storm::settings::counterexampleGeneratorSettings().isMinimalCommandSetGenerationSet()) {
@@ -341,7 +349,9 @@ namespace storm {
                         options = typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options(*formula.get());
                     }
                     options.addConstantDefinitionsFromString(program, settings.getConstantDefinitionString());
-                    
+                    options.buildRewards = buildRewards;
+                    options.rewardModelName = rewardModelName;
+
                     result = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
                 }
                 
@@ -503,6 +513,18 @@ namespace storm {
                     if (modelchecker.canHandle(*formula.get())) {
                         result = modelchecker.check(*formula.get());
                     }
+                } else if (model->getType() == storm::models::ModelType::Ctmc) {
+                    std::shared_ptr<storm::models::symbolic::Ctmc<DdType>> ctmc = model->template as<storm::models::symbolic::Ctmc<DdType>>();
+                    storm::modelchecker::HybridCtmcCslModelChecker<DdType, double> modelchecker(*ctmc);
+                    if (modelchecker.canHandle(*formula.get())) {
+                        result = modelchecker.check(*formula.get());
+                    }
+                } else if (model->getType() == storm::models::ModelType::Mdp) {
+                    std::shared_ptr<storm::models::symbolic::Mdp<DdType>> mdp = model->template as<storm::models::symbolic::Mdp<DdType>>();
+                    storm::modelchecker::HybridMdpPrctlModelChecker<DdType, double> modelchecker(*mdp);
+                    if (modelchecker.canHandle(*formula.get())) {
+                        result = modelchecker.check(*formula.get());
+                    }
                 } else {
                     STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This functionality is not yet implemented.");
                 }
@@ -607,7 +629,9 @@ namespace storm {
 							while (inputFileStream.good()) {
 								std::string prop;
 								std::getline(inputFileStream, prop);
-								properties.push_back(prop);
+								if (!prop.empty()) {
+									properties.push_back(prop);
+								}
 							}
 						}
 						catch (std::exception& e) {
@@ -621,14 +645,19 @@ namespace storm {
 
 					for (std::string prop : properties) {
 						boost::optional<std::shared_ptr<storm::logic::Formula>> formula;
-						if (program) {
-							storm::parser::FormulaParser formulaParser(program.get().getManager().getSharedPointer());
-							formula = formulaParser.parseFromString(prop);
-						} else {
-							storm::parser::FormulaParser formulaParser;
-							formula = formulaParser.parseFromString(prop);
+						try {
+							if (program) {
+								storm::parser::FormulaParser formulaParser(program.get().getManager().getSharedPointer());
+								formula = formulaParser.parseFromString(prop);
+							} else {
+								storm::parser::FormulaParser formulaParser;
+								formula = formulaParser.parseFromString(prop);
+							}
+							formulas.push_back(formula);
+						}
+						catch (storm::exceptions::WrongFormatException &e) {
+							STORM_LOG_WARN("Unable to parse line as formula: " << prop);
 						}
-						formulas.push_back(formula);
 					}
 					std::cout << "Parsed " << formulas.size() << " properties from file " << settings.getPropertiesFilename() << std::endl;
 				}
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 3c04feb86..fb8bdb158 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -947,7 +947,7 @@ namespace storm {
                 result.second = performProb1E(model, transitionMatrix, phiStates, psiStates, !result.first && model.getReachableStates());
                 return result;
             }
-
+            
             template <storm::dd::DdType Type>
             std::pair<storm::dd::Bdd<Type>, storm::dd::Bdd<Type>> performProb01Min(storm::models::symbolic::NondeterministicModel<Type> const& model, storm::dd::Bdd<Type> const& phiStates, storm::dd::Bdd<Type> const& psiStates) {
                 std::pair<storm::dd::Bdd<Type>, storm::dd::Bdd<Type>> result;
diff --git a/test/functional/builder/leader4.nm b/test/functional/builder/leader4.nm
new file mode 100644
index 000000000..9fbf259c4
--- /dev/null
+++ b/test/functional/builder/leader4.nm
@@ -0,0 +1,88 @@
+// asynchronous leader election
+// 4 processes
+// gxn/dxp 29/01/01
+
+mdp
+
+const int N = 4; // number of processes
+
+module process1
+	
+	// COUNTER
+	c1 : [0..N-1];
+	
+	// STATES
+	s1 : [0..4];
+	// 0  make choice
+	// 1 have not received neighbours choice
+	// 2 active
+	// 3 inactive
+	// 4 leader
+	
+	// PREFERENCE
+	p1 : [0..1];
+	
+	// VARIABLES FOR SENDING AND RECEIVING
+	receive1 : [0..2];
+	// not received anything
+	// received choice
+	// received counter
+	sent1 : [0..2];
+	// not send anything
+	// sent choice
+	// sent counter
+	
+	// pick value
+	[] (s1=0) -> 0.5 : (s1'=1) & (p1'=0) + 0.5 : (s1'=1) & (p1'=1);
+	
+	// send preference
+	[p12] (s1=1) & (sent1=0) -> (sent1'=1);
+	// receive preference
+	// stay active
+	[p41] (s1=1) & (receive1=0) & !( (p1=0) & (p4=1) ) -> (s1'=2) & (receive1'=1);
+	// become inactive
+	[p41] (s1=1) & (receive1=0) & (p1=0) & (p4=1) -> (s1'=3) & (receive1'=1);
+	
+	// send preference (can now reset preference)
+	[p12] (s1=2) & (sent1=0) -> (sent1'=1) & (p1'=0);
+	// send counter (already sent preference)
+	// not received counter yet
+	[c12] (s1=2) & (sent1=1) & (receive1=1) -> (sent1'=2);
+	// received counter (pick again)
+	[c12] (s1=2) & (sent1=1) & (receive1=2) -> (s1'=0) & (p1'=0) & (c1'=0) & (sent1'=0) & (receive1'=0);
+	
+	// receive counter and not sent yet (note in this case do not pass it on as will send own counter)
+	[c41] (s1=2) & (receive1=1) & (sent1<2) -> (receive1'=2);
+	// receive counter and sent counter
+	// only active process (decide)
+	[c41] (s1=2) & (receive1=1) & (sent1=2) & (c4=N-1) -> (s1'=4) & (p1'=0) & (c1'=0) & (sent1'=0) & (receive1'=0);
+	// other active process (pick again)
+	[c41] (s1=2) & (receive1=1) & (sent1=2) & (c4<N-1) -> (s1'=0) & (p1'=0) & (c1'=0) & (sent1'=0) & (receive1'=0);
+	
+	// send preference (must have received preference) and can now reset
+	[p12] (s1=3) & (receive1>0) & (sent1=0) -> (sent1'=1) & (p1'=0);
+	// send counter (must have received counter first) and can now reset
+	[c12] (s1=3) & (receive1=2) & (sent1=1) ->  (s1'=3) & (p1'=0) & (c1'=0) & (sent1'=0) & (receive1'=0);
+	
+	// receive preference
+	[p41] (s1=3) & (receive1=0) -> (p1'=p4) & (receive1'=1);
+	// receive counter
+	[c41] (s1=3) & (receive1=1) & (c4<N-1) -> (c1'=c4+1) & (receive1'=2);
+		
+	// done
+	[done] (s1=4) -> (s1'=s1);
+	// add loop for processes who are inactive
+	[done] (s1=3) -> (s1'=s1);
+
+endmodule
+
+module process2=process1[s1=s2,p1=p2,c1=c2,sent1=sent2,receive1=receive2,p12=p23,p41=p12,c12=c23,c41=c12,p4=p1,c4=c1] endmodule
+module process3=process1[s1=s3,p1=p3,c1=c3,sent1=sent3,receive1=receive3,p12=p34,p41=p23,c12=c34,c41=c23,p4=p2,c4=c2] endmodule
+module process4=process1[s1=s4,p1=p4,c1=c4,sent1=sent4,receive1=receive4,p12=p41,p41=p34,c12=c41,c41=c34,p4=p3,c4=c3] endmodule
+
+// reward - expected number of rounds (equals the number of times a process receives a counter)
+rewards "rounds"
+	[c12] true : 1;
+endrewards
+
+label "elected" = s1=4|s2=4|s3=4|s4=4;
diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
index 5c637afcb..4c79aa929 100644
--- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
@@ -52,26 +52,40 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) {
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult2 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2[initialState], storm::settings::generalSettings().getPrecision());
     
-    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
+    formula = formulaParser.parseFromString("P=? [ F[100,2000] !\"minimum\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult3 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(1, quantitativeCheckResult3[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3[initialState], storm::settings::generalSettings().getPrecision());
     
-    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult4 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(0, quantitativeCheckResult4[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult4[initialState], storm::settings::generalSettings().getPrecision());
     
-    formula = formulaParser.parseFromString("R=? [C<=100]");
+    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U[1,inf] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=100]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
@@ -134,7 +148,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(2.7720429852369946, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(GmmxxCtmcCslModelCheckerTest, Polling) {
diff --git a/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp
new file mode 100644
index 000000000..522e61691
--- /dev/null
+++ b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp
@@ -0,0 +1,279 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+#include "src/settings/SettingMemento.h"
+#include "src/parser/PrismParser.h"
+#include "src/parser/FormulaParser.h"
+#include "src/logic/Formulas.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/storage/dd/DdType.h"
+
+#include "src/utility/solver.h"
+#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+
+#include "src/settings/SettingsManager.h"
+
+TEST(GmmxxHybridCtmcCslModelCheckerTest, Cluster) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/cluster2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "num_repairs";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=100 !\"minimum\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(5.5461254704419085E-5, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(5.5461254704419085E-5, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F[100,100] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F[100,2000] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U[1,inf] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=100]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/embedded2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "up";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=10000 \"down\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.0019216435246119591, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.0019216435246119591, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_actuators\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(3.7079151806696567E-6, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.7079151806696567E-6, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_io\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.001556839327673734, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.001556839327673734, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_sensors\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(4.429620626755424E-5, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.429620626755424E-5, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=10000]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/polling2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=?[ F<=10 \"target\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(GmmxxHybridCtmcCslModelCheckerTest, Fms) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // No properties to check at this point.
+}
+
+TEST(GmmxxHybridCtmcCslModelCheckerTest, Tandem) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/tandem5.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model with the customers reward structure.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "customers";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=10 \"network_full\" ]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.015446370562428037, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.015446370562428037, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F<=10 \"first_queue_full\" ]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.999999837225515, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.999999837225515, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [\"second_queue_full\" U<=1 !\"second_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [I=10]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(5.679243850315877, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(5.679243850315877, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=10]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(55.44792186036232, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(55.44792186036232, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [F \"first_queue_full\"&\"second_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(262.85103498583413, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
+}
diff --git a/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp
new file mode 100644
index 000000000..3a7cd5402
--- /dev/null
+++ b/test/functional/modelchecker/GmmxxHybridMdpPrctlModelCheckerTest.cpp
@@ -0,0 +1,190 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/settings/SettingsManager.h"
+
+TEST(GmmxxHybridMdpPrctlModelCheckerTest, Dice) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/two_dice.nm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "coinflips";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(169, model->getNumberOfStates());
+    EXPECT_EQ(436, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxMinMaxLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("four");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+    
+    result = checker.check(*minRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult7 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+    
+    result = checker.check(*maxRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult8 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(GmmxxHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader4.nm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "rounds";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(3172, model->getNumberOfStates());
+    EXPECT_EQ(7144, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxMinMaxLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 25);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, boundedUntilFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, boundedUntilFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+    
+    result = checker.check(*minRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+    
+    result = checker.check(*maxRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
diff --git a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
similarity index 92%
rename from test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
rename to test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
index 9db619115..688acda91 100644
--- a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
@@ -51,27 +51,41 @@ TEST(SparseCtmcCslModelCheckerTest, Cluster) {
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult2 = checkResult->asExplicitQuantitativeCheckResult<double>();
     EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2[initialState], storm::settings::generalSettings().getPrecision());
-
-    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
+    
+    formula = formulaParser.parseFromString("P=? [ F[100,2000] !\"minimum\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult3 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(1, quantitativeCheckResult3[initialState], storm::settings::generalSettings().getPrecision());
-
-    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult4 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(0, quantitativeCheckResult4[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult4[initialState], storm::settings::generalSettings().getPrecision());
     
-    formula = formulaParser.parseFromString("R=? [C<=100]");
+    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
     checkResult = modelchecker.check(*formula);
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U[1,inf] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult6 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6[initialState], storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=100]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
+    storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult7 = checkResult->asExplicitQuantitativeCheckResult<double>();
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(SparseCtmcCslModelCheckerTest, Embedded) {
@@ -133,7 +147,7 @@ TEST(SparseCtmcCslModelCheckerTest, Embedded) {
     
     ASSERT_TRUE(checkResult->isExplicitQuantitativeCheckResult());
     storm::modelchecker::ExplicitQuantitativeCheckResult<double> quantitativeCheckResult5 = checkResult->asExplicitQuantitativeCheckResult<double>();
-    EXPECT_NEAR(2.7720429852369946, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5[initialState], storm::settings::generalSettings().getPrecision());
 }
 
 TEST(SparseCtmcCslModelCheckerTest, Polling) {
diff --git a/test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
similarity index 100%
rename from test/functional/modelchecker/SparseDtmcPrctlModelCheckerTest.cpp
rename to test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
diff --git a/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp
new file mode 100644
index 000000000..ab90e8960
--- /dev/null
+++ b/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp
@@ -0,0 +1,279 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+#include "src/settings/SettingMemento.h"
+#include "src/parser/PrismParser.h"
+#include "src/parser/FormulaParser.h"
+#include "src/logic/Formulas.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/storage/dd/DdType.h"
+
+#include "src/utility/solver.h"
+#include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+
+#include "src/settings/SettingsManager.h"
+
+TEST(NativeHybridCtmcCslModelCheckerTest, Cluster) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/cluster2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "num_repairs";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=100 !\"minimum\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(5.5461254704419085E-5, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(5.5461254704419085E-5, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F[100,100] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.3397873548343415E-6, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F[100,2000] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.001105335651670241, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U<=10 \"premium\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"minimum\" U[1,inf] \"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ \"minimum\" U[1,inf] !\"minimum\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9999999033633374, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=100]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>  quantitativeCheckResult7 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.8602815057967503, quantitativeCheckResult7.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(NativeHybridCtmcCslModelCheckerTest, Embedded) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/embedded2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "up";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=10000 \"down\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.0019216435246119591, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.0019216435246119591, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_actuators\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(3.7079151806696567E-6, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.7079151806696567E-6, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_io\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.001556839327673734, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.001556839327673734, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ !\"down\" U<=10000 \"fail_sensors\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(4.429620626755424E-5, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.429620626755424E-5, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=10000]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.7745274082080154, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(NativeHybridCtmcCslModelCheckerTest, Polling) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/polling2.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model.
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=?[ F<=10 \"target\"]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+}
+
+TEST(NativeHybridCtmcCslModelCheckerTest, Fms) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // No properties to check at this point.
+}
+
+TEST(NativeHybridCtmcCslModelCheckerTest, Tandem) {
+    // Set the PRISM compatibility mode temporarily. It is set to its old value once the returned object is destructed.
+    std::unique_ptr<storm::settings::SettingMemento> enablePrismCompatibility = storm::settings::mutableGeneralSettings().overridePrismCompatibilityMode(true);
+    
+    // Parse the model description.
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/tandem5.sm");
+    storm::parser::FormulaParser formulaParser(program.getManager().getSharedPointer());
+    std::shared_ptr<storm::logic::Formula> formula(nullptr);
+    
+    // Build the model with the customers reward structure.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "customers";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    ASSERT_EQ(storm::models::ModelType::Ctmc, model->getType());
+    std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
+    
+    // Create model checker.
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    // Start checking properties.
+    formula = formulaParser.parseFromString("P=? [ F<=10 \"network_full\" ]");
+    std::unique_ptr<storm::modelchecker::CheckResult> checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult1 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.015446370562428037, quantitativeCheckResult1.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.015446370562428037, quantitativeCheckResult1.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [ F<=10 \"first_queue_full\" ]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult2 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(0.999999837225515, quantitativeCheckResult2.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.999999837225515, quantitativeCheckResult2.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("P=? [\"second_queue_full\" U<=1 !\"second_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult3 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(1, quantitativeCheckResult3.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeCheckResult3.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [I=10]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult4 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(5.679243850315877, quantitativeCheckResult4.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(5.679243850315877, quantitativeCheckResult4.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [C<=10]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult5 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(55.44792186036232, quantitativeCheckResult5.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(55.44792186036232, quantitativeCheckResult5.getMax(), storm::settings::generalSettings().getPrecision());
+    
+    formula = formulaParser.parseFromString("R=? [F \"first_queue_full\"&\"second_queue_full\"]");
+    checkResult = modelchecker.check(*formula);
+    
+    ASSERT_TRUE(checkResult->isHybridQuantitativeCheckResult());
+    checkResult->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(ctmc->getReachableStates(), ctmc->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD> quantitativeCheckResult6 = checkResult->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    EXPECT_NEAR(262.78571505691389, quantitativeCheckResult6.getMin(), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(262.78571505691389, quantitativeCheckResult6.getMax(), storm::settings::generalSettings().getPrecision());
+}
diff --git a/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp
new file mode 100644
index 000000000..8cd46bcb1
--- /dev/null
+++ b/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp
@@ -0,0 +1,165 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/settings/SettingsManager.h"
+
+TEST(NativeHybridDtmcPrctlModelCheckerTest, Die) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/die.pm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "coin_flips";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(13, model->getNumberOfStates());
+    EXPECT_EQ(20, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("one");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult1.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult2.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0/6.0, quantitativeResult3.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    auto done = std::make_shared<storm::logic::AtomicLabelFormula>("done");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(done);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(3.6666646003723145, quantitativeResult4.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(3.6666646003723145, quantitativeResult4.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+}
+
+TEST(NativeHybridDtmcPrctlModelCheckerTest, Crowds) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/crowds-5-5.pm");
+    
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program);
+    EXPECT_EQ(8607, model->getNumberOfStates());
+    EXPECT_EQ(15113, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observe0Greater1");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.33288205191646525, quantitativeResult1.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.33288205191646525, quantitativeResult1.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeIGreater1");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.15222066094730619, quantitativeResult2.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.15222066094730619, quantitativeResult2.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observeOnlyTrueSender");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.32153900158185761, quantitativeResult3.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.32153900158185761, quantitativeResult3.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+}
+
+TEST(NativeHybridDtmcPrctlModelCheckerTest, SynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader-3-5.pm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "num_rounds";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(273, model->getNumberOfStates());
+    EXPECT_EQ(397, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Dtmc);
+    
+    std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*eventuallyFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0, quantitativeResult1.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0, quantitativeResult1.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 20);
+    
+    result = checker.check(*boundedUntilFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.99999989760000074, quantitativeResult2.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    
+    result = checker.check(*reachabilityRewardFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMin(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1.0416666666666643, quantitativeResult3.getMax(), storm::settings::gmmxxEquationSolverSettings().getPrecision());
+}
+
diff --git a/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp
new file mode 100644
index 000000000..1b1af244d
--- /dev/null
+++ b/test/functional/modelchecker/NativeHybridMdpPrctlModelCheckerTest.cpp
@@ -0,0 +1,190 @@
+#include "gtest/gtest.h"
+#include "storm-config.h"
+
+#include "src/logic/Formulas.h"
+#include "src/utility/solver.h"
+#include "src/modelchecker/prctl/HybridMdpPrctlModelChecker.h"
+#include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
+#include "src/modelchecker/results/SymbolicQuantitativeCheckResult.h"
+#include "src/parser/PrismParser.h"
+#include "src/builder/DdPrismModelBuilder.h"
+#include "src/models/symbolic/Dtmc.h"
+#include "src/settings/SettingsManager.h"
+
+TEST(NativeHybridMdpPrctlModelCheckerTest, Dice) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/two_dice.nm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "coinflips";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(169, model->getNumberOfStates());
+    EXPECT_EQ(436, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0277777612209320068, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("three");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0555555224418640136, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("four");
+    eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.083333283662796020508, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+    
+    result = checker.check(*minRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult7 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult7.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+    
+    result = checker.check(*maxRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult8 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(7.3333283960819244, quantitativeResult8.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
+
+TEST(NativeHybridMdpPrctlModelCheckerTest, AsynchronousLeader) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/leader4.nm");
+    
+    // Build the die model with its reward model.
+    typename storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::Options options;
+    options.buildRewards = true;
+    options.rewardModelName = "rounds";
+    std::shared_ptr<storm::models::symbolic::Model<storm::dd::DdType::CUDD>> model = storm::builder::DdPrismModelBuilder<storm::dd::DdType::CUDD>::translateProgram(program, options);
+    EXPECT_EQ(3172, model->getNumberOfStates());
+    EXPECT_EQ(7144, model->getNumberOfTransitions());
+    
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    
+    std::shared_ptr<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>> mdp = model->as<storm::models::symbolic::Mdp<storm::dd::DdType::CUDD>>();
+    
+    storm::modelchecker::HybridMdpPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*mdp, std::unique_ptr<storm::utility::solver::MinMaxLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeMinMaxLinearEquationSolverFactory<double>()));
+    
+    auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
+    auto minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, eventuallyFormula);
+    
+    std::unique_ptr<storm::modelchecker::CheckResult> result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult1 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult1.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult1.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, eventuallyFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::SymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult2 = result->asSymbolicQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(1, quantitativeResult2.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(1, quantitativeResult2.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto trueFormula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
+    auto boundedUntilFormula = std::make_shared<storm::logic::BoundedUntilFormula>(trueFormula, labelFormula, 25);
+    minProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Minimize, boundedUntilFormula);
+    
+    result = checker.check(*minProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult3 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult3.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    maxProbabilityOperatorFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(storm::logic::OptimalityType::Maximize, boundedUntilFormula);
+    
+    result = checker.check(*maxProbabilityOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult4 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(0.0625, quantitativeResult4.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
+    auto reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
+    auto minRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Minimize, reachabilityRewardFormula);
+    
+    result = checker.check(*minRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult5 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856925589077264, quantitativeResult5.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    
+    auto maxRewardOperatorFormula = std::make_shared<storm::logic::RewardOperatorFormula>(storm::logic::OptimalityType::Maximize, reachabilityRewardFormula);
+    
+    result = checker.check(*maxRewardOperatorFormula);
+    result->filter(storm::modelchecker::SymbolicQualitativeCheckResult<storm::dd::DdType::CUDD>(model->getReachableStates(), model->getInitialStates()));
+    storm::modelchecker::HybridQuantitativeCheckResult<storm::dd::DdType::CUDD>& quantitativeResult6 = result->asHybridQuantitativeCheckResult<storm::dd::DdType::CUDD>();
+    
+    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMin(), storm::settings::nativeEquationSolverSettings().getPrecision());
+    EXPECT_NEAR(4.2856953906798676, quantitativeResult6.getMax(), storm::settings::nativeEquationSolverSettings().getPrecision());
+}
diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeMdpPrctlModelCheckerTest.cpp
similarity index 100%
rename from test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
rename to test/functional/modelchecker/NativeMdpPrctlModelCheckerTest.cpp
diff --git a/test/functional/storage/CuddDdTest.cpp b/test/functional/storage/CuddDdTest.cpp
index 6255e912c..f8b149c2d 100644
--- a/test/functional/storage/CuddDdTest.cpp
+++ b/test/functional/storage/CuddDdTest.cpp
@@ -351,7 +351,7 @@ TEST(CuddDd, AddOddTest) {
     EXPECT_EQ(25, matrix.getNonzeroEntryCount());
     
     dd = manager->getRange(x.first).toAdd() * manager->getRange(x.second).toAdd() * manager->getEncoding(a.first, 0).toAdd().ite(dd, dd + manager->getConstant(1));
-    ASSERT_NO_THROW(matrix = dd.toMatrix({x.first}, {x.second}, {a.first}, rowOdd, columnOdd));
+    ASSERT_NO_THROW(matrix = dd.toMatrix({a.first}, rowOdd, columnOdd));
     EXPECT_EQ(18, matrix.getRowCount());
     EXPECT_EQ(9, matrix.getRowGroupCount());
     EXPECT_EQ(9, matrix.getColumnCount());
@@ -398,7 +398,7 @@ TEST(CuddDd, BddOddTest) {
     EXPECT_EQ(25, matrix.getNonzeroEntryCount());
     
     dd = manager->getRange(x.first).toAdd() * manager->getRange(x.second).toAdd() * manager->getEncoding(a.first, 0).toAdd().ite(dd, dd + manager->getConstant(1));
-    ASSERT_NO_THROW(matrix = dd.toMatrix({x.first}, {x.second}, {a.first}, rowOdd, columnOdd));
+    ASSERT_NO_THROW(matrix = dd.toMatrix({a.first}, rowOdd, columnOdd));
     EXPECT_EQ(18, matrix.getRowCount());
     EXPECT_EQ(9, matrix.getRowGroupCount());
     EXPECT_EQ(9, matrix.getColumnCount());