diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
index 66ecc60f8..f6c52095f 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
@@ -204,7 +204,7 @@ namespace storm {
 
                 if(this->isResultConstant()){
                     STORM_LOG_DEBUG("Checking a region although the result is constant, i.e., independent of the region. This makes sense none.");
-                    if(this->valueIsInBoundOfFormula(this->getReachabilityValue(region.getSomePoint()))){
+                    if(this->checkFormulaOnSamplingPoint(region.getSomePoint())){
                         region.setCheckResult(RegionCheckResult::ALLSAT);
                     }
                     else{
@@ -308,30 +308,26 @@ namespace storm {
                          proveAllSat=true;
                 }
 
-                bool formulaSatisfied = this->valueIsInBoundOfFormula(this->getApproximativeReachabilityValue(region, proveAllSat));
-
-                //check if approximation was conclusive
-                if(proveAllSat && formulaSatisfied){
-                    region.setCheckResult(RegionCheckResult::ALLSAT);
-                    return true;
-                }
-                if(!proveAllSat && !formulaSatisfied){
-                    region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                if(this->checkRegionWithApproximation(region, proveAllSat)){
+                    //approximation was conclusive
+                    if(proveAllSat){
+                        region.setCheckResult(RegionCheckResult::ALLSAT);
+                    } else {
+                        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;
-                    formulaSatisfied = this->valueIsInBoundOfFormula(this->getApproximativeReachabilityValue(region, proveAllSat));
-
-                    //check if approximation was conclusive
-                    if(proveAllSat && formulaSatisfied){
-                        region.setCheckResult(RegionCheckResult::ALLSAT);
-                        return true;
-                    }
-                    if(!proveAllSat && !formulaSatisfied){
-                        region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                    if(this->checkRegionWithApproximation(region, proveAllSat)){
+                        //approximation was conclusive
+                        if(proveAllSat){
+                            region.setCheckResult(RegionCheckResult::ALLSAT);
+                        } else {
+                            region.setCheckResult(RegionCheckResult::ALLVIOLATED);
+                        }
                         return true;
                     }
                 }
@@ -340,11 +336,15 @@ namespace storm {
             }
             
             template<typename ParametricSparseModelType, typename ConstantType>
-            ConstantType AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getApproximativeReachabilityValue(ParameterRegion<ParametricType> const& region, bool proveAllSat){
+            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::checkRegionWithApproximation(ParameterRegion<ParametricType> const& region, bool proveAllSat){
+                if(this->isResultConstant()){
+                    return (proveAllSat==this->checkFormulaOnSamplingPoint(region.getSomePoint()));
+                }
                 bool computeLowerBounds = (this->specifiedFormulaHasLowerBound() && proveAllSat) || (!this->specifiedFormulaHasLowerBound() && !proveAllSat);
-                return this->getApproximationModel()->computeInitialStateValue(region, computeLowerBounds);
+                bool formulaSatisfied = this->getApproximationModel()->checkFormulaOnRegion(region, computeLowerBounds);
+                return (proveAllSat==formulaSatisfied);
             }
-            
+                
             template<typename ParametricSparseModelType, typename ConstantType>
             std::shared_ptr<ApproximationModel<ParametricSparseModelType, ConstantType>> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getApproximationModel() {
                 if(this->approximationModel==nullptr){
@@ -373,6 +373,14 @@ namespace storm {
                 return this->getSamplingModel()->computeInitialStateValue(point);
             }
             
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::checkFormulaOnSamplingPoint(std::map<VariableType, CoefficientType> const& point) {
+                if(this->isResultConstant()){
+                    return this->valueIsInBoundOfFormula(this->constantResult.get());
+                }
+                return this->getSamplingModel()->checkFormulaOnSamplingPoint(point);
+            }
+            
             template<typename ParametricSparseModelType, typename ConstantType>
             std::shared_ptr<SamplingModel<ParametricSparseModelType, ConstantType>> const& AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::getSamplingModel() {
                 if(this->samplingModel==nullptr){
diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.h b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
index d4cfedc78..a4b3f6fb5 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.h
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
@@ -82,13 +82,22 @@ namespace storm {
                 ConstantType getReachabilityValue(std::map<VariableType, CoefficientType>const& point);
                 
                 /*!
-                 * Returns the approximative Value for the given region by instantiating and checking the approximation model. 
+                 * Computes the reachability Value at the specified point by instantiating and checking the sampling model. 
+                 * @param point The point (i.e. parameter evaluation) at which to compute the reachability value.
+                 * @return true iff the specified formula is satisfied
+                 */
+                bool checkFormulaOnSamplingPoint(std::map<VariableType, CoefficientType>const& point);
+                
+                /*!
+                 * Computes the approximative Value for the given region by instantiating and checking the approximation model. 
+                 * returns true iff the provided formula is satisfied w.r.t. the approximative value
                  * 
                  * @param region The region for which to compute the approximative value
-                 * @param proveAllSat if set to true, the returned value can be used to prove that the property is SATISFIED for all parameters in the given region. (false for VIOLATED)
-                 * @return a lower or upper bound of the actual reachability value (depending on the given flag 'proveAllSat' and whether the specified formula has a lower or an upper bound)
+                 * @param proveAllSat if set to true, it is checked whether the property is satisfied for all parameters in the given region. Otherwise, it is checked
+                          whether the property is violated for all parameters.
+                 * @return true iff the objective (given by the proveAllSat flag) was accomplished.
                  */
-                ConstantType getApproximativeReachabilityValue(ParameterRegion<ParametricType> const& region, bool proveAllSat);
+                bool checkRegionWithApproximation(ParameterRegion<ParametricType> const& region, bool proveAllSat);
                 
                 /*!
                  * Returns true iff the given value satisfies the bound given by the specified property
diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp
index 887ad08e2..5a0ec1316 100644
--- a/src/modelchecker/region/ApproximationModel.cpp
+++ b/src/modelchecker/region/ApproximationModel.cpp
@@ -276,8 +276,8 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             std::vector<ConstantType>  ApproximationModel<ParametricSparseModelType, ConstantType>::computeValues(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
                 instantiate(region, computeLowerBounds);
-                std::vector<std::size_t> policy;
-                invokeSolver(computeLowerBounds, policy);
+                Policy& policy = computeLowerBounds ? this->solverData.lastMinimizingPolicy : this->solverData.lastMaximizingPolicy;
+                invokeSolver(computeLowerBounds, policy, false);
                 
                 std::vector<ConstantType> result(this->maybeStates.size());
                 storm::utility::vector::setVectorValues(result, this->maybeStates, this->solverData.result);
@@ -291,63 +291,18 @@ namespace storm {
             ConstantType  ApproximationModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
                 instantiate(region, computeLowerBounds);
                 Policy& policy = computeLowerBounds ? this->solverData.lastMinimizingPolicy : this->solverData.lastMaximizingPolicy;
-                //TODO: at this point, set policy to the one stored in the region.
-                invokeSolver(computeLowerBounds, policy);
-    /*
-                //TODO: (maybe) when a few parameters are mapped to another value, build a "nicer" scheduler and check whether it induces values that are more optimal.
-                //Get the set of parameters which are (according to the policy) always mapped to the same region boundary.
-                //First, collect all (relevant) parameters, i.e., the ones that are set to the lower or upper boundary.
-                std::map<VariableType, RegionBoundary> fixedVariables;
-                std::map<VariableType, std::pair<std::size_t, std::size_t>> VarCount;
-                std::size_t substitutionCount =0;
-                for(auto const& substitution : this->funcSubData.substitutions){
-                    for( auto const& sub : substitution){
-                        if(sub.second!= RegionBoundary::UNSPECIFIED){
-                            fixedVariables.insert(typename std::map<VariableType, RegionBoundary>::value_type(sub.first, RegionBoundary::UNSPECIFIED));
-                            VarCount.insert(typename std::map<VariableType, std::pair<std::size_t, std::size_t>>::value_type(sub.first, std::pair<std::size_t, std::size_t>(0,0)));
-                        }
-                    }
-                }
-                //Now erase the parameters that are mapped to different boundaries.
-                for(std::size_t rowGroup=0; rowGroup<this->matrixData.matrix.getRowGroupCount(); ++rowGroup){
-                    std::size_t row = this->matrixData.matrix.getRowGroupIndices()[rowGroup] + policy[rowGroup];
-                    for(std::pair<VariableType, RegionBoundary> const& sub : this->funcSubData.substitutions[this->matrixData.rowSubstitutions[row]]){
-                        auto fixedVarIt = fixedVariables.find(sub.first);
-                        if(fixedVarIt != fixedVariables.end() && fixedVarIt->second != sub.second){
-                            if(fixedVarIt->second == RegionBoundary::UNSPECIFIED){
-                                fixedVarIt->second = sub.second;
-                            } else {
-                                // the solution maps  the current variable at least once to lower boundary and at least once to the upper boundary.
-                                fixedVariables.erase(fixedVarIt);
-                            }
-                        }
-                        auto varcountIt = VarCount.find(sub.first);
-                        if(sub.second==RegionBoundary::LOWER){
-                            ++(varcountIt->second.first);
-                        } else if (sub.second==RegionBoundary::UPPER){
-                            ++(varcountIt->second.second);
-                        }
-                        ++substitutionCount;
-                    }
-                    if (fixedVariables.empty()){
-                       // break;
-                    }
-                }
-       //         std::cout << "Used Approximation" << std::endl;
-                for (auto const& varcount : VarCount){
-                    if(varcount.second.first > 0 && varcount.second.second > 0){
-                        std::cout << "  Variable " << varcount.first << " has been set to lower " << varcount.second.first << " times and to upper " << varcount.second.second << " times. (total: " << substitutionCount << ", " << (computeLowerBounds ? "MIN" : "MAX") << ")" << std::endl;
-                    }
-                }
-                for (auto const& fixVar : fixedVariables){
-                    //std::cout << "  APPROXMODEL: variable " << fixVar.first << " is always mapped to " << fixVar.second << std::endl;
-                }
-                    
-        //        std::cout << "    Result is " << this->solverData.result[this->solverData.initialStateIndex] << std::endl;
-     */
+                invokeSolver(computeLowerBounds, policy, false);
                 return this->solverData.result[this->solverData.initialStateIndex];
             }
             
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool  ApproximationModel<ParametricSparseModelType, ConstantType>::checkFormulaOnRegion(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
+                instantiate(region, computeLowerBounds);
+                Policy& policy = computeLowerBounds ? this->solverData.lastMinimizingPolicy : this->solverData.lastMaximizingPolicy;
+                invokeSolver(computeLowerBounds, policy, true); //allow early termination
+                return this->solverData.player1Goal->achieved(this->solverData.result);
+            }
+            
             template<typename ParametricSparseModelType, typename ConstantType>
             void ApproximationModel<ParametricSparseModelType, ConstantType>::instantiate(const ParameterRegion<ParametricType>& region, bool computeLowerBounds) {
                 //Instantiate the substitutions
diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h
index 21322f026..b3c61a18b 100644
--- a/src/modelchecker/region/ApproximationModel.h
+++ b/src/modelchecker/region/ApproximationModel.h
@@ -56,6 +56,14 @@ namespace storm {
                  */
                 ConstantType computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
 
+                /*!
+                 * Instantiates the approximation model w.r.t. the given region.
+                 * Then computes the approximated reachability probabilities or reward value for the initial state.
+                 * If computeLowerBounds is true, the computed value will be a lower bound for the actual value. Otherwise, we get an upper bound.
+                 * Returns true iff the computed bound satisfies the formula (given upon construction of *this)
+                 */
+                bool checkFormulaOnRegion(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
+
             private:
 
                 typedef std::pair<ParametricType, std::size_t> FunctionSubstitution;
@@ -76,7 +84,7 @@ namespace storm {
                 void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
                 void initializePlayer1Matrix(ParametricSparseModelType const& parametricModel);
                 void instantiate(ParameterRegion<ParametricType> const& region, bool computeLowerBounds);
-                void invokeSolver(bool computeLowerBounds, Policy& policy, bool allowEarlyTermination=true);
+                void invokeSolver(bool computeLowerBounds, Policy& policy, bool allowEarlyTermination);
 
                 //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 dc361129e..58baed2b9 100644
--- a/src/modelchecker/region/SamplingModel.cpp
+++ b/src/modelchecker/region/SamplingModel.cpp
@@ -223,7 +223,7 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             std::vector<ConstantType> SamplingModel<ParametricSparseModelType, ConstantType>::computeValues(std::map<VariableType, CoefficientType>const& point) {
                 instantiate(point);
-                invokeSolver();
+                invokeSolver(false); //no early termination
                 std::vector<ConstantType> result(this->maybeStates.size());
                 storm::utility::vector::setVectorValues(result, this->maybeStates, this->solverData.result);
                 storm::utility::vector::setVectorValues(result, this->targetStates, this->computeRewards ? storm::utility::zero<ConstantType>() : storm::utility::one<ConstantType>());
@@ -235,10 +235,17 @@ namespace storm {
             template<typename ParametricSparseModelType, typename ConstantType>
             ConstantType SamplingModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue(std::map<VariableType, CoefficientType>const& point) {
                 instantiate(point);
-                invokeSolver();
+                invokeSolver(false); //no early termination
                 return this->solverData.result[this->solverData.initialStateIndex];
             }
             
+            template<typename ParametricSparseModelType, typename ConstantType>
+            bool SamplingModel<ParametricSparseModelType, ConstantType>::checkFormulaOnSamplingPoint(std::map<VariableType, CoefficientType>const& point) {
+                instantiate(point);
+                invokeSolver(true); //allow early termination
+                return this->solverData.solveGoal->achieved(this->solverData.result);
+            }
+            
             template<typename ParametricSparseModelType, typename ConstantType>
             void SamplingModel<ParametricSparseModelType, ConstantType>::instantiate(std::map<VariableType, CoefficientType>const& point) {
                 //write results into the placeholders
diff --git a/src/modelchecker/region/SamplingModel.h b/src/modelchecker/region/SamplingModel.h
index e5c65ff49..83fdb6b0e 100644
--- a/src/modelchecker/region/SamplingModel.h
+++ b/src/modelchecker/region/SamplingModel.h
@@ -46,9 +46,14 @@ namespace storm {
                 /*!
                  * Instantiates the underlying model according to the given point
                  * Returns the reachability probability (or the expected rewards) of the initial state.
-                 * Undefined behavior if model has not been instantiated first!
                  */
                 ConstantType computeInitialStateValue(std::map<VariableType, CoefficientType>const& point);
+                
+                /*!
+                 * Instantiates the underlying model according to the given point
+                 * Returns true iff the formula (given upon construction of *this) is true in the initial state of the instantiated model
+                 */
+                bool checkFormulaOnSamplingPoint(std::map<VariableType, CoefficientType>const& point);
 
             private:
                 typedef typename std::unordered_map<ParametricType, ConstantType>::value_type FunctionEntry;
@@ -57,7 +62,7 @@ namespace storm {
                 void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
                 void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices);
                 void instantiate(std::map<VariableType, CoefficientType>const& point);
-                void invokeSolver(bool allowEarlyTermination=true);
+                void invokeSolver(bool allowEarlyTermination);
                 
                 //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 fee98db91..cfe45c005 100644
--- a/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
+++ b/src/modelchecker/region/SparseDtmcRegionModelChecker.cpp
@@ -429,7 +429,7 @@ namespace storm {
                 }
                 else{
                     //instantiate the sampling model
-                    valueInBoundOfFormula = this->valueIsInBoundOfFormula(this->getReachabilityValue(point));
+                    valueInBoundOfFormula = this->checkFormulaOnSamplingPoint(point);
                 }
 
                 if(valueInBoundOfFormula){
diff --git a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp
index 542c85e4c..3653c8e97 100644
--- a/src/modelchecker/region/SparseMdpRegionModelChecker.cpp
+++ b/src/modelchecker/region/SparseMdpRegionModelChecker.cpp
@@ -247,7 +247,7 @@ namespace storm {
 
             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(this->checkFormulaOnSamplingPoint(point)){
                     if (region.getCheckResult()!=RegionCheckResult::EXISTSSAT){
                         region.setSatPoint(point);
                         if(region.getCheckResult()==RegionCheckResult::EXISTSVIOLATED){
diff --git a/src/solver/SolveGoal.h b/src/solver/SolveGoal.h
index fe573dd3a..3da781e75 100644
--- a/src/solver/SolveGoal.h
+++ b/src/solver/SolveGoal.h
@@ -54,6 +54,26 @@ namespace storm {
             VT thresholdValue() const { return threshold; }
             storm::storage::BitVector relevantColumns() const { return relevantColumnVector; }
             
+            bool achieved(std::vector<VT> const& result) const{
+                for(std::size_t i : relevantColumnVector){
+                    switch(boundType) {
+                    case storm::logic::ComparisonType::Greater:
+                        if( result[i] <= threshold) return false;
+                        break;
+                    case storm::logic::ComparisonType::GreaterEqual:
+                        if( result[i] < threshold) return false;
+                        break;
+                    case storm::logic::ComparisonType::Less:
+                        if( result[i] >= threshold) return false;
+                        break;
+                    case storm::logic::ComparisonType::LessEqual:
+                        if( result[i] > threshold) return false;
+                        break;
+                    }
+                }
+                return true;
+            }
+            
             storm::logic::ComparisonType boundType;
             VT threshold;
             storm::storage::BitVector relevantColumnVector;
diff --git a/src/utility/policyguessing.cpp b/src/utility/policyguessing.cpp
index 1981bc42d..018d83085 100644
--- a/src/utility/policyguessing.cpp
+++ b/src/utility/policyguessing.cpp
@@ -186,9 +186,9 @@ namespace storm {
                     storm::utility::vector::addVectors(copyX, subB, copyX); // = Ax + b
                     ++iterations;
                     newPrecision *= 0.5;
-                } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision*0.5, relative) && iterations<50);
+                } while(!storm::utility::vector::equalModuloPrecision(subX, copyX, precision*0.5, relative) && iterations<60);
                 
-                STORM_LOG_WARN_COND(iterations<50, "Solving linear equation system did not yield a precise result");
+                STORM_LOG_WARN_COND(iterations<60, "Solving linear equation system did not yield a precise result");
                 
                 STORM_LOG_DEBUG("Required to increase the precision " << iterations << " times in order to obtain a precise result");
                 //fill in the result
diff --git a/src/utility/storm.h b/src/utility/storm.h
index 643c826a0..91f12e1fb 100644
--- a/src/utility/storm.h
+++ b/src/utility/storm.h
@@ -368,12 +368,12 @@ namespace storm {
      * @param upperBoundaries maps every variable to its highest possible value within the region. (corresponds to the top right corner point in the 2D case)
      * @param proveAllSat if set to true, it is checked whether the property is satisfied for all parameters in the given region. Otherwise, it is checked
      *                    whether the property is violated for all parameters.
-     * @return true iff the objective was accomplished.
+     * @return true iff the objective (given by the proveAllSat flag) was accomplished.
      * 
      * So there are the following cases:
      * proveAllSat=true,  return=true  ==> the property is SATISFIED for all parameters in the given region
      * proveAllSat=true,  return=false ==> the approximative value is NOT within the bound of the formula (either the approximation is too bad or there are points in the region that violate the property)
-     * proveAllSat=false, return=true  ==> the property is VIOLATED for all parameters in teh given region
+     * proveAllSat=false, return=true  ==> the property is VIOLATED for all parameters in the given region
      * proveAllSat=false, return=false ==> the approximative value IS within the bound of the formula (either the approximation is too bad or there are points in the region that satisfy the property)
      */
     inline bool checkRegionApproximation(std::shared_ptr<storm::modelchecker::region::AbstractSparseRegionModelChecker<storm::models::sparse::Model<storm::RationalFunction>, double>> regionModelChecker,
@@ -381,14 +381,7 @@ namespace storm {
                                          std::map<storm::Variable, storm::RationalNumber> const& upperBoundaries,
                                          bool proveAllSat){
         storm::modelchecker::region::ParameterRegion<storm::RationalFunction> region(lowerBoundaries, upperBoundaries);
-        double result=regionModelChecker->getApproximativeReachabilityValue(region, proveAllSat);
-        if(proveAllSat && regionModelChecker->valueIsInBoundOfFormula(result)){
-            return true;
-        } else if (!proveAllSat && !regionModelChecker->valueIsInBoundOfFormula(result)){
-            return true;
-        } else {
-            return false;
-        }
+        return regionModelChecker->checkRegionWithApproximation(region, proveAllSat);
     }
     
             
diff --git a/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp b/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
index 9f0b774ca..59551522d 100644
--- a/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseDtmcRegionModelCheckerTest.cpp
@@ -30,6 +30,9 @@ 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.73,0.2<=pK<=0.715");
 
+    EXPECT_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(dtmcModelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
+    
     EXPECT_NEAR(0.8369631407, dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.8369631407, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getLowerBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.0476784174, dtmcModelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
@@ -89,6 +92,9 @@ 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_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(dtmcModelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
+    
     EXPECT_NEAR(4.367791292, dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(4.367791292, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getLowerBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(3.044795147, dtmcModelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
@@ -167,7 +173,8 @@ TEST(SparseDtmcRegionModelCheckerTest, Brp_Rew_Infty) {
     
     //start testing
     auto allSatRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("");
-
+    
+    EXPECT_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
     EXPECT_EQ(storm::utility::infinity<double>(), dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()));
     
     //test approximative method
@@ -206,6 +213,9 @@ 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_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(dtmcModelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
+    
     EXPECT_NEAR(4.834779705, dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(4.834779705, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getLowerBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(4.674651623, dtmcModelchecker->getReachabilityValue(exBothRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
@@ -260,6 +270,9 @@ 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_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(dtmcModelchecker->checkFormulaOnSamplingPoint(allVioHardRegion.getSomePoint()));
+    
     EXPECT_NEAR(0.1734086422, dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.1734086422, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getLowerBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.47178, dtmcModelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
@@ -337,6 +350,9 @@ 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_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(dtmcModelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
+    
     EXPECT_NEAR(0.8430128158, dtmcModelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.8430128158, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getUpperBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.7731321947, dtmcModelchecker->getReachabilityValue(exBothRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
@@ -390,6 +406,8 @@ TEST(SparseDtmcRegionModelCheckerTest, Crowds_Prob_Const) {
     //start testing
     auto allSatRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("");
     
+    EXPECT_TRUE(dtmcModelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    
     EXPECT_NEAR(0.6119660237, dtmcModelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.6119660237, storm::utility::region::convertNumber<double>(dtmcModelchecker->evaluateReachabilityFunction(allSatRegion.getUpperBoundaries())),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.6119660237, dtmcModelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()),  storm::settings::generalSettings().getPrecision());
diff --git a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp
index 6a6db1078..bedb3048b 100644
--- a/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseMdpRegionModelCheckerTest.cpp
@@ -27,6 +27,9 @@ TEST(SparseMdpRegionModelCheckerTest, two_dice_Prob) {
     auto exBothRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.45<=p1<=0.55,0.45<=p2<=0.55");
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.6<=p1<=0.7,0.6<=p2<=0.6");
 
+    EXPECT_TRUE(modelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(modelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
+    
     //Test the methods provided in storm.h
     EXPECT_TRUE(storm::checkSamplingPoint(modelchecker,allSatRegion.getLowerBoundaries()));
     EXPECT_TRUE(storm::checkSamplingPoint(modelchecker,allSatRegion.getUpperBoundaries()));
@@ -82,8 +85,11 @@ TEST(SparseMdpRegionModelCheckerTest, coin_Prob) {
     auto allSatRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.3<=p1<=0.45,0.2<=p2<=0.54");
     auto exBothRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.4<=p1<=0.65,0.5<=p2<=0.7");
     auto allVioRegion=storm::modelchecker::region::ParameterRegion<storm::RationalFunction>::parseRegion("0.4<=p1<=0.7,0.55<=p2<=0.6");
+    
+    EXPECT_TRUE(modelchecker->checkFormulaOnSamplingPoint(allSatRegion.getSomePoint()));
+    EXPECT_FALSE(modelchecker->checkFormulaOnSamplingPoint(allVioRegion.getSomePoint()));
 
-    EXPECT_NEAR(0.95127874851, modelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
+    EXPECT_NEAR(0.95128124239, modelchecker->getReachabilityValue(allSatRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.26787251126, modelchecker->getReachabilityValue(allSatRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.41880006098, modelchecker->getReachabilityValue(exBothRegion.getLowerBoundaries()), storm::settings::generalSettings().getPrecision());
     EXPECT_NEAR(0.01535089684, modelchecker->getReachabilityValue(exBothRegion.getUpperBoundaries()),  storm::settings::generalSettings().getPrecision());