diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
index 4e464fd84..d39e90051 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
@@ -5,7 +5,8 @@
  * Created on September 9, 2015, 12:34 PM
  */
 
-#include "AbstractSparseRegionModelChecker.h"
+#include "src/modelchecker/region/AbstractSparseRegionModelChecker.h"
+
 #include "src/adapters/CarlAdapter.h"
 #include "src/modelchecker/region/RegionCheckResult.h"
 #include "src/logic/Formulas.h"
@@ -22,6 +23,7 @@
 #include "src/exceptions/InvalidSettingsException.h"
 #include "src/exceptions/NotImplementedException.h"
 #include "src/exceptions/UnexpectedException.h"
+#include "utility/ConversionHelper.h"
 
 namespace storm {
     namespace modelchecker {
@@ -45,39 +47,40 @@ namespace storm {
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            bool const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::isResultConstant() const {
-                return this->resultConstant;
+            std::shared_ptr<storm::logic::OperatorFormula> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSpecifiedFormula() const {
+                return specifiedFormula;
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            std::shared_ptr<ParametricSparseModelType> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSimpleModel() const {
-                return this->simpleModel;
+            ConstantType AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSpecifiedFormulaBound() const {
+                return storm::utility::region::convertNumber<ConstantType>(this->getSpecifiedFormula()->getBound());
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            bool const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::isComputeRewards() const {
-                return computeRewards;
+            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::specifiedFormulaHasUpperBound() const {
+                return !storm::logic::isLowerBound(this->getSpecifiedFormula()->getComparisonType());
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            std::shared_ptr<storm::logic::Formula> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSpecifiedFormula() const {
-                return specifiedFormula;
+            bool const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::isComputeRewards() const {
+                return computeRewards;
             }
-
+            
             template<typename ParametricSparseModelType, typename ConstantType>
-            ConstantType const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSpecifiedFormulaBound() const {
-                return specifiedFormulaBound;
+            bool const AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::isResultConstant() const {
+                return this->constantResult.operator bool();
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            storm::logic::ComparisonType const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSpecifiedFormulaCompType() const {
-                return specifiedFormulaCompType;
+            std::shared_ptr<ParametricSparseModelType> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSimpleModel() const {
+                return this->simpleModel;
             }
             
             template<typename ParametricSparseModelType, typename ConstantType>
-            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::specifiedFormulaHasUpperBound() const {
-                return !storm::logic::isLowerBound(this->getSpecifiedFormulaCompType());
+            std::shared_ptr<storm::logic::OperatorFormula> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSimpleFormula() const {
+                return this->simpleFormula;
             }
+
             
 
             
@@ -85,14 +88,23 @@ namespace storm {
             void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::specifyFormula(std::shared_ptr<storm::logic::Formula> formula) {
                 std::chrono::high_resolution_clock::time_point timeSpecifyFormulaStart = std::chrono::high_resolution_clock::now();
                 STORM_LOG_THROW(this->canHandle(*formula), storm::exceptions::InvalidArgumentException, "Tried to specify a formula that can not be handled.");
-
                 //Initialize the context for this formula
-                this->specifiedFormula = formula;
-                this->resultConstant=false;
-                this->simpleFormula=nullptr;
-                this->isApproximationApplicable=false;
-                this->approximationModel=nullptr;
-                this->samplingModel=nullptr;
+                if (formula->isProbabilityOperatorFormula()) {
+                    this->specifiedFormula = std::make_shared<storm::logic::ProbabilityOperatorFormula>(formula->asProbabilityOperatorFormula());
+                    this->computeRewards = false;
+                }
+                else if (formula->isRewardOperatorFormula()) {
+                    this->specifiedFormula =  std::make_shared<storm::logic::RewardOperatorFormula>(formula->asRewardOperatorFormula());
+                    this->computeRewards=true;
+                }
+                else {
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The specified property " << this->getSpecifiedFormula() << "is not supported");
+                }
+                this->constantResult = boost::none;
+                this->simpleFormula = nullptr;
+                this->isApproximationApplicable = false;
+                this->approximationModel = nullptr;
+                this->samplingModel = nullptr;
                 //stuff for statistics:
                 this->numOfCheckedRegions=0;
                 this->numOfRegionsSolvedThroughSampling=0;
@@ -109,33 +121,25 @@ namespace storm {
                 this->timeComputeReachabilityFunction=std::chrono::high_resolution_clock::duration::zero();
                 this->timeApproxModelInstantiation=std::chrono::high_resolution_clock::duration::zero();
 
-                // set some information regarding the formula.
-                if (this->getSpecifiedFormula()->isProbabilityOperatorFormula()) {
-                    this->computeRewards = false;
-                    this->specifiedFormulaCompType = this->getSpecifiedFormula()->asProbabilityOperatorFormula().getComparisonType();
-                    this->specifiedFormulaBound = storm::utility::region::convertNumber<ConstantType>(this->getSpecifiedFormula()->asProbabilityOperatorFormula().getBound());
-                }
-                else if (this->getSpecifiedFormula()->isRewardOperatorFormula()) {
-                    this->computeRewards=true;
-                    this->specifiedFormulaCompType = this->getSpecifiedFormula()->asRewardOperatorFormula().getComparisonType();
-                    this->specifiedFormulaBound = this->getSpecifiedFormula()->asRewardOperatorFormula().getBound();
-                }
-                else {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidPropertyException, "The specified property " << this->getSpecifiedFormula() << "is not supported");
-                }
                 
                 std::chrono::high_resolution_clock::time_point timePreprocessingStart = std::chrono::high_resolution_clock::now();
-                this->preprocess(this->simpleModel, this->simpleFormula, isApproximationApplicable, resultConstant);
+                this->preprocess(this->simpleModel, this->simpleFormula, isApproximationApplicable, constantResult);
                 std::chrono::high_resolution_clock::time_point timePreprocessingEnd = std::chrono::high_resolution_clock::now();
                 //Check if the approximation and the sampling model needs to be computed
                 if(!this->isResultConstant()){
                     if(this->isApproximationApplicable && storm::settings::regionSettings().doApprox()){
-                        initializeApproximationModel(*this->getSimpleModel(), this->simpleFormula);
+                        initializeApproximationModel(*this->getSimpleModel(), this->getSimpleFormula());
                     }
                     if(storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE ||
                             (!storm::settings::regionSettings().doSample() && storm::settings::regionSettings().getApproxMode()==storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST)){
-                        initializeSamplingModel(*this->getSimpleModel(), this->simpleFormula);
+                        initializeSamplingModel(*this->getSimpleModel(), this->getSimpleFormula());
                     }
+                } else if (this->isResultConstant() && this->constantResult.get() == storm::utility::region::convertNumber<ConstantType>(-1.0)){
+                    //In this case, the result is constant but has not been computed yet. so do it now!
+                    initializeSamplingModel(*this->getSimpleModel(), this->getSimpleFormula());
+                    std::map<VariableType, CoefficientType> emptySubstitution;
+                    this->getSamplingModel()->instantiate(emptySubstitution);
+                    this->constantResult = this->getSamplingModel()->computeValues()[*this->getSamplingModel()->getModel()->getInitialStates().begin()];
                 }
 
                 //some more information for statistics...
@@ -145,7 +149,7 @@ namespace storm {
             }
             
             template<typename ParametricSparseModelType, typename ConstantType>
-            void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::Formula> formula) {
+            void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula) {
                 std::chrono::high_resolution_clock::time_point timeInitApproxModelStart = std::chrono::high_resolution_clock::now();
                 STORM_LOG_THROW(this->isApproximationApplicable, storm::exceptions::UnexpectedException, "Approximation model requested but approximation is not applicable");
                 this->approximationModel=std::make_shared<ApproximationModel<ParametricSparseModelType, ConstantType>>(model, formula);
@@ -155,7 +159,7 @@ namespace storm {
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeSamplingModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::Formula> formula) {
+            void AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::initializeSamplingModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula) {
                 std::chrono::high_resolution_clock::time_point timeInitSamplingModelStart = std::chrono::high_resolution_clock::now();
                 this->samplingModel=std::make_shared<SamplingModel<ParametricSparseModelType, ConstantType>>(model, formula);
                 std::chrono::high_resolution_clock::time_point timeInitSamplingModelEnd = std::chrono::high_resolution_clock::now();
@@ -269,7 +273,7 @@ namespace storm {
             std::shared_ptr<ApproximationModel<ParametricSparseModelType, ConstantType>> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getApproximationModel() {
                 if(this->approximationModel==nullptr){
                     STORM_LOG_WARN("Approximation model requested but it has not been initialized when specifying the formula. Will initialize it now.");
-                    initializeApproximationModel(*this->getSimpleModel(), this->simpleFormula);
+                    initializeApproximationModel(*this->getSimpleModel(), this->getSimpleFormula());
                 }
                 return this->approximationModel;
             }
@@ -285,11 +289,20 @@ namespace storm {
                 return false;
             }
             
+            template<typename ParametricSparseModelType, typename ConstantType>
+            ConstantType AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getReachabilityValue(std::map<VariableType, CoefficientType> const& point) {
+                if(this->isResultConstant()){
+                    return this->constantResult.get();
+                }
+                this->getSamplingModel()->instantiate(point);
+                return this->getSamplingModel()->computeValues()[*this->getSamplingModel()->getModel()->getInitialStates().begin()];
+            }
+            
             template<typename ParametricSparseModelType, typename ConstantType>
             std::shared_ptr<SamplingModel<ParametricSparseModelType, ConstantType>> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSamplingModel() {
                 if(this->samplingModel==nullptr){
                     STORM_LOG_WARN("Sampling model requested but it has not been initialized when specifying the formula. Will initialize it now.");
-                    initializeSamplingModel(*this->getSimpleModel(), this->simpleFormula);
+                    initializeSamplingModel(*this->getSimpleModel(), this->getSimpleFormula());
                 }
                 return this->samplingModel;
             }
@@ -297,7 +310,7 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::valueIsInBoundOfFormula(ConstantType const& value){
                 STORM_LOG_THROW(this->getSpecifiedFormula()!=nullptr, storm::exceptions::InvalidStateException, "Tried to compare a value to the bound of a formula, but no formula specified.");
-                switch (this->getSpecifiedFormulaCompType()) {
+                switch (this->getSpecifiedFormula()->getComparisonType()) {
                     case storm::logic::ComparisonType::Greater:
                         return (value > this->getSpecifiedFormulaBound());
                     case storm::logic::ComparisonType::GreaterEqual:
@@ -381,8 +394,10 @@ namespace storm {
             }
         
         
+        //note: for other template instantiations, add rules for the typedefs of VariableType and CoefficientType in utility/regions.h
 #ifdef STORM_HAVE_CARL
         template class AbstractSparseRegionModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
+        template class AbstractSparseRegionModelChecker<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
 #endif
         } // namespace region
     } //namespace modelchecker
diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.h b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
index e4bb9b08c..70e5d70bc 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.h
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
@@ -9,6 +9,7 @@
 #define	STORM_MODELCHECKER_REGION_ABSTRACTSPARSEREGIONMODELCHECKER_H
 
 #include <ostream>
+#include <boost/optional.hpp>
 
 #include "src/utility/region.h"
 #include "src/modelchecker/region/ParameterRegion.h"
@@ -18,6 +19,7 @@
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/models/sparse/Model.h"
 #include "src/models/sparse/Dtmc.h"
+#include "src/models/sparse/Mdp.h"
 #include "src/logic/Formulas.h"
 
 namespace storm {
@@ -77,9 +79,8 @@ namespace storm {
                  * Returns the reachability Value at the specified point by instantiating the sampling model. 
                  * 
                  * @param point The point (i.e. parameter evaluation) at which to compute the reachability value.
-                 * @param evaluateFunction If set, the reachability function is evaluated. Otherwise, the sampling model is instantiated.
                  */
-                virtual ConstantType getReachabilityValue(std::map<VariableType, CoefficientType>const& point, bool evaluateFunction=false) = 0;
+                ConstantType getReachabilityValue(std::map<VariableType, CoefficientType>const& point);
                 
                 /*!
                  * Returns true iff the given value satisfies the bound given by the specified property
@@ -100,25 +101,25 @@ namespace storm {
             protected:
                 
                 /*!
-                 * returns the considered model
+                 * some trivial getters
                  */
                 ParametricSparseModelType const& getModel() const;
-                bool const& isResultConstant() const;
-                std::shared_ptr<ParametricSparseModelType> const& getSimpleModel() const;
-                bool const& isComputeRewards() const;
-                std::shared_ptr<storm::logic::Formula> const& getSpecifiedFormula() const;
-                ConstantType const& getSpecifiedFormulaBound() const;
-                storm::logic::ComparisonType const& getSpecifiedFormulaCompType() const;
+                std::shared_ptr<storm::logic::OperatorFormula> const& getSpecifiedFormula() const;
+                ConstantType getSpecifiedFormulaBound() const;
                 bool specifiedFormulaHasUpperBound() const;
+                bool const& isComputeRewards() const;
+                bool const isResultConstant() const;
+                std::shared_ptr<ParametricSparseModelType> const& getSimpleModel() const;
+                std::shared_ptr<storm::logic::OperatorFormula> const& getSimpleFormula() const;
                 
                 /*!
                  * Makes the required preprocessing steps for the specified model and formula
                  * Computes a simplified version of the model and formula that can be analyzed more efficiently.
                  * Also checks whether the approximation technique is applicable and whether the result is constant.
-                 * 
+                 * In the latter case, the result is already computed and set to the given parameter. (otherwise the parameter is not touched).
                  * @note this->specifiedFormula and this->computeRewards has to be set accordingly, before calling this function
                  */
-                virtual void preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::Formula>& simpleFormula, bool& isApproximationApplicable, bool& isResultConstant) = 0;
+                virtual void preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult) = 0;
                 
                 /*!
                  * Instantiates the approximation model to compute bounds on the maximal/minimal reachability probability (or reachability reward).
@@ -151,11 +152,9 @@ namespace storm {
                  * May set the satPoint and violatedPoint of the regions if thy are not yet specified and such point is given.
                  * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
                  * 
-                 * @param favorViaFunction if not stated otherwise (e.g. in the settings), the sampling will be done via the
-                 *                          reachabilityFunction if this flag is true. If the flag is false, sampling will be 
-                 *                          done via instantiation of the samplingmodel. Note that this argument is ignored,
-                 *                          unless sampling has been turned of in the settings
-                 * 
+                 * @param favorViaFunction If sampling has been turned off in the settings and a computation via evaluating
+                 *                         the reachability function is possible, this flag decides whether to instantiate the
+                 *                         sampling model or evaluate the function.
                  * @return true if an violated point as well as a sat point has been found, i.e., the check result is changed to EXISTSOTH
                  */
                 virtual bool checkPoint(ParameterRegion<ParametricType>& region, std::map<VariableType, CoefficientType>const& point, bool favorViaFunction=false) = 0;
@@ -182,25 +181,23 @@ namespace storm {
                  * 
                  * @note does not check whether approximation can be applied
                  */
-                void initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::Formula> formula);
+                void initializeApproximationModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula);
 
                 /*!
                  * initializes the Sampling Model
                  */
-                void initializeSamplingModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::Formula> formula);
+                void initializeSamplingModel(ParametricSparseModelType const& model, std::shared_ptr<storm::logic::OperatorFormula> formula);
                 
                 // The model this model checker is supposed to analyze.
                 ParametricSparseModelType const& model;
                 //The currently specified formula
-                std::shared_ptr<storm::logic::Formula> specifiedFormula;
-                storm::logic::ComparisonType specifiedFormulaCompType;
-                ConstantType specifiedFormulaBound;
+                std::shared_ptr<storm::logic::OperatorFormula> specifiedFormula;
                 //A flag that is true iff we are interested in rewards
                 bool computeRewards;
                 // the original model after states with constant transitions have been eliminated
                 std::shared_ptr<ParametricSparseModelType> simpleModel;
                 // a formula that can be checked on the simplified model
-                std::shared_ptr<storm::logic::Formula> simpleFormula;
+                std::shared_ptr<storm::logic::OperatorFormula> simpleFormula;
                 // a flag that is true if approximation is applicable, i.e., there are only linear functions at transitions of the model
                 bool isApproximationApplicable;
                 // the model that  is used to approximate the reachability values
@@ -208,7 +205,7 @@ namespace storm {
                 // the model that can be instantiated to check the value at a certain point
                 std::shared_ptr<SamplingModel<ParametricSparseModelType, ConstantType>> samplingModel;
                 // a flag that is true iff the resulting reachability function is constant
-                bool resultConstant;
+                boost::optional<ConstantType> constantResult;
                 
                 // runtimes and other information for statistics. 
                 uint_fast64_t numOfCheckedRegions;
diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp
index 546f7f556..aba216cc8 100644
--- a/src/modelchecker/region/ApproximationModel.cpp
+++ b/src/modelchecker/region/ApproximationModel.cpp
@@ -20,10 +20,10 @@ namespace storm {
         namespace region {
         
             template<typename ParametricSparseModelType, typename ConstantType>
-            ApproximationModel<ParametricSparseModelType, ConstantType>::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::Formula> formula) : formula(formula){
-                if(this->formula->isEventuallyFormula()){
+            ApproximationModel<ParametricSparseModelType, ConstantType>::ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula) : formula(formula){
+                if(this->formula->isProbabilityOperatorFormula()){
                     this->computeRewards=false;
-                } else if(this->formula->isReachabilityRewardFormula()){
+                } else if(this->formula->isRewardOperatorFormula()){
                     this->computeRewards=true;
                     STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should be unique");
                     STORM_LOG_THROW(parametricModel.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the approximation model should have state rewards only");
@@ -328,11 +328,11 @@ namespace storm {
                     }
                     //perform model checking on the mdp
                     boost::optional<std::string> noRewardModelName; //it should be uniquely given
-                    resultPtr = modelChecker.computeReachabilityRewards(this->formula->asReachabilityRewardFormula(), noRewardModelName, false, optDir);
+                    resultPtr = modelChecker.computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula(), noRewardModelName, false, optDir);
                 }
                 else {
                     //perform model checking on the mdp
-                    resultPtr = modelChecker.computeEventuallyProbabilities(this->formula->asEventuallyFormula(), false, optDir);
+                    resultPtr = modelChecker.computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula(), false, optDir);
                 }
                 return resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
             }
@@ -340,6 +340,7 @@ namespace storm {
 
 #ifdef STORM_HAVE_CARL
             template class ApproximationModel<storm::models::sparse::Dtmc<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
+            template class ApproximationModel<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
 #endif
         } //namespace region
     }
diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h
index d861209ff..656efe1a8 100644
--- a/src/modelchecker/region/ApproximationModel.h
+++ b/src/modelchecker/region/ApproximationModel.h
@@ -31,7 +31,7 @@ namespace storm {
                 /*!
                  * @note this will not check whether approximation is applicable
                  */
-                ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::Formula> formula);
+                ApproximationModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula);
                 virtual ~ApproximationModel();
 
                 /*!
@@ -129,7 +129,7 @@ namespace storm {
                 //The Model with which we work
                 std::shared_ptr<storm::models::sparse::Mdp<ConstantType>> model;
                 //The formula for which we will compute the values
-                std::shared_ptr<storm::logic::Formula> formula;
+                std::shared_ptr<storm::logic::OperatorFormula> formula;
                 //A flag that denotes whether we compute probabilities or rewards
                 bool computeRewards;
 
diff --git a/src/modelchecker/region/SamplingModel.cpp b/src/modelchecker/region/SamplingModel.cpp
index 6880d24fc..c50b74d26 100644
--- a/src/modelchecker/region/SamplingModel.cpp
+++ b/src/modelchecker/region/SamplingModel.cpp
@@ -6,7 +6,11 @@
  */
 
 #include "src/modelchecker/region/SamplingModel.h"
+
+#include "src/models/sparse/Dtmc.h"
+#include "src/models/sparse/Mdp.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
+#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/utility/macros.h"
 #include "src/utility/region.h"
@@ -14,16 +18,18 @@
 #include "src/exceptions/UnexpectedException.h"
 #include "src/exceptions/InvalidArgumentException.h"
 #include "models/sparse/StandardRewardModel.h"
+#include "cuddObj.hh"
+#include "exceptions/IllegalArgumentException.h"
 
 namespace storm {
     namespace modelchecker {
         namespace region {
         
             template<typename ParametricSparseModelType, typename ConstantType>
-            SamplingModel<ParametricSparseModelType, ConstantType>::SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::Formula> formula) : formula(formula){
-                if(this->formula->isEventuallyFormula()){
+            SamplingModel<ParametricSparseModelType, ConstantType>::SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula) : formula(formula){
+                if(this->formula->isProbabilityOperatorFormula()){
                     this->computeRewards=false;
-                } else if(this->formula->isReachabilityRewardFormula()){
+                } else if(this->formula->isRewardOperatorFormula()){
                     this->computeRewards=true;
                     STORM_LOG_THROW(parametricModel.hasUniqueRewardModel(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should be unique");
                     STORM_LOG_THROW(parametricModel.getUniqueRewardModel()->second.hasOnlyStateRewards(), storm::exceptions::InvalidArgumentException, "The rewardmodel of the sampling model should have state rewards only");
@@ -48,9 +54,17 @@ namespace storm {
 
                 //Obtain other model ingredients and the sampling model itself
                 storm::models::sparse::StateLabeling labeling(parametricModel.getStateLabeling());
-                boost::optional<storm::storage::SparseMatrix<ConstantType>> noTransitionRewards;
                 boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> noChoiceLabeling;  
-                this->model=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
+                switch(parametricModel.getType()){
+                    case storm::models::ModelType::Dtmc:
+                        this->model=std::make_shared<storm::models::sparse::Dtmc<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
+                        break;
+                    case storm::models::ModelType::Mdp:
+                        this->model=std::make_shared<storm::models::sparse::Mdp<ConstantType>>(std::move(probabilityMatrix), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
+                        break;
+                    default:
+                        STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Tried to build a sampling model for an unsupported model type");
+                }
 
                 //translate the matrixEntryToEvalTableMapping into the actual probability mapping
                 auto sampleModelEntry = this->model->getTransitionMatrix().begin();
@@ -83,10 +97,15 @@ namespace storm {
                                                                                                                     std::vector<TableEntry*>& matrixEntryToEvalTableMapping,
                                                                                                                     TableEntry* constantEntry) {
                 // Run through the rows of the original model and obtain the probability evaluation table, a matrix with dummy entries and the mapping between the two.
-                storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getNumberOfStates(), parametricModel.getNumberOfStates(), parametricModel.getTransitionMatrix().getEntryCount());
+                bool addRowGroups = parametricModel.getTransitionMatrix().hasNontrivialRowGrouping();
+                auto currRowGroup = parametricModel.getTransitionMatrix().getRowGroupIndices().begin();
+                storm::storage::SparseMatrixBuilder<ConstantType> matrixBuilder(parametricModel.getTransitionMatrix().getRowCount(), parametricModel.getTransitionMatrix().getColumnCount(), parametricModel.getTransitionMatrix().getEntryCount(), addRowGroups);
                 matrixEntryToEvalTableMapping.reserve(parametricModel.getTransitionMatrix().getEntryCount());
                 std::size_t numOfNonConstEntries=0;
                 for(typename storm::storage::SparseMatrix<ParametricType>::index_type row=0; row < parametricModel.getTransitionMatrix().getRowCount(); ++row ){
+                    if(addRowGroups && row==*currRowGroup){
+                        matrixBuilder.newRowGroup(row);
+                    }
                     ConstantType dummyEntry=storm::utility::one<ConstantType>();
                     for(auto const& entry : parametricModel.getTransitionMatrix().getRow(row)){
                         if(storm::utility::isConstant(entry.getValue())){
@@ -133,7 +152,7 @@ namespace storm {
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> const& SamplingModel<ParametricSparseModelType, ConstantType>::getModel() const {
+            std::shared_ptr<storm::models::sparse::Model<ConstantType>> const& SamplingModel<ParametricSparseModelType, ConstantType>::getModel() const {
                 return this->model;
             }
 
@@ -162,20 +181,31 @@ namespace storm {
 
             template<typename ParametricSparseModelType, typename ConstantType>
             std::vector<ConstantType> const& SamplingModel<ParametricSparseModelType, ConstantType>::computeValues() {
-                storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<ConstantType>> modelChecker(*this->model);
+                std::unique_ptr<storm::modelchecker::AbstractModelChecker> modelChecker;
+                switch(this->getModel()->getType()){
+                    case storm::models::ModelType::Dtmc:
+                        modelChecker = std::make_unique<storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<ConstantType>>>(*this->model->template as<storm::models::sparse::Dtmc<ConstantType>>());
+                        break;
+                    case storm::models::ModelType::Mdp:
+                        modelChecker = std::make_unique<storm::modelchecker::SparseMdpPrctlModelChecker<storm::models::sparse::Mdp<ConstantType>>>(*this->model->template as<storm::models::sparse::Mdp<ConstantType>>());
+                        break;
+                    default:
+                        STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Tried to build a sampling model for an unsupported model type");
+                }
                 std::unique_ptr<storm::modelchecker::CheckResult> resultPtr;
-                //perform model checking on the dtmc
+                //perform model checking
                 if(this->computeRewards){
-                    resultPtr = modelChecker.computeReachabilityRewards(this->formula->asReachabilityRewardFormula());
+                    resultPtr = modelChecker->computeReachabilityRewards(this->formula->asRewardOperatorFormula().getSubformula().asReachabilityRewardFormula());
                 }
                 else {
-                    resultPtr = modelChecker.computeEventuallyProbabilities(this->formula->asEventuallyFormula());
+                    resultPtr = modelChecker->computeEventuallyProbabilities(this->formula->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula());
                 }
                 return resultPtr->asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
             }
         
 #ifdef STORM_HAVE_CARL
             template class SamplingModel<storm::models::sparse::Dtmc<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
+            template class SamplingModel<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
 #endif
         } //namespace region
     }
diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h
index 5b1464f91..26bdce429 100644
--- a/src/modelchecker/region/SamplingModel.h
+++ b/src/modelchecker/region/SamplingModel.h
@@ -14,7 +14,7 @@
 #include "src/utility/region.h"
 
 #include "src/logic/Formulas.h"
-#include "src/models/sparse/Dtmc.h"
+#include "src/models/sparse/Model.h"
 #include "src/storage/SparseMatrix.h"
 
 namespace storm {
@@ -26,14 +26,17 @@ namespace storm {
                 typedef typename ParametricSparseModelType::ValueType ParametricType;
                 typedef typename storm::utility::region::VariableType<ParametricType> VariableType;
                 typedef typename storm::utility::region::CoefficientType<ParametricType> CoefficientType;
-
-                SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::Formula> formula);
+                
+                /*!
+                 * Creates a sampling model.
+                 */
+                SamplingModel(ParametricSparseModelType const& parametricModel, std::shared_ptr<storm::logic::OperatorFormula> formula);
                 virtual ~SamplingModel();
 
                 /*!
                  * returns the underlying model
                  */
-                std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> const& getModel() const;
+                std::shared_ptr<storm::models::sparse::Model<ConstantType>> const& getModel() const;
 
                 /*!
                  * Instantiates the underlying model according to the given point
@@ -55,9 +58,9 @@ namespace storm {
                 void initializeRewards(ParametricSparseModelType const& parametricModel, boost::optional<std::vector<ConstantType>>& stateRewards, std::vector<TableEntry*>& rewardEntryToEvalTableMapping, TableEntry* constantEntry);
 
                 //The model with which we work
-                std::shared_ptr<storm::models::sparse::Dtmc<ConstantType>> model;
+                std::shared_ptr<storm::models::sparse::Model<ConstantType>> model;
                 //The formula for which we will compute the values
-                std::shared_ptr<storm::logic::Formula> formula;
+                std::shared_ptr<storm::logic::OperatorFormula> formula;
                 //A flag that denotes whether we compute probabilities or rewards
                 bool computeRewards;
 
diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
index 8f53175f6..c9778c5ad 100644
--- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
@@ -8,6 +8,7 @@
 #include "src/logic/Formulas.h"
 #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
 #include "src/modelchecker/region/RegionCheckResult.h"
+#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/RegionSettings.h"
@@ -24,7 +25,7 @@
 #include "src/exceptions/InvalidSettingsException.h"
 #include "src/exceptions/NotImplementedException.h"
 #include "src/exceptions/UnexpectedException.h"
-#include "exceptions/NotSupportedException.h"
+#include "src/exceptions/NotSupportedException.h"
 
 
 namespace storm {
@@ -62,21 +63,18 @@ namespace storm {
                     if (reachabilityRewardFormula.getSubformula().isPropositionalFormula()) {
                         return true;
                     }
-                } else if (formula.isConditionalPathFormula()) {
-                    storm::logic::ConditionalPathFormula conditionalPathFormula = formula.asConditionalPathFormula();
-                    if (conditionalPathFormula.getLeftSubformula().isEventuallyFormula() && conditionalPathFormula.getRightSubformula().isEventuallyFormula()) {
-                        return this->canHandle(conditionalPathFormula.getLeftSubformula()) && this->canHandle(conditionalPathFormula.getRightSubformula());
-                    }
                 }
                  STORM_LOG_DEBUG("Region Model Checker could not handle (sub)formula " << formula);
                 return false;
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::Formula>& simpleFormula, bool& isApproximationApplicable, bool& isResultConstant){
+            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel,
+                                                                                                   std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula,
+                                                                                                   bool& isApproximationApplicable,
+                                                                                                   boost::optional<ConstantType>& constantResult){
                 STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidArgumentException, "Input model is required to have exactly one initial state.");
                 //Reset some data
-                this->isResultInfinity=false;
                 this->smtSolver=nullptr;
                 this->reachabilityFunction=nullptr;
 
@@ -85,10 +83,10 @@ namespace storm {
                 boost::optional<std::vector<ParametricType>> stateRewards;
                 if (this->isComputeRewards()) {
                     std::vector<ParametricType> stateRewardsAsVector;
-                    preprocessForRewards(maybeStates, targetStates, stateRewardsAsVector, isApproximationApplicable, isResultConstant);
+                    preprocessForRewards(maybeStates, targetStates, stateRewardsAsVector, isApproximationApplicable, constantResult);
                     stateRewards=std::move(stateRewardsAsVector);
                 } else {
-                    preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, isResultConstant);
+                    preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, constantResult);
                 }
                 // Determine the set of states that is reachable from the initial state without jumping over a target state.
                 storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates(this->getModel().getTransitionMatrix(), this->getModel().getInitialStates(), maybeStates, targetStates);
@@ -220,22 +218,26 @@ namespace storm {
                 // the final model
                 simpleModel = std::make_shared<ParametricSparseModelType>(matrixBuilder.build(), std::move(labeling), std::move(rewardModels), std::move(noChoiceLabeling));
                 // the corresponding formula
-                std::shared_ptr<storm::logic::Formula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
+                std::shared_ptr<storm::logic::AtomicLabelFormula> targetFormulaPtr(new storm::logic::AtomicLabelFormula("target"));
                 if(this->isComputeRewards()){
-                    simpleFormula = std::shared_ptr<storm::logic::Formula>(new storm::logic::ReachabilityRewardFormula(targetFormulaPtr));
+                    std::shared_ptr<storm::logic::ReachabilityRewardFormula> reachRewFormula(new storm::logic::ReachabilityRewardFormula(targetFormulaPtr));
+                    simpleFormula = std::shared_ptr<storm::logic::OperatorFormula>(new storm::logic::RewardOperatorFormula(boost::optional<std::string>(), this->getSpecifiedFormula()->getComparisonType(), this->getSpecifiedFormulaBound(), reachRewFormula));
                 } else {
-                    simpleFormula = std::shared_ptr<storm::logic::Formula>(new storm::logic::EventuallyFormula(targetFormulaPtr));
+                    std::shared_ptr<storm::logic::EventuallyFormula> eventuallyFormula(new storm::logic::EventuallyFormula(targetFormulaPtr));
+                    simpleFormula = std::shared_ptr<storm::logic::OperatorFormula>(new storm::logic::ProbabilityOperatorFormula(this->getSpecifiedFormula()->getComparisonType(), this->getSpecifiedFormulaBound(), eventuallyFormula));
                 }
                 //Check if the reachability function needs to be computed
-                if((isResultConstant && this->reachabilityFunction==nullptr) || //If the result is constant, we can already compute it now.
-                        (storm::settings::regionSettings().getSmtMode()==storm::settings::modules::RegionSettings::SmtMode::FUNCTION) || //The smt solver will need it
-                        (storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::EVALUATE)){ //The sampling will need it 
+                if((storm::settings::regionSettings().getSmtMode()==storm::settings::modules::RegionSettings::SmtMode::FUNCTION) || 
+                        (storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::EVALUATE)){
                     this->computeReachabilityFunction(*simpleModel);
                 }
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable /* = 0 */, bool& isResultConstant /* = 0 */) {
+            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForProbabilities(storm::storage::BitVector& maybeStates,
+                                                                                                                   storm::storage::BitVector& targetStates,
+                                                                                                                   bool& isApproximationApplicable,
+                                                                                                                   boost::optional<ConstantType>& constantResult) {
                 //Get Target States
                 storm::logic::AtomicLabelFormula const& labelFormula = this->getSpecifiedFormula()->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula().getSubformula().asAtomicLabelFormula();
                 storm::modelchecker::SparsePropositionalModelChecker<ParametricSparseModelType> modelChecker(this->getModel());
@@ -249,14 +251,14 @@ namespace storm {
                 if (!maybeStates.get(initialState)) {
                     STORM_LOG_WARN("The probability of the initial state is constant (zero or one)");
                     this->reachabilityFunction = std::make_shared<ParametricType>(statesWithProbability01.first.get(initialState) ? storm::utility::zero<ParametricType>() : storm::utility::one<ParametricType>());
-                    isResultConstant=true;
+                    constantResult = statesWithProbability01.first.get(initialState) ? storm::utility::zero<ConstantType>() : storm::utility::one<ConstantType>();
                     return; //nothing else to do...
                 }
                 //extend target states
                 targetStates=statesWithProbability01.second;
                 //check if approximation is applicable and whether the result is constant
                 isApproximationApplicable=true;
-                isResultConstant=true;
+                bool isResultConstant=true;
                 for (auto state=maybeStates.begin(); (state!=maybeStates.end()) && isApproximationApplicable; ++state) {
                     for(auto const& entry : this->getModel().getTransitionMatrix().getRow(*state)){
                         if(!storm::utility::isConstant(entry.getValue())){
@@ -268,12 +270,19 @@ namespace storm {
                         }
                     }
                 }
-                STORM_LOG_WARN_COND(!isResultConstant, "For the given property, the reachability Value is constant, i.e., independent of the region");
+                if(isResultConstant){
+                    STORM_LOG_WARN("For the given property, the reachability Value is constant, i.e., independent of the region");
+                    constantResult = storm::utility::region::convertNumber<ConstantType>(-1.0);
+                }
             }
 
 
             template<typename ParametricSparseModelType, typename ConstantType>
-            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForRewards(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, std::vector<ParametricType>& stateRewards, bool& isApproximationApplicable /* = true */, bool& isResultConstant /* = false */) {
+            void SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForRewards(storm::storage::BitVector& maybeStates,
+                                                                                                             storm::storage::BitVector& targetStates,
+                                                                                                             std::vector<ParametricType>& stateRewards,
+                                                                                                             bool& isApproximationApplicable,
+                                                                                                             boost::optional<ConstantType>& constantResult) {
                 //get the correct reward model
                 ParametricRewardModelType const* rewardModel;
                 if(this->getSpecifiedFormula()->asRewardOperatorFormula().hasRewardModelName()){
@@ -299,18 +308,15 @@ namespace storm {
                 storm::storage::sparse::state_type initialState = *this->getModel().getInitialStates().begin();
                 if (!maybeStates.get(initialState)) {
                     STORM_LOG_WARN("The expected reward of the initial state is constant (infinity or zero)");
-                    if(statesWithProbability1.get(initialState)){
-                        this->reachabilityFunction = std::make_shared<ParametricType>(storm::utility::zero<ParametricType>());
-                    } else {
-                        this->reachabilityFunction = std::make_shared<ParametricType>(storm::utility::one<ParametricType>());
-                        this->isResultInfinity=true;
-                    }
-                    isResultConstant=true;
+                    // Note: storm::utility::infinity<storm::RationalFunction> does not work at this moment.
+                    // In that case, we are going to throw in exception if the function is accessed (i.e. in getReachabilityFunction);
+                    this->reachabilityFunction = statesWithProbability1.get(initialState) ? std::make_shared<ParametricType>(storm::utility::zero<ParametricType>()) : nullptr;
+                    constantResult = statesWithProbability1.get(initialState) ? storm::utility::zero<ConstantType>() : storm::utility::infinity<ConstantType>();
                     return; //nothing else to do...
                 }
                  //check if approximation is applicable and whether the result is constant
                 isApproximationApplicable=true;
-                isResultConstant=true;
+                bool isResultConstant=true;
                 std::set<VariableType> rewardPars; //the set of parameters that occur on a reward function
                 std::set<VariableType> probPars;   //the set of parameters that occur on a probability function
                 for (auto state=maybeStates.begin(); state!=maybeStates.end() && isApproximationApplicable; ++state) {
@@ -357,7 +363,10 @@ namespace storm {
                         break;
                     }
                 }
-                STORM_LOG_WARN_COND(!isResultConstant, "For the given property, the reachability Value is constant, i.e., independent of the region");
+                if(isResultConstant){
+                    STORM_LOG_WARN("For the given property, the reachability Value is constant, i.e., independent of the region");
+                    constantResult = storm::utility::region::convertNumber<ConstantType>(-1.0);
+                }
             }
 
             template<typename ParametricSparseModelType, typename ConstantType>
@@ -400,7 +409,7 @@ namespace storm {
                 std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now();
                 this->timeApproxModelInstantiation += timeMDPBuildEnd-timeMDPBuildStart;
 
-                // Decide whether to prove allsat or allviolated. (Hence, there are 4 cases)
+                // Decide whether to prove allsat or allviolated. 
                 bool proveAllSat;
                 switch (region.getCheckResult()){
                     case RegionCheckResult::UNKNOWN: 
@@ -494,11 +503,11 @@ namespace storm {
                 if((storm::settings::regionSettings().getSampleMode()==storm::settings::modules::RegionSettings::SampleMode::EVALUATE) ||
                         (!storm::settings::regionSettings().doSample() && favorViaFunction)){
                     //evaluate the reachability function
-                    valueInBoundOfFormula = this->valueIsInBoundOfFormula(this->getReachabilityValue(point, true));
+                    valueInBoundOfFormula = this->valueIsInBoundOfFormula(this->evaluateReachabilityFunction(point));
                 }
                 else{
                     //instantiate the sampling model
-                    valueInBoundOfFormula = this->valueIsInBoundOfFormula(this->getReachabilityValue(point, false));
+                    valueInBoundOfFormula = this->valueIsInBoundOfFormula(this->getReachabilityValue(point));
                 }
 
                 if(valueInBoundOfFormula){
@@ -527,38 +536,21 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             std::shared_ptr<typename SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::ParametricType> const& SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::getReachabilityFunction() {
                 if(this->reachabilityFunction==nullptr){
+                    //Todo: remove workaround (infinity<storm::RationalNumber>() does not work)
+                    std::map<VariableType, CoefficientType> emptySubstitution;
+                    if(this->isResultConstant() && this->getReachabilityValue(emptySubstitution)==storm::utility::infinity<ConstantType>()){
+                        STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Requested the reachability function but it can not be represented (The value is infinity)");
+                        return this->reachabilityFunction;
+                    }
                     STORM_LOG_WARN("Reachability Function requested but it has not been computed when specifying the formula. Will compute it now.");
                     computeReachabilityFunction(*this->getSimpleModel());
                 }
+                STORM_LOG_THROW((!this->isResultConstant() || storm::utility::isConstant(*this->reachabilityFunction)), storm::exceptions::UnexpectedException, "The result was assumed to be constant but it isn't.");
                 return this->reachabilityFunction;
             }
 
-            template<typename ParametricSparseModelType, typename ConstantType>
-            ConstantType SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::getReachabilityValue(std::map<VariableType, CoefficientType> const& point, bool evaluateFunction) {
-                if(this->isResultConstant()){
-                    //Todo: remove workaround (infinity<storm::RationalNumber>() does not work)
-                    if(this->isResultInfinity){
-                        return storm::utility::infinity<ConstantType>();
-                    }
-                    STORM_LOG_THROW(storm::utility::isConstant(*getReachabilityFunction()), storm::exceptions::UnexpectedException, "The result was assumed to be constant but it isn't.");
-                    return storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(*getReachabilityFunction()));
-                }
-                if(evaluateFunction){
-                    return storm::utility::region::convertNumber<ConstantType>(this->evaluateReachabilityFunction(point));
-                } else {
-                    this->getSamplingModel()->instantiate(point);
-                    return storm::utility::region::convertNumber<ConstantType>(this->getSamplingModel()->computeValues()[*this->getSamplingModel()->getModel()->getInitialStates().begin()]);
-                }
-            }
-            
             template<typename ParametricSparseModelType, typename ConstantType>
             typename SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::CoefficientType SparseDtmcRegionModelChecker<ParametricSparseModelType, ConstantType>::evaluateReachabilityFunction(std::map<VariableType, CoefficientType> const& point) {
-                if(this->isResultConstant()){
-                    //Todo: remove workaround (infinity<storm::RationalNumber>() does not work)
-                    STORM_LOG_THROW(!this->isResultInfinity, storm::exceptions::NotSupportedException, "Requested the reachability value but can not represent it as a Coefficient type (The value is infinity)");
-                    STORM_LOG_THROW(storm::utility::isConstant(*getReachabilityFunction()), storm::exceptions::UnexpectedException, "The result was assumed to be constant but it isn't.");
-                    return storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(*getReachabilityFunction()));
-                }
                 return storm::utility::region::evaluateFunction(*getReachabilityFunction(), point);
             }
 
@@ -678,7 +670,7 @@ namespace storm {
                 // Hence: If  f(x) > p is unsat, the property is satisfied for all parameters.
 
                 storm::logic::ComparisonType proveAllSatRel; //the relation from the property needs to be inverted
-                switch (this->getSpecifiedFormulaCompType()) {
+                switch (this->getSpecifiedFormula()->getComparisonType()) {
                     case storm::logic::ComparisonType::Greater:
                         proveAllSatRel=storm::logic::ComparisonType::LessEqual;
                         break;
@@ -705,7 +697,7 @@ namespace storm {
                 //Property:    P<=p [ F 'target' ] holds iff...
                 // f(x)         <= p
                 // Hence: If f(x)  <= p is unsat, the property is violated for all parameters. 
-                storm::logic::ComparisonType proveAllViolatedRel = this->getSpecifiedFormulaCompType();
+                storm::logic::ComparisonType proveAllViolatedRel = this->getSpecifiedFormula()->getComparisonType();
                 storm::utility::region::addGuardedConstraintToSmtSolver(this->smtSolver, proveAllViolatedVar, *getReachabilityFunction(), proveAllViolatedRel, bound);          
             }
 
@@ -714,7 +706,6 @@ namespace storm {
 #ifdef STORM_HAVE_CARL
             template class SparseDtmcRegionModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
 #endif
-            //note: for other template instantiations, add rules for the typedefs of VariableType and CoefficientType in utility/regions.h
         } // namespace region 
     } // namespace modelchecker
 } // namespace storm
diff --git a/src/modelchecker/region/SparseDtmcRegionModelChecker.h b/src/modelchecker/region/SparseDtmcRegionModelChecker.h
index f3b993137..a36b22b80 100644
--- a/src/modelchecker/region/SparseDtmcRegionModelChecker.h
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.h
@@ -36,16 +36,6 @@ namespace storm {
                  * If it is not yet available, it is computed.
                  */
                 std::shared_ptr<ParametricType> const& getReachabilityFunction();
-
-                /*!
-                 * Returns the reachability Value at the specified point. 
-                 * The given flag decides whether to initialize a sampling model or to evaluate a reachability function.
-                 * Might invoke sampling model initialization or the computation of the reachability function (if these are not available yet)
-                 * 
-                 * @param point The point (i.e. parameter evaluation) at which to compute the reachability value.
-                 * @param evaluateFunction If set, the reachability function is evaluated. Otherwise, the sampling model is instantiated.
-                 */
-                virtual ConstantType getReachabilityValue(std::map<VariableType, CoefficientType>const& point, bool evaluateFunction=false);
                 
                 /*!
                  * Evaluates the reachability function with the given substitution.
@@ -57,41 +47,41 @@ namespace storm {
             protected:
                 
                 /*!
-                 * Checks whether the approximation technique is applicable and whether the model checking result is independent of parameters (i.e., constant)
-                 * 
+                 * Checks whether the approximation technique is applicable and whether the model checking result is independent of parameters (i.e., constant).
+                 * In the latter case, the given parameter is set either to the correct result or -1
                  * Computes a model with a single target and at most one sink state.
                  * Eliminates all states for which the outgoing transitions are constant.
                  * If rewards are relevant, transition rewards are transformed to state rewards
                  * 
                  * @note this->specifiedFormula and this->computeRewards has to be set accordingly before calling this function
                  */
-                virtual void preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::Formula>& simpleFormula, bool& isApproximationApplicable, bool& isResultConstant);
+                virtual void preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult);
 
                 private:
                 /*!
                  * Does some sanity checks and preprocessing steps on the currently specified model and 
                  * reachability probability formula, i.e., 
-                 * * Sets some formula data and that we do not compute rewards
                  * * Computes maybeStates and targetStates
-                 * * Sets the flags that state whether the result is constant and approximation is applicable
+                 * * Sets whether approximation is applicable
+                 * * If the result is constant, it is already computed or set to -1
                  * 
                  * @note The returned set of target states also includes states where an 'actual' target state is reached with probability 1
                  * 
                  */
-                void preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable, bool& isResultConstant);
+                void preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult);
 
                 /*!
                  * Does some sanity checks and preprocessing steps on the currently specified model and 
                  * reachability reward formula, i.e.
-                 * * Sets some formula data and that we do compute rewards
                  * * Computes maybeStates, targetStates
                  * * Computes a new stateReward vector that considers state+transition rewards of the original model. (in a sense that we can abstract away from transition rewards)
-                 * * Sets the flags that state whether the result is constant and approximation is applicable
+                 * * Sets whether approximation is applicable
+                 * * If the result is constant, it is already computed or set to -1
                  * 
                  * @note stateRewards.size will equal to maybeStates.numberOfSetBits
                  * 
                  */
-                void preprocessForRewards(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, std::vector<ParametricType>& stateRewards, bool& isApproximationApplicable, bool& isResultConstant);
+                void preprocessForRewards(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, std::vector<ParametricType>& stateRewards, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult);
 
                 /*!
                  * Computes the reachability function via state elimination
@@ -139,8 +129,6 @@ namespace storm {
 
                 // The  function for the reachability probability (or: reachability reward) in the initial state 
                 std::shared_ptr<ParametricType> reachabilityFunction;
-                // workaround to represent that the result is infinity (utility::infinity<storm::RationalFunction>() does not work at this moment)
-                bool isResultInfinity;
                 
                 // Instance of an elimination model checker to access its functions
                 storm::modelchecker::SparseDtmcEliminationModelChecker<storm::models::sparse::Dtmc<ParametricType>> eliminationModelChecker;
diff --git a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp
new file mode 100644
index 000000000..b2f0966f6
--- /dev/null
+++ b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp
@@ -0,0 +1,254 @@
+#include "src/modelchecker/region/SparseMdpRegionModelChecker.h"
+
+#include <chrono>
+#include <memory>
+#include <boost/optional.hpp>
+
+#include "src/adapters/CarlAdapter.h"
+#include "src/logic/Formulas.h"
+#include "src/modelchecker/results/ExplicitQualitativeCheckResult.h"
+#include "src/modelchecker/region/RegionCheckResult.h"
+#include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
+#include "src/models/sparse/StandardRewardModel.h"
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/RegionSettings.h"
+#include "src/solver/OptimizationDirection.h"
+#include "src/storage/sparse/StateType.h"
+#include "src/utility/constants.h"
+#include "src/utility/graph.h"
+#include "src/utility/macros.h"
+#include "src/utility/vector.h"
+
+#include "src/exceptions/InvalidArgumentException.h"
+#include "src/exceptions/InvalidPropertyException.h"
+#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidSettingsException.h"
+#include "src/exceptions/NotImplementedException.h"
+#include "src/exceptions/UnexpectedException.h"
+#include "src/exceptions/NotSupportedException.h"
+
+
+namespace storm {
+    namespace modelchecker {
+        namespace region {
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::SparseMdpRegionModelChecker(ParametricSparseModelType const& model) : 
+                    AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>(model){
+                //intentionally left empty
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::~SparseMdpRegionModelChecker(){
+                //intentionally left empty
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::canHandle(storm::logic::Formula const& formula) const {
+                 //for simplicity we only support state formulas with eventually (e.g. P<0.5 [ F "target" ])
+                if (formula.isProbabilityOperatorFormula()) {
+                    storm::logic::ProbabilityOperatorFormula const& probabilityOperatorFormula = formula.asProbabilityOperatorFormula();
+                    return probabilityOperatorFormula.hasBound() &&  this->canHandle(probabilityOperatorFormula.getSubformula());
+                //} else if (formula.isRewardOperatorFormula()) {
+                //    storm::logic::RewardOperatorFormula const& rewardOperatorFormula = formula.asRewardOperatorFormula();
+                //    return rewardOperatorFormula.hasBound() && this->canHandle(rewardOperatorFormula.getSubformula());
+                } else if (formula.isEventuallyFormula()) {
+                    storm::logic::EventuallyFormula const& eventuallyFormula = formula.asEventuallyFormula();
+                    if (eventuallyFormula.getSubformula().isPropositionalFormula()) {
+                        return true;
+                    }
+                } else if (formula.isReachabilityRewardFormula()) {
+                    storm::logic::ReachabilityRewardFormula reachabilityRewardFormula = formula.asReachabilityRewardFormula();
+                    if (reachabilityRewardFormula.getSubformula().isPropositionalFormula()) {
+                        return true;
+                    }
+                }
+                 STORM_LOG_DEBUG("Region Model Checker could not handle (sub)formula " << formula);
+                return false;
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            void SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel,
+                                                                                                  std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula,
+                                                                                                  bool& isApproximationApplicable,
+                                                                                                  boost::optional<ConstantType>& constantResult){
+                STORM_LOG_THROW(this->getModel().getInitialStates().getNumberOfSetBits() == 1, storm::exceptions::InvalidArgumentException, "Input model is required to have exactly one initial state.");
+                storm::storage::BitVector maybeStates, targetStates;
+                preprocessForProbabilities(maybeStates, targetStates, isApproximationApplicable, constantResult);
+                //TODO: Actually get a more simple model. This makes a deep copy of the model...
+                simpleModel = std::make_shared<ParametricSparseModelType>(this->getModel());
+                simpleFormula = this->getSpecifiedFormula();
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            void SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult) {
+                //Get Target States
+                storm::logic::AtomicLabelFormula const& labelFormula = this->getSpecifiedFormula()->asProbabilityOperatorFormula().getSubformula().asEventuallyFormula().getSubformula().asAtomicLabelFormula();
+                storm::modelchecker::SparsePropositionalModelChecker<ParametricSparseModelType> modelChecker(this->getModel());
+                std::unique_ptr<CheckResult> targetStatesResultPtr = modelChecker.checkAtomicLabelFormula(labelFormula);
+                targetStates = std::move(targetStatesResultPtr->asExplicitQualitativeCheckResult().getTruthValuesVector());
+                
+                //maybeStates: Compute the subset of states that have a probability of 0 or 1, respectively and reduce the considered states accordingly.
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01;
+                if (this->specifiedFormulaHasUpperBound()){
+                    statesWithProbability01 = storm::utility::graph::performProb01Max(this->getModel(), storm::storage::BitVector(this->getModel().getNumberOfStates(),true), targetStates);
+                } else {
+                    statesWithProbability01 = storm::utility::graph::performProb01Min(this->getModel(), storm::storage::BitVector(this->getModel().getNumberOfStates(),true), targetStates);
+                }
+                maybeStates = ~(statesWithProbability01.first | statesWithProbability01.second);
+                // If the initial state is known to have either probability 0 or 1, we can directly set the reachProbFunction.
+                storm::storage::sparse::state_type initialState = *this->getModel().getInitialStates().begin();
+                if (!maybeStates.get(initialState)) {
+                    STORM_LOG_WARN("The probability of the initial state is constant (zero or one)");
+                    constantResult = statesWithProbability01.first.get(initialState) ? storm::utility::zero<ConstantType>() : storm::utility::one<ConstantType>();
+                    return; //nothing else to do...
+                }
+                //extend target states
+                targetStates=statesWithProbability01.second;
+                //check if approximation is applicable and whether the result is constant
+                isApproximationApplicable=true;
+                bool isResultConstant=true;
+                for (auto state=maybeStates.begin(); (state!=maybeStates.end()) && isApproximationApplicable; ++state) {
+                    for(auto const& entry : this->getModel().getTransitionMatrix().getRowGroup(*state)){
+                        if(!storm::utility::isConstant(entry.getValue())){
+                            isResultConstant=false;
+                            if(!storm::utility::region::functionIsLinear(entry.getValue())){
+                                isApproximationApplicable=false;
+                                break;
+                            }
+                        }
+                    }
+                }
+                if(isResultConstant){
+                    STORM_LOG_WARN("For the given property, the reachability Value is constant, i.e., independent of the region");
+                    constantResult = storm::utility::region::convertNumber<ConstantType>(-1.0); //-1 denotes that the result is constant but not yet computed
+                }
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::checkApproximativeValues(ParameterRegion<ParametricType>& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds) {
+                std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now();
+                STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "The approximation model for mdps has not been implemented yet");
+                this->getApproximationModel()->instantiate(region);
+                std::chrono::high_resolution_clock::time_point timeMDPBuildEnd = std::chrono::high_resolution_clock::now();
+                this->timeApproxModelInstantiation += timeMDPBuildEnd-timeMDPBuildStart;
+
+                // Decide whether to prove allsat or allviolated. 
+                bool proveAllSat;
+                switch (region.getCheckResult()){
+                    case RegionCheckResult::UNKNOWN: 
+                        switch(storm::settings::regionSettings().getApproxMode()){
+                            case storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST:
+                                //Sample a single point to know whether we should try to prove ALLSAT or ALLVIOLATED
+                                checkPoint(region,region.getSomePoint(), false);
+                                proveAllSat= (region.getCheckResult()==RegionCheckResult::EXISTSSAT);
+                                break;
+                            case storm::settings::modules::RegionSettings::ApproxMode::GUESSALLSAT:
+                                proveAllSat=true;
+                                break;
+                            case storm::settings::modules::RegionSettings::ApproxMode::GUESSALLVIOLATED:
+                                proveAllSat=false;
+                                break;
+                            default:
+                                STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The specified approxmode is not supported");
+                        }
+                        break;
+                    case RegionCheckResult::ALLSAT:
+                         STORM_LOG_WARN("The checkresult of the current region should not be conclusive (ALLSAT)");
+                         //Intentionally no break;
+                    case RegionCheckResult::EXISTSSAT:
+                        proveAllSat=true;
+                        break;
+                    case RegionCheckResult::ALLVIOLATED:
+                         STORM_LOG_WARN("The checkresult of the current region should not be conclusive (ALLViolated)");
+                         //Intentionally no break;
+                    case RegionCheckResult::EXISTSVIOLATED:
+                        proveAllSat=false;
+                        break;
+                    default:
+                         STORM_LOG_WARN("The checkresult of the current region should not be conclusive, i.e. it should be either EXISTSSAT or EXISTSVIOLATED or UNKNOWN in order to apply approximative values");
+                         proveAllSat=true;
+                }
+
+                bool formulaSatisfied;
+                if((this->specifiedFormulaHasUpperBound() && proveAllSat) || (!this->specifiedFormulaHasUpperBound() && !proveAllSat)){
+                    //these are the cases in which we need to compute upper bounds
+                    upperBounds = this->getApproximationModel()->computeValues(storm::solver::OptimizationDirection::Maximize);
+                    lowerBounds = std::vector<ConstantType>();
+                    formulaSatisfied = this->valueIsInBoundOfFormula(upperBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
+                }
+                else{
+                    //for the remaining cases we compute lower bounds
+                    lowerBounds = this->getApproximationModel()->computeValues(storm::solver::OptimizationDirection::Minimize);
+                    upperBounds = std::vector<ConstantType>();
+                    formulaSatisfied = this->valueIsInBoundOfFormula(lowerBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
+                }
+
+                //check if approximation was conclusive
+                if(proveAllSat && formulaSatisfied){
+                    region.setCheckResult(RegionCheckResult::ALLSAT);
+                    return true;
+                }
+                if(!proveAllSat && !formulaSatisfied){
+                    region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                    return true;
+                }
+
+                if(region.getCheckResult()==RegionCheckResult::UNKNOWN){
+                    //In this case, it makes sense to try to prove the contrary statement
+                    proveAllSat=!proveAllSat;
+
+                    if(lowerBounds.empty()){
+                        lowerBounds = this->getApproximationModel()->computeValues(storm::solver::OptimizationDirection::Minimize);
+                        formulaSatisfied=this->valueIsInBoundOfFormula(lowerBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
+                    }
+                    else{
+                        upperBounds = this->getApproximationModel()->computeValues(storm::solver::OptimizationDirection::Maximize);
+                        formulaSatisfied=this->valueIsInBoundOfFormula(upperBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
+                    }
+
+                    //check if approximation was conclusive
+                    if(proveAllSat && formulaSatisfied){
+                        region.setCheckResult(RegionCheckResult::ALLSAT);
+                        return true;
+                    }
+                    if(!proveAllSat && !formulaSatisfied){
+                        region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                        return true;
+                    }
+                }
+                //if we reach this point than the result is still inconclusive.
+                return false;            
+            }
+
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool SparseMdpRegionModelChecker<ParametricSparseModelType, ConstantType>::checkPoint(ParameterRegion<ParametricType>& region, std::map<VariableType, CoefficientType>const& point, bool favorViaFunction) {
+                if(this->valueIsInBoundOfFormula(this->getReachabilityValue(point))){
+                    if (region.getCheckResult()!=RegionCheckResult::EXISTSSAT){
+                        region.setSatPoint(point);
+                        if(region.getCheckResult()==RegionCheckResult::EXISTSVIOLATED){
+                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                            return true;
+                        }
+                        region.setCheckResult(RegionCheckResult::EXISTSSAT);
+                    }
+                }
+                else{
+                    if (region.getCheckResult()!=RegionCheckResult::EXISTSVIOLATED){
+                        region.setViolatedPoint(point);
+                        if(region.getCheckResult()==RegionCheckResult::EXISTSSAT){
+                            region.setCheckResult(RegionCheckResult::EXISTSBOTH);
+                            return true;
+                        }
+                        region.setCheckResult(RegionCheckResult::EXISTSVIOLATED);
+                    }
+                }
+                return false;
+            }
+
+#ifdef STORM_HAVE_CARL
+            template class SparseMdpRegionModelChecker<storm::models::sparse::Mdp<storm::RationalFunction, storm::models::sparse::StandardRewardModel<storm::RationalFunction>>, double>;
+#endif
+        } // namespace region 
+    } // namespace modelchecker
+} // namespace storm
diff --git a/src/modelchecker/region/SparseMdpRegionModelChecker.h b/src/modelchecker/region/SparseMdpRegionModelChecker.h
new file mode 100644
index 000000000..ce1f9be4f
--- /dev/null
+++ b/src/modelchecker/region/SparseMdpRegionModelChecker.h
@@ -0,0 +1,88 @@
+#ifndef STORM_MODELCHECKER_REACHABILITY_SPARSEMDPREGIONMODELCHECKER_H_
+#define STORM_MODELCHECKER_REACHABILITY_SPARSEMDPREGIONMODELCHECKER_H_
+
+#include "src/modelchecker/region/AbstractSparseRegionModelChecker.h"
+
+#include "src/models/sparse/StandardRewardModel.h"
+#include "src/models/sparse/Mdp.h"
+#include "src/utility/region.h"
+
+namespace storm {
+    namespace modelchecker {
+        namespace region {
+            template<typename ParametricSparseModelType, typename ConstantType>
+            class SparseMdpRegionModelChecker : public AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType> {
+            public:
+
+                typedef typename ParametricSparseModelType::ValueType ParametricType;
+                typedef typename ParametricSparseModelType::RewardModelType ParametricRewardModelType;
+                typedef typename storm::utility::region::VariableType<ParametricType> VariableType;
+                typedef typename storm::utility::region::CoefficientType<ParametricType> CoefficientType;
+
+                explicit SparseMdpRegionModelChecker(ParametricSparseModelType const& model);
+
+                virtual ~SparseMdpRegionModelChecker();
+
+                /*!
+                 * Checks if the given formula can be handled by This region model checker
+                 * @param formula the formula to be checked
+                 */
+                virtual bool canHandle(storm::logic::Formula const& formula) const;
+
+            protected:
+                
+                /*!
+                 * Checks whether the approximation technique is applicable and whether the model checking result is independent of parameters (i.e., constant).
+                 * In the latter case, the given parameter is set either to the correct result or -1
+                 * Computes a model with a single target and at most one sink state.
+                 * Eliminates all states for which the outgoing transitions are constant.
+                 * If rewards are relevant, transition rewards are transformed to state rewards
+                 * 
+                 * @note this->specifiedFormula and this->computeRewards has to be set accordingly before calling this function
+                 */
+                virtual void preprocess(std::shared_ptr<ParametricSparseModelType>& simpleModel, std::shared_ptr<storm::logic::OperatorFormula>& simpleFormula, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult);
+
+                private:
+                /*!
+                 * Does some sanity checks and preprocessing steps on the currently specified model and 
+                 * reachability probability formula, i.e., 
+                 * * Computes maybeStates and targetStates
+                 * * Sets whether approximation is applicable
+                 * * If the result is constant, it is already computed or set to -1
+                 * 
+                 * @note The returned set of target states also includes states where an 'actual' target state is reached with probability 1
+                 * 
+                 */
+                void preprocessForProbabilities(storm::storage::BitVector& maybeStates, storm::storage::BitVector& targetStates, bool& isApproximationApplicable, boost::optional<ConstantType>& constantResult);
+
+                /*!
+                 * Instantiates the approximation model to compute bounds on the maximal/minimal reachability probability (or reachability reward).
+                 * If the current region result is EXISTSSAT (or EXISTSVIOLATED), then this function tries to prove ALLSAT (or ALLVIOLATED).
+                 * If this succeeded, then the region check result is changed accordingly.
+                 * If the current region result is UNKNOWN, then this function first tries to prove ALLSAT and if that failed, it tries to prove ALLVIOLATED.
+                 * In any case, the computed bounds are written to the given lowerBounds/upperBounds.
+                 * However, if only the lowerBounds (or upperBounds) have been computed, the other vector is set to a vector of size 0.
+                 * True is returned iff either ALLSAT or ALLVIOLATED could be proved.
+                 */
+                virtual bool checkApproximativeValues(ParameterRegion<ParametricType>& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds); 
+
+                /*!
+                 * Checks the value of the function at the given sampling point.
+                 * May set the satPoint and violatedPoint of the regions if thy are not yet specified and such point is given.
+                 * Also changes the regioncheckresult of the region to EXISTSSAT, EXISTSVIOLATED, or EXISTSBOTH
+                 * 
+                 * @param favorViaFunction if not stated otherwise (e.g. in the settings), the sampling will be done via the
+                 *                          reachabilityFunction if this flag is true. If the flag is false, sampling will be 
+                 *                          done via instantiation of the samplingmodel. Note that this argument is ignored,
+                 *                          unless sampling has been turned of in the settings
+                 * 
+                 * @return true if an violated point as well as a sat point has been found, i.e., the check result is changed to EXISTSOTH
+                 */
+                virtual bool checkPoint(ParameterRegion<ParametricType>& region, std::map<VariableType, CoefficientType>const& point, bool favorViaFunction=false);
+
+            };
+        } //namespace region
+    } // namespace modelchecker
+} // namespace storm
+
+#endif /* STORM_MODELCHECKER_REACHABILITY_SPARSEMDPREGIONMODELCHECKER_H_ */
diff --git a/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp b/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
index 1034811de..e5f3a213a 100644
--- a/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
@@ -46,18 +46,18 @@ TEST(SparseDtmcRegionModelCheckerTest, Brp_Prob) {
     auto exBothRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.4<=pL<=0.65,0.75<=pK<=0.95");
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.1<=pL<=0.9,0.2<=pK<=0.5");
 
-    EXPECT_NEAR(0.8369631407, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.8369631407, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.0476784174, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.0476784174, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.9987948367, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.9987948367, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6020480995, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6020480995, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(1.0000000000, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(1.0000000000, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.9456084185, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.9456084185, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(0.8369631407, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.8369631407, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.0476784174, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.0476784174, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9987948367, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9987948367, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6020480995, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6020480995, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.0000000000, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.0000000000, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9456084185, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.9456084185, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -116,22 +116,22 @@ TEST(SparseDtmcRegionModelCheckerTest, Brp_Rew) {
     auto exBothHardRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.5<=pK<=0.75,0.3<=TOMsg<=0.4"); //this region has a local maximum!
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.1<=pK<=0.3,0.2<=TOMsg<=0.3");
 
-    EXPECT_NEAR(4.367791292, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(4.367791292, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(3.044795147, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(3.044795147, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(3.182535759, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(3.182535759, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(2.609602197, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(2.609602197, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(1.842551039, modelchecker.getReachabilityValue(exBothHardRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(1.842551039, modelchecker.getReachabilityValue(exBothHardRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(2.453500364, modelchecker.getReachabilityValue(exBothHardRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(2.453500364, modelchecker.getReachabilityValue(exBothHardRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6721974438, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.6721974438, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(1.308324558, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(1.308324558, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(4.367791292, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.367791292, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.044795147, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.044795147, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.182535759, modelchecker.getReachabilityValue(exBothRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(3.182535759, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.609602197, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.609602197, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.842551039, modelchecker.getReachabilityValue(exBothHardRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.842551039, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothHardRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.453500364, modelchecker.getReachabilityValue(exBothHardRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(2.453500364, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothHardRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6721974438, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6721974438, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.308324558, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(1.308324558, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -206,8 +206,7 @@ TEST(SparseDtmcRegionModelCheckerTest, Brp_Rew_Infty) {
     //start testing
     auto allSatRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("");
 
-    EXPECT_EQ(storm::utility::infinity<double>(), modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false)); //instantiation of sampling model
-    EXPECT_EQ(storm::utility::infinity<double>(), modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true)); //instantiation of sampling model
+    EXPECT_EQ(storm::utility::infinity<double>(), modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()));
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -255,12 +254,12 @@ TEST(SparseDtmcRegionModelCheckerTest, Brp_Rew_4Par) {
     auto exBothRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.1<=pK<=0.7,0.2<=pL<=0.8,0.15<=TOMsg<=0.65,0.3<=TOAck<=0.9");
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.1<=pK<=0.4,0.2<=pL<=0.3,0.15<=TOMsg<=0.3,0.1<=TOAck<=0.2");
 
-    EXPECT_NEAR(4.834779705, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(4.834779705, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(4.674651623, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(4.674651623, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.4467496536, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.4467496536, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(4.834779705, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.834779705, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.674651623, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(4.674651623, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.4467496536, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.4467496536, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -319,18 +318,18 @@ TEST(SparseDtmcRegionModelCheckerTest, Crowds_Prob) {
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.8<=PF<=0.95,0.2<=badC<=0.2");
     auto allVioHardRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.8<=PF<=0.95,0.2<=badC<=0.9");
 
-    EXPECT_NEAR(0.1734086422, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.1734086422, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.47178, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.47178, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.5095205203, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.5095205203, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6819701472, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.6819701472, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.999895897, modelchecker.getReachabilityValue(allVioHardRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.999895897, modelchecker.getReachabilityValue(allVioHardRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(0.1734086422, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.1734086422, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.47178, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.47178, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7085157883, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.5095205203, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.5095205203, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6819701472, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6819701472, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.999895897, modelchecker.getReachabilityValue(allVioHardRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.999895897, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioHardRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -406,14 +405,14 @@ TEST(SparseDtmcRegionModelCheckerTest, Crowds_Prob_1Par) {
     auto exBothRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.8<=PF<=0.9");
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.01<=PF<=0.8");
 
-    EXPECT_NEAR(0.8430128158, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.8430128158, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.7731321947, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.7731321947, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.4732302663, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.4732302663, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), false), storm::settings::generalSettings().getPrecision()); //instantiation of sampling model
-    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(0.8430128158, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.8430128158, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7731321947, modelchecker.getReachabilityValue(exBothRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7731321947, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(exBothRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.4732302663, modelchecker.getReachabilityValue(allVioRegion.getLowerBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.4732302663, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7085157883, modelchecker.getReachabilityValue(allVioRegion.getUpperBounds()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.7085157883, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allVioRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -465,14 +464,13 @@ TEST(SparseDtmcRegionModelCheckerTest, Crowds_Prob_Const) {
     storm::modelchecker::region::SparseDtmcRegionModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>, double> modelchecker(*dtmc);
     ASSERT_TRUE(modelchecker.canHandle(*formulas[0]));
     modelchecker.specifyFormula(formulas[0]);
-    
     //start testing
     auto allSatRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("");
     
-    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), false),  storm::settings::generalSettings().getPrecision()); //evaluation of function
-    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds(), true),  storm::settings::generalSettings().getPrecision()); //evaluation of function
+    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getUpperBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6119660237, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getUpperBounds())),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6119660237, modelchecker.getReachabilityValue(allSatRegion.getLowerBounds()),  storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.6119660237, storm::utility::region::convertNumber<double>(modelchecker.evaluateReachabilityFunction(allSatRegion.getLowerBounds())),  storm::settings::generalSettings().getPrecision());
     
     //test approximative method
     storm::settings::mutableRegionSettings().modifyModes(storm::settings::modules::RegionSettings::ApproxMode::TESTFIRST, storm::settings::modules::RegionSettings::SampleMode::INSTANTIATE, storm::settings::modules::RegionSettings::SmtMode::OFF);
@@ -491,6 +489,7 @@ TEST(SparseDtmcRegionModelCheckerTest, Crowds_Prob_Const) {
     modelchecker.checkRegion(allSatRegionSmt);
 //smt    EXPECT_EQ((storm::modelchecker::region::RegionCheckResult::ALLSAT), allSatRegionSmt.getCheckResult());
 
+    std::cout << "End" << std::endl;
     storm::settings::mutableRegionSettings().resetModes();
 }