diff --git a/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp b/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp
index dbbd9d7c2..6c2f15d21 100644
--- a/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp
+++ b/src/storm-pars/modelchecker/region/SparseDtmcParameterLiftingModelChecker.cpp
@@ -10,7 +10,8 @@
 #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h"
 #include "storm/models/sparse/Dtmc.h"
 #include "storm/models/sparse/StandardRewardModel.h"
-#include "storm/solver/StandardMinMaxLinearEquationSolver.h"
+#include "storm/solver/MinMaxLinearEquationSolver.h"
+#include "storm/solver/Multiplier.h"
 #include "storm/utility/vector.h"
 #include "storm/utility/graph.h"
 #include "storm/utility/NumberTraits.h"
@@ -249,49 +250,50 @@ namespace storm {
             
             parameterLifter->specifyRegion(region, dirForParameters);
             
-            auto solver = solverFactory->create(env, parameterLifter->getMatrix());
-            solver->setHasUniqueSolution();
-            if (lowerResultBound) solver->setLowerBound(lowerResultBound.get());
-            if (upperResultBound) {
-                solver->setUpperBound(upperResultBound.get());
-            } else if (solvingRequiresUpperRewardBounds) {
-                // For the min-case, we use DS-MPI, for the max-case variant 2 of the Baier et al. paper (CAV'17).
-                std::vector<ConstantType> oneStepProbs;
-                oneStepProbs.reserve(parameterLifter->getMatrix().getRowCount());
-                for (uint64_t row = 0; row < parameterLifter->getMatrix().getRowCount(); ++row) {
-                    oneStepProbs.push_back(storm::utility::one<ConstantType>() - parameterLifter->getMatrix().getRowSum(row));
-                }
-                if (dirForParameters == storm::OptimizationDirection::Minimize) {
-                    storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer<ConstantType> dsmpi(parameterLifter->getMatrix(), parameterLifter->getVector(), oneStepProbs);
-                    solver->setUpperBounds(dsmpi.computeUpperBounds());
-                } else {
-                    storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ConstantType> baier(parameterLifter->getMatrix(), parameterLifter->getVector(), oneStepProbs);
-                    solver->setUpperBound(baier.computeUpperBound());
-                }
-            }
-            if (!stepBound) solver->setTrackScheduler(true);
-            if (storm::solver::minimize(dirForParameters) && minSchedChoices && !stepBound) solver->setInitialScheduler(std::move(minSchedChoices.get()));
-            if (storm::solver::maximize(dirForParameters) && maxSchedChoices && !stepBound) solver->setInitialScheduler(std::move(maxSchedChoices.get()));
-            if (this->currentCheckTask->isBoundSet() && solver->hasInitialScheduler()) {
-                // If we reach this point, we know that after applying the hint, the x-values can only become larger (if we maximize) or smaller (if we minimize).
-                std::unique_ptr<storm::solver::TerminationCondition<ConstantType>> termCond;
-                storm::storage::BitVector relevantStatesInSubsystem = this->currentCheckTask->isOnlyInitialStatesRelevantSet() ? this->parametricModel->getInitialStates() % maybeStates : storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true);
-                if (storm::solver::minimize(dirForParameters)) {
-                    // Terminate if the value for ALL relevant states is already below the threshold
-                    termCond = std::make_unique<storm::solver::TerminateIfFilteredExtremumBelowThreshold<ConstantType>> (relevantStatesInSubsystem, true, this->currentCheckTask->getBoundThreshold(), false);
-                } else {
-                    // Terminate if the value for ALL relevant states is already above the threshold
-                    termCond = std::make_unique<storm::solver::TerminateIfFilteredExtremumExceedsThreshold<ConstantType>> (relevantStatesInSubsystem, true, this->currentCheckTask->getBoundThreshold(), true);
-                }
-                solver->setTerminationCondition(std::move(termCond));
-            }
-            
-            // Invoke the solver
             if (stepBound) {
                 assert(*stepBound > 0);
                 x = std::vector<ConstantType>(maybeStates.getNumberOfSetBits(), storm::utility::zero<ConstantType>());
-                solver->repeatedMultiply(env, dirForParameters, x, &parameterLifter->getVector(), *stepBound);
+                auto multiplier = storm::solver::MultiplierFactory<ConstantType>().create();
+                multiplier->repeatedMultiply(env, dirForParameters, x, &parameterLifter->getVector(), *stepBound);
             } else {
+                auto solver = solverFactory->create(env, parameterLifter->getMatrix());
+                solver->setHasUniqueSolution();
+                if (lowerResultBound) solver->setLowerBound(lowerResultBound.get());
+                if (upperResultBound) {
+                    solver->setUpperBound(upperResultBound.get());
+                } else if (solvingRequiresUpperRewardBounds) {
+                    // For the min-case, we use DS-MPI, for the max-case variant 2 of the Baier et al. paper (CAV'17).
+                    std::vector<ConstantType> oneStepProbs;
+                    oneStepProbs.reserve(parameterLifter->getMatrix().getRowCount());
+                    for (uint64_t row = 0; row < parameterLifter->getMatrix().getRowCount(); ++row) {
+                        oneStepProbs.push_back(storm::utility::one<ConstantType>() - parameterLifter->getMatrix().getRowSum(row));
+                    }
+                    if (dirForParameters == storm::OptimizationDirection::Minimize) {
+                        storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer<ConstantType> dsmpi(parameterLifter->getMatrix(), parameterLifter->getVector(), oneStepProbs);
+                        solver->setUpperBounds(dsmpi.computeUpperBounds());
+                    } else {
+                        storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ConstantType> baier(parameterLifter->getMatrix(), parameterLifter->getVector(), oneStepProbs);
+                        solver->setUpperBound(baier.computeUpperBound());
+                    }
+                }
+                solver->setTrackScheduler(true);
+                if (storm::solver::minimize(dirForParameters) && minSchedChoices) solver->setInitialScheduler(std::move(minSchedChoices.get()));
+                if (storm::solver::maximize(dirForParameters) && maxSchedChoices) solver->setInitialScheduler(std::move(maxSchedChoices.get()));
+                if (this->currentCheckTask->isBoundSet() && solver->hasInitialScheduler()) {
+                    // If we reach this point, we know that after applying the hint, the x-values can only become larger (if we maximize) or smaller (if we minimize).
+                    std::unique_ptr<storm::solver::TerminationCondition<ConstantType>> termCond;
+                    storm::storage::BitVector relevantStatesInSubsystem = this->currentCheckTask->isOnlyInitialStatesRelevantSet() ? this->parametricModel->getInitialStates() % maybeStates : storm::storage::BitVector(maybeStates.getNumberOfSetBits(), true);
+                    if (storm::solver::minimize(dirForParameters)) {
+                        // Terminate if the value for ALL relevant states is already below the threshold
+                        termCond = std::make_unique<storm::solver::TerminateIfFilteredExtremumBelowThreshold<ConstantType>> (relevantStatesInSubsystem, true, this->currentCheckTask->getBoundThreshold(), false);
+                    } else {
+                        // Terminate if the value for ALL relevant states is already above the threshold
+                        termCond = std::make_unique<storm::solver::TerminateIfFilteredExtremumExceedsThreshold<ConstantType>> (relevantStatesInSubsystem, true, this->currentCheckTask->getBoundThreshold(), true);
+                    }
+                    solver->setTerminationCondition(std::move(termCond));
+                }
+            
+                // Invoke the solver
                 x.resize(maybeStates.getNumberOfSetBits(), storm::utility::zero<ConstantType>());
                 solver->solveEquations(env, dirForParameters, x, parameterLifter->getVector());
                 if(storm::solver::minimize(dirForParameters)) {