diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
index f2afacbde..c8dc72c54 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.cpp
@@ -218,11 +218,9 @@ namespace storm {
                 }
 
                 std::chrono::high_resolution_clock::time_point timeApproximationStart = std::chrono::high_resolution_clock::now();
-                std::vector<ConstantType> lowerBounds;
-                std::vector<ConstantType> upperBounds;
                 if(!done && doApproximation){
                     STORM_LOG_DEBUG("Checking approximative values...");
-                    if(this->checkApproximativeValues(region, lowerBounds, upperBounds)){
+                    if(this->checkApproximativeValues(region)){
                         ++this->numOfRegionsSolvedThroughApproximation;
                         STORM_LOG_DEBUG("Result '" << region.getCheckResult() <<"' obtained through approximation.");
                         done=true;
@@ -274,11 +272,8 @@ namespace storm {
             }
             
             template<typename ParametricSparseModelType, typename ConstantType>
-            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::checkApproximativeValues(ParameterRegion<ParametricType>& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds) {
+            bool AbstractSparseRegionModelChecker<ParametricSparseModelType, ConstantType>::checkApproximativeValues(ParameterRegion<ParametricType>& region) {
                 std::chrono::high_resolution_clock::time_point timeMDPBuildStart = std::chrono::high_resolution_clock::now();
-                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;
@@ -317,19 +312,9 @@ namespace storm {
                          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)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
-                    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)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
-                    upperBounds = std::vector<ConstantType>();
-                    formulaSatisfied = this->valueIsInBoundOfFormula(lowerBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
-                }
+                bool computeLowerBounds = (!this->specifiedFormulaHasUpperBound() && proveAllSat) || (this->specifiedFormulaHasUpperBound() && !proveAllSat);
+                STORM_LOG_DEBUG("Approximation for " << (computeLowerBounds ? "lower" : "upper") << "bounds");
+                bool formulaSatisfied = this->valueIsInBoundOfFormula(this->getApproximationModel()->computeInitialStateValue(region, computeLowerBounds));
 
                 //check if approximation was conclusive
                 if(proveAllSat && formulaSatisfied){
@@ -344,15 +329,8 @@ namespace storm {
                 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)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
-                        formulaSatisfied=this->valueIsInBoundOfFormula(lowerBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
-                    }
-                    else{
-                        upperBounds = this->getApproximationModel()->computeValues(storm::solver::OptimizationDirection::Maximize)->template asExplicitQuantitativeCheckResult<ConstantType>().getValueVector();
-                        formulaSatisfied=this->valueIsInBoundOfFormula(upperBounds[*this->getApproximationModel()->getModel()->getInitialStates().begin()]);
-                    }
+                    computeLowerBounds=!computeLowerBounds;
+                    formulaSatisfied = this->valueIsInBoundOfFormula(this->getApproximationModel()->computeInitialStateValue(region, computeLowerBounds));
 
                     //check if approximation was conclusive
                     if(proveAllSat && formulaSatisfied){
diff --git a/src/modelchecker/region/AbstractSparseRegionModelChecker.h b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
index 415145e69..0f1fc8dd2 100644
--- a/src/modelchecker/region/AbstractSparseRegionModelChecker.h
+++ b/src/modelchecker/region/AbstractSparseRegionModelChecker.h
@@ -130,7 +130,7 @@ namespace storm {
                  * 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.
                  */
-                bool checkApproximativeValues(ParameterRegion<ParametricType>& region, std::vector<ConstantType>& lowerBounds, std::vector<ConstantType>& upperBounds); 
+                bool checkApproximativeValues(ParameterRegion<ParametricType>& region); 
 
                 /*!
                  * Returns the approximation model.
diff --git a/src/modelchecker/region/ApproximationModel.cpp b/src/modelchecker/region/ApproximationModel.cpp
index 2f72d61e7..42b2a1181 100644
--- a/src/modelchecker/region/ApproximationModel.cpp
+++ b/src/modelchecker/region/ApproximationModel.cpp
@@ -35,7 +35,7 @@ namespace storm {
                     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");
                 } else {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << this->formula << ". Approximation model only supports eventually or reachability reward formulae.");
+                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Invalid formula: " << formula << ". Approximation model only supports eventually or reachability reward formulae.");
                 }
                 STORM_LOG_THROW(parametricModel.hasLabel("target"), storm::exceptions::InvalidArgumentException, "The given Model has no \"target\"-statelabel.");
                 this->targetStates = parametricModel.getStateLabeling().getStates("target");
@@ -63,7 +63,7 @@ namespace storm {
                 this->matrixData.assignment.shrink_to_fit();
                 this->vectorData.assignment.shrink_to_fit();
                 if(parametricModel.getType()==storm::models::ModelType::Mdp){
-                    initializePlayer1Matrix(parametricModel, newIndices);
+                    initializePlayer1Matrix(parametricModel);
                 }
                 
                 this->eqSysResult = std::vector<ConstantType>(maybeStates.getNumberOfSetBits(), this->computeRewards ? storm::utility::one<ConstantType>() : ConstantType(0.5));
@@ -86,7 +86,7 @@ namespace storm {
                                                                                 false, // no force dimensions
                                                                                 true, //will have custom row grouping
                                                                                 numOfMaybeStates); //exact number of rowgroups is unknown at this point, but at least this many
-                rowSubstitutions.push_back(numOfMaybeStates);
+                rowSubstitutions.reserve(numOfMaybeStates);
                 std::size_t curRow = 0;
                 for (auto oldRowGroup : this->maybeStates){
                     for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){
@@ -123,8 +123,8 @@ namespace storm {
                                     matrixBuilder.addNextValue(curRow, newIndices[oldEntry.getColumn()], dummyValue);
                                 }
                             }
+                            ++curRow;
                         }
-                        ++curRow;
                     }
                 }
                 this->matrixData.matrix=matrixBuilder.build();               
@@ -198,7 +198,7 @@ namespace storm {
                     if(storm::utility::isConstant(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState])){
                         ConstantType reward = storm::utility::region::convertNumber<ConstantType>(storm::utility::region::getConstantPart(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState]));
                         //Add one of these entries for every row in the row group of oldState
-                        for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[state+1]; ++matrixRow){
+                        for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[oldState+1]; ++matrixRow){
                             *vectorIt = reward;
                             ++vectorIt;
                         }
@@ -209,7 +209,7 @@ namespace storm {
                         // We might find out that the reward is independent of the probability parameters (and will thus be independent of nondeterministic choices)
                         // In that case, the reward function and the substitution will not change and thus we can use the same FunctionSubstitution
                         bool rewardDependsOnProbVars=true;
-                        std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::iterator functionsIt;
+                        typename std::unordered_map<FunctionSubstitution, ConstantType, FuncSubHash>::iterator functionsIt;
                         for(auto matrixRow=this->matrixData.matrix.getRowGroupIndices()[oldState]; matrixRow<this->matrixData.matrix.getRowGroupIndices()[oldState+1]; ++matrixRow){
                             auto probabilitySub = this->funcSubData.substitutions[rowSubstitutions[matrixRow]];
                             if(rewardDependsOnProbVars){ //always executed in first iteration
@@ -227,7 +227,7 @@ namespace storm {
                                 }
                                 // insert the substitution and the FunctionSubstitution
                                 std::size_t substitutionIndex=storm::utility::vector::findOrInsert(this->funcSubData.substitutions, std::move(substitution));
-                                functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[state], substitutionIndex), dummyValue)).first;
+                                functionsIt = this->funcSubData.functions.insert(FunctionEntry(FunctionSubstitution(parametricModel.getUniqueRewardModel()->second.getStateRewardVector()[oldState], substitutionIndex), dummyValue)).first;
                             }
                             //insert assignment and dummy data
                             this->vectorData.assignment.emplace_back(std::make_pair(vectorIt, &(functionsIt->second)));
@@ -239,7 +239,25 @@ namespace storm {
                 STORM_LOG_THROW(vectorIt==this->vectorData.vector.end(), storm::exceptions::UnexpectedException, "initRewards: The size of the eq-sys vector is not as it was expected");
             }
 
-
+            template<typename ParametricSparseModelType, typename ConstantType>
+            void ApproximationModel<ParametricSparseModelType, ConstantType>::initializePlayer1Matrix(ParametricSparseModelType const& parametricModel){
+                std::size_t p1MatrixSize = matrixData.matrix.getRowGroupCount();
+                storm::storage::SparseMatrixBuilder<storm::storage::sparse::state_type> matrixBuilder(p1MatrixSize, //rows
+                                                                                p1MatrixSize, //columns
+                                                                                p1MatrixSize, //entries
+                                                                                true, // force dimensions
+                                                                                true, //will have custom row grouping
+                                                                                this->maybeStates.getNumberOfSetBits()); // number of rowgroups
+                std::size_t curRow = 0;
+                for (auto oldRowGroup : this->maybeStates){
+                    matrixBuilder.newRowGroup(curRow);
+                    for (std::size_t oldRow = parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup]; oldRow < parametricModel.getTransitionMatrix().getRowGroupIndices()[oldRowGroup+1]; ++oldRow){
+                        matrixBuilder.addNextValue(curRow,curRow, storm::utility::one<storm::storage::sparse::state_type>());
+                        ++curRow;
+                    }
+                }
+                this->player1Matrix = matrixBuilder.build();
+            }
 
             template<typename ParametricSparseModelType, typename ConstantType>
             ApproximationModel<ParametricSparseModelType, ConstantType>::~ApproximationModel() {
@@ -263,6 +281,7 @@ namespace storm {
             ConstantType  ApproximationModel<ParametricSparseModelType, ConstantType>::computeInitialStateValue(ParameterRegion<ParametricType> const& region, bool computeLowerBounds) {
                 instantiate(region, computeLowerBounds);
                 invokeSolver(computeLowerBounds);
+  //              std::cout << "initialStateValue is " << this->eqSysResult[this->eqSysInitIndex] << std::endl;
                 return this->eqSysResult[this->eqSysInitIndex];
             }
             
@@ -297,7 +316,7 @@ namespace storm {
                     auto& result = functionResult.second;
                     result = computeLowerBounds ? storm::utility::infinity<ConstantType>() : storm::utility::zero<ConstantType>();
                      //Iterate over the different combinations of lower and upper bounds and update the min and max values
-                    auto const& vertices=region.getVerticesOfRegion(this->choseOptimalParameters[funcSub.getSubstitution()]);
+                    auto const& vertices=region.getVerticesOfRegion(choseOptimalParameters[funcSub.getSubstitution()]);
                     for(auto const& vertex : vertices){
                         //extend the substitution
                         for(auto const& vertexSub : vertex){
diff --git a/src/modelchecker/region/ApproximationModel.h b/src/modelchecker/region/ApproximationModel.h
index 4b11f22f2..ca1938c8f 100644
--- a/src/modelchecker/region/ApproximationModel.h
+++ b/src/modelchecker/region/ApproximationModel.h
@@ -127,7 +127,7 @@ namespace storm {
 
                 void initializeProbabilities(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t>& rowSubstitutions);
                 void initializeRewards(ParametricSparseModelType const& parametricModel, std::vector<std::size_t> const& newIndices, std::vector<std::size_t> const& rowSubstitutions);
-                void initializePlayer1Matrix(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);
 
@@ -142,7 +142,7 @@ namespace storm {
                 bool computeRewards;
                 //Player 1 represents the nondeterminism of the given mdp (so, this is irrelevant if we approximate values of a DTMC)
                 storm::solver::SolveGoal player1Goal;
-                storm::storage::SparseMatrix<storm::storage::sparse::state_type> const& player1Matrix;
+                storm::storage::SparseMatrix<storm::storage::sparse::state_type> player1Matrix;
                 
                 /* The data required for the equation system, i.e., a matrix and a vector.
                  * 
diff --git a/src/utility/region.cpp b/src/utility/region.cpp
index 2592f3a68..128762b2a 100644
--- a/src/utility/region.cpp
+++ b/src/utility/region.cpp
@@ -17,6 +17,7 @@
 
 #ifdef STORM_HAVE_CARL
 #include<carl/numbers/numbers.h>
+#include<carl/core/VariablePool.h>
 #endif
 
 namespace storm {