diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 14e133014..2a39daf3c 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -13,6 +13,9 @@ #include "storm/environment/Environment.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/environment/solver/TopologicalSolverEnvironment.h" +#include "storm/environment/solver/LongRunAverageSolverEnvironment.h" +#include "storm/environment/solver/EigenSolverEnvironment.h" #include "storm/utility/macros.h" #include "storm/utility/vector.h" @@ -703,7 +706,9 @@ namespace storm { auto underlyingSolverEnvironment = env; if (env.solver().isForceSoundness()) { // For sound computations, the error in the MECS plus the error in the remaining system should be less then the user defined precsion. - underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().minMax().getPrecision() / storm::utility::convertNumber(2)); + underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber(2)); + underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion()); + underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber(2)); } for (uint64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { @@ -877,11 +882,11 @@ namespace storm { } // Solve MEC with the method specified in the settings - storm::solver::LraMethod method = env.solver().minMax().getLraMethod(); - if (storm::NumberTraits::IsExact && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) { + storm::solver::LraMethod method = env.solver().lra().getNondetLraMethod(); + if (storm::NumberTraits::IsExact && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) { STORM_LOG_INFO("Selecting 'LP' as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method."); method = storm::solver::LraMethod::LinearProgramming; - } else if (env.solver().isForceSoundness() && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { + } else if (env.solver().isForceSoundness() && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { STORM_LOG_INFO("Selecting 'VI' as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method."); method = storm::solver::LraMethod::ValueIteration; } @@ -982,7 +987,7 @@ namespace storm { ValueType uniformizationRate = storm::utility::vector::max_if(exitRateVector, markovianMecStates); // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop. // Hence, we increase the uniformization rate a little. - uniformizationRate += storm::utility::one(); // Todo: try other values such as *=1.01 + uniformizationRate *= (storm::utility::one() + storm::utility::convertNumber(env.solver().lra().getAperiodicFactor())); // Get the transitions of the submodel, that is // * a matrix aMarkovian with all (uniformized) transitions from Markovian mec states to all Markovian mec states. @@ -1042,24 +1047,35 @@ namespace storm { // start the iterations - ValueType precision = storm::utility::convertNumber(env.solver().minMax().getPrecision()) / uniformizationRate; - bool relative = env.solver().minMax().getRelativeTerminationCriterion(); + ValueType precision = storm::utility::convertNumber(env.solver().lra().getPrecision()) / uniformizationRate; + bool relative = env.solver().lra().getRelativeTerminationCriterion(); std::vector v(aMarkovian.getRowCount(), storm::utility::zero()); std::vector w = v; std::vector x, b; + auto solverEnv = env; std::unique_ptr> solver; if (hasProbabilisticStates) { + if (env.solver().isForceSoundness()) { + // To get correct results, the inner equation systems are solved exactly. + // TODO investigate how an error would propagate + solverEnv.solver().minMax().setMethod(storm::solver::MinMaxMethod::Topological); + solverEnv.solver().topological().setUnderlyingMinMaxMethod(storm::solver::MinMaxMethod::PolicyIteration); + solverEnv.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Topological); + solverEnv.solver().topological().setUnderlyingEquationSolverType(storm::solver::EquationSolverType::Eigen); + solverEnv.solver().eigen().setMethod(storm::solver::EigenLinearEquationSolverMethod::SparseLU); + } + x.resize(aProbabilistic.getRowGroupCount(), storm::utility::zero()); b = probabilisticChoiceRewards; // Check for requirements of the solver. // The solution is unique as we assume non-zeno MAs. storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxLinearEquationSolverFactory; - storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(env, true, true, dir); + storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(solverEnv, true, true, dir); requirements.clearLowerBounds(); STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); - solver = minMaxLinearEquationSolverFactory.create(env, std::move(aProbabilistic)); + solver = minMaxLinearEquationSolverFactory.create(solverEnv, std::move(aProbabilistic)); solver->setLowerBound(storm::utility::zero()); solver->setHasUniqueSolution(true); solver->setHasNoEndComponents(true); @@ -1067,10 +1083,16 @@ namespace storm { solver->setCachingEnabled(true); } - while (true) { + uint64_t iter = 0; + boost::optional maxIter; + if (env.solver().lra().isMaximalIterationCountSet()) { + maxIter = env.solver().lra().getMaximalIterationCount(); + } + while (!maxIter.is_initialized() || iter < maxIter.get()) { + ++iter; // Compute the expected total rewards for the probabilistic states if (hasProbabilisticStates) { - solver->solveEquations(env, dir, x, b); + solver->solveEquations(solverEnv, dir, x, b); } // now compute the values for the markovian states. We also keep track of the maximal and minimal difference between two values (for convergence checking) auto vIt = v.begin(); @@ -1106,8 +1128,12 @@ namespace storm { storm::utility::vector::addVectors(b, probabilisticChoiceRewards, b); } } + if (maxIter.is_initialized() && iter == maxIter.get()) { + STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations."); + } else { + STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations."); + } return v.front() * uniformizationRate; - } template std::vector SparseMarkovAutomatonCslHelper::computeBoundedUntilProbabilities(Environment const& env, OptimizationDirection dir, storm::storage::SparseMatrix const& transitionMatrix, std::vector const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::pair const& boundsPair); diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index f970244d1..0bcc31726 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -38,6 +38,7 @@ #include "storm/transformer/EndComponentEliminator.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/environment/solver/LongRunAverageSolverEnvironment.h" #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidPropertyException.h" @@ -1251,8 +1252,10 @@ namespace storm { auto underlyingSolverEnvironment = env; if (env.solver().isForceSoundness()) { - // For sound computations, the error in the MECS plus the error in the remaining system should be less then the user defined precsion. - underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().minMax().getPrecision() / storm::utility::convertNumber(2)); + // For sound computations, the error in the MECS plus the error in the remaining system should be less than the user defined precsion. + underlyingSolverEnvironment.solver().lra().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber(2)); + underlyingSolverEnvironment.solver().minMax().setPrecision(env.solver().lra().getPrecision() / storm::utility::convertNumber(2)); + underlyingSolverEnvironment.solver().minMax().setRelativeTerminationCriterion(env.solver().lra().getRelativeTerminationCriterion()); } for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { @@ -1473,11 +1476,11 @@ namespace storm { } // Solve MEC with the method specified in the settings - storm::solver::LraMethod method = env.solver().minMax().getLraMethod(); - if (storm::NumberTraits::IsExact && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) { + storm::solver::LraMethod method = env.solver().lra().getNondetLraMethod(); + if (storm::NumberTraits::IsExact && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::LinearProgramming) { STORM_LOG_INFO("Selecting 'LP' as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method."); method = storm::solver::LraMethod::LinearProgramming; - } else if (env.solver().isForceSoundness() && env.solver().minMax().isLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { + } else if (env.solver().isForceSoundness() && env.solver().lra().isNondetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { STORM_LOG_INFO("Selecting 'VI' as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method."); method = storm::solver::LraMethod::ValueIteration; } @@ -1518,7 +1521,7 @@ namespace storm { std::vector choiceRewards; choiceRewards.reserve(mecChoices.getNumberOfSetBits()); uint64_t currRow = 0; - ValueType selfLoopProb = storm::utility::convertNumber(0.1); // todo try other values + ValueType selfLoopProb = storm::utility::convertNumber(env.solver().lra().getAperiodicFactor()); ValueType scalingFactor = storm::utility::one() - selfLoopProb; for (auto const& mecState : mecStates) { mecTransitionBuilder.newRowGroup(currRow); @@ -1553,14 +1556,21 @@ namespace storm { STORM_LOG_ASSERT(mecTransitions.isProbabilistic(), "The MEC-Matrix is not probabilistic."); // start the iterations - ValueType precision = storm::utility::convertNumber(env.solver().minMax().getPrecision()) / scalingFactor; - bool relative = env.solver().minMax().getRelativeTerminationCriterion(); + ValueType precision = storm::utility::convertNumber(env.solver().lra().getPrecision()) / scalingFactor; + bool relative = env.solver().lra().getRelativeTerminationCriterion(); std::vector x(mecTransitions.getRowGroupCount(), storm::utility::zero()); std::vector xPrime = x; auto multiplier = storm::solver::MultiplierFactory().create(env, mecTransitions); ValueType maxDiff, minDiff; - while (true) { + + uint64_t iter = 0; + boost::optional maxIter; + if (env.solver().lra().isMaximalIterationCountSet()) { + maxIter = env.solver().lra().getMaximalIterationCount(); + } + while (!maxIter.is_initialized() || iter < maxIter.get()) { + ++iter; // Compute the obtained rewards for the next step multiplier->multiplyAndReduce(env, dir, x, &choiceRewards, x); @@ -1585,6 +1595,11 @@ namespace storm { break; } } + if (maxIter.is_initialized() && iter == maxIter.get()) { + STORM_LOG_WARN("LRA computation did not converge within " << iter << " iterations."); + } else { + STORM_LOG_TRACE("LRA computation converged after " << iter << " iterations."); + } if (scheduler) { std::vector localMecChoices(mecTransitions.getRowGroupCount(), 0);