From c674de5893b5d7d19c667e4079c573490d005afc Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Tue, 11 Aug 2020 16:53:05 +0200 Subject: [PATCH] Deterministic infinite horizon: Added gain-bias and lra-distribution based solution methods --- ...arseDeterministicInfiniteHorizonHelper.cpp | 260 +++++++++++++++++- ...SparseDeterministicInfiniteHorizonHelper.h | 12 +- 2 files changed, 259 insertions(+), 13 deletions(-) diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp index d6dc8eb3c..f7327d8c0 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp @@ -17,9 +17,10 @@ #include "storm/utility/vector.h" #include "storm/environment/solver/LongRunAverageSolverEnvironment.h" -#include "storm/environment/solver/MinMaxSolverEnvironment.h" +#include "storm/environment/solver/TopologicalSolverEnvironment.h" #include "storm/exceptions/UnmetRequirementException.h" +#include "storm/exceptions/NotSupportedException.h" namespace storm { namespace modelchecker { @@ -46,12 +47,12 @@ namespace storm { } template - ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) { + ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForComponent(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, storm::storage::StronglyConnectedComponent const& component) { // For deterministic models, we compute the LRA for a BSCC STORM_LOG_ASSERT(!this->isProduceSchedulerSet(), "Scheduler production enabled for deterministic model."); - auto trivialResult = computeLraForTrivialBscc(env, stateRewardsGetter, actionRewardsGetter, component); + auto trivialResult = computeLraForTrivialBscc(env, stateValueGetter, actionValueGetter, component); if (trivialResult.first) { return trivialResult.second; } @@ -67,18 +68,18 @@ namespace storm { } STORM_LOG_TRACE("Computing LRA for BSCC of size " << component.size() << " using '" << storm::solver::toString(method) << "'."); if (method == storm::solver::LraMethod::ValueIteration) { - return computeLraForBsccVi(env, stateRewardsGetter, actionRewardsGetter, component); - }/* else if (method == storm::solver::LraMethod::LraDistributionEquations) { + return computeLraForBsccVi(env, stateValueGetter, actionValueGetter, component); + } else if (method == storm::solver::LraMethod::LraDistributionEquations) { // We only need the first element of the pair as the lra distribution is not relevant at this point. - return computeLongRunAveragesForBsccLraDistr(env, bscc, rateMatrix, valueGetter, exitRateVector).first; + return computeLraForBsccLraDistr(env, stateValueGetter, actionValueGetter, component).first; } STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to " << storm::solver::toString(storm::solver::LraMethod::GainBiasEquations) << "."); // We don't need the bias values - return computeLongRunAveragesForBsccGainBias(env, bscc, rateMatrix, valueGetter, exitRateVector).first;*/ + return computeLraForBsccGainBias(env, stateValueGetter, actionValueGetter, component).first; } template - std::pair SparseDeterministicInfiniteHorizonHelper::computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) { + std::pair SparseDeterministicInfiniteHorizonHelper::computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, storm::storage::StronglyConnectedComponent const& component) { // For deterministic models, we can catch the case where all values are the same. This includes the special case where the BSCC consist only of just one state. bool first = true; @@ -86,8 +87,9 @@ namespace storm { for (auto const& element : component) { auto state = internal::getComponentElementState(element); STORM_LOG_ASSERT(state == *internal::getComponentElementChoicesBegin(element), "Unexpected choice index at state " << state << " of deterministic model."); - ValueType curr = stateRewardsGetter(state) + (this->isContinuousTime() ? (*this->_exitRates)[state] * actionRewardsGetter(state) : actionRewardsGetter(state)); + ValueType curr = stateValueGetter(state) + (this->isContinuousTime() ? (*this->_exitRates)[state] * actionValueGetter(state) : actionValueGetter(state)); if (first) { + val = curr; first = false; } else if (val != curr) { return {false, storm::utility::zero()}; @@ -98,8 +100,12 @@ namespace storm { } + template <> + storm::RationalFunction SparseDeterministicInfiniteHorizonHelper::computeLraForBsccVi(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, storm::storage::StronglyConnectedComponent const& bscc) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The requested Method for LRA computation is not supported for parametric models."); + } template - ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForBsccVi(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& bscc) { + ValueType SparseDeterministicInfiniteHorizonHelper::computeLraForBsccVi(Environment const& env, ValueGetter const& stateValueGetter, ValueGetter const& actionValueGetter, storm::storage::StronglyConnectedComponent const& bscc) { // Collect parameters of the computation ValueType aperiodicFactor = storm::utility::convertNumber(env.solver().lra().getAperiodicFactor()); @@ -108,14 +114,243 @@ namespace storm { if (this->isContinuousTime()) { // We assume a CTMC (with deterministic timed states and no instant states) storm::modelchecker::helper::internal::LraViHelper viHelper(bscc, this->_transitionMatrix, aperiodicFactor, this->_markovianStates, this->_exitRates); - return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter, this->_exitRates); + return viHelper.performValueIteration(env, stateValueGetter, actionValueGetter, this->_exitRates); } else { // We assume a DTMC (with deterministic timed states and no instant states) storm::modelchecker::helper::internal::LraViHelper viHelper(bscc, this->_transitionMatrix, aperiodicFactor); - return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter); + return viHelper.performValueIteration(env, stateValueGetter, actionValueGetter); } } + template + std::pair> SparseDeterministicInfiniteHorizonHelper::computeLraForBsccGainBias(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) { + // We build the equation system as in Line 3 of Algorithm 3 from + // Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017) + // https://doi.org/10.1007/978-3-319-68167-2_25 + // The first variable corresponds to the gain of the bscc whereas the subsequent variables yield the bias for each state s_1, s_2, .... + // No bias variable for s_0 is needed since it is always set to zero, yielding an nxn equation system matrix + // To make this work for CTMC, we could uniformize the model. This preserves LRA and ensures that we can compute the + // LRA as for a DTMC (the soujourn time in each state is the same). If we then multiply the equations with the uniformization rate, + // the uniformization rate cancels out. Hence, we obtain the equation system below. + + // Get a mapping from global state indices to local ones. + std::unordered_map toLocalIndexMap; + uint64_t localIndex = 0; + for (auto const& globalIndex : bscc) { + toLocalIndexMap[globalIndex] = localIndex; + ++localIndex; + } + + // Prepare an environment for the underlying equation solver + auto subEnv = env; + if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) { + // Topological solver does not make any sense since the BSCC is connected. + subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType(), subEnv.solver().topological().isUnderlyingEquationSolverTypeSetFromDefault()); + } + subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); + + // Build the equation system matrix and vector. + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + storm::storage::SparseMatrixBuilder builder(bscc.size(), bscc.size()); + std::vector eqSysVector; + eqSysVector.reserve(bscc.size()); + // The first row asserts that the weighted bias variables and the reward at s_0 sum up to the gain + uint64_t row = 0; + ValueType entryValue; + for (auto const& globalState : bscc) { + ValueType rateAtState = this->_exitRates ? (*this->_exitRates)[globalState] : storm::utility::one(); + // Coefficient for the gain variable + if (isEquationSystemFormat) { + // '1-0' in row 0 and -(-1) in other rows + builder.addNextValue(row, 0, storm::utility::one()); + } else if (row > 0) { + // No coeficient in row 0, othwerise substract the gain + builder.addNextValue(row, 0, -storm::utility::one()); + } + // Compute weighted sum over successor state. As this is a BSCC, each successor state will again be in the BSCC. + auto diagonalValue = storm::utility::zero(); + if (row > 0) { + if (isEquationSystemFormat) { + diagonalValue = rateAtState; + } else { + diagonalValue = storm::utility::one() - rateAtState; + } + } + bool needDiagonalEntry = !storm::utility::isZero(diagonalValue); + for (auto const& entry : this->_transitionMatrix.getRow(globalState)) { + uint64_t col = toLocalIndexMap[entry.getColumn()]; + if (col == 0) { + //Skip transition to state_0. This corresponds to setting the bias of state_0 to zero + continue; + } + entryValue = entry.getValue() * rateAtState; + if (isEquationSystemFormat) { + entryValue = -entryValue; + } + if (needDiagonalEntry && col >= row) { + if (col == row) { + entryValue += diagonalValue; + } else { // col > row + builder.addNextValue(row, row, diagonalValue); + } + needDiagonalEntry = false; + } + builder.addNextValue(row, col, entryValue); + } + if (needDiagonalEntry) { + builder.addNextValue(row, row, diagonalValue); + } + + eqSysVector.push_back(stateValuesGetter(globalState) + rateAtState * actionValuesGetter(globalState)); + ++row; + } + + // Create a linear equation solver + auto solver = linearEquationSolverFactory.create(subEnv, builder.build()); + // Check solver requirements. + auto requirements = solver->getRequirements(subEnv); + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + // Todo: Find bounds on the bias variables. Just inserting the maximal value from the vector probably does not work. + + std::vector eqSysSol(bscc.size(), storm::utility::zero()); + // Take the mean of the rewards as an initial guess for the gain + //eqSysSol.front() = std::accumulate(eqSysVector.begin(), eqSysVector.end(), storm::utility::zero()) / storm::utility::convertNumber(bscc.size()); + solver->solveEquations(subEnv, eqSysSol, eqSysVector); + + ValueType gain = eqSysSol.front(); + // insert bias value for state 0 + eqSysSol.front() = storm::utility::zero(); + // Return the gain and the bias values + return std::pair>(std::move(gain), std::move(eqSysSol)); + } + + template + std::pair> SparseDeterministicInfiniteHorizonHelper::computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc) { + + // Let A be ab auxiliary Matrix with A[s,s] = R(s,s) - r(s) & A[s,s'] = R(s,s') for s,s' in BSCC and s!=s'. + // We build and solve the equation system for + // x*A=0 & x_0+...+x_n=1 <=> A^t*x=0=x-x & x_0+...+x_n=1 <=> (1+A^t)*x = x & 1-x_0-...-x_n-1=x_n + // Then, x[i] will be the fraction of the time we are in state i. + + // This method assumes that this BSCC consist of more than one state + if (bscc.size() == 1) { + ValueType lraValue = stateValuesGetter(*bscc.begin()) + (this->isContinuousTime() ? (*this->_exitRates)[*bscc.begin()] * actionValuesGetter(*bscc.begin()) : actionValuesGetter(*bscc.begin())); + return { lraValue, {storm::utility::one()} }; + } + + // Prepare an environment for the underlying linear equation solver + auto subEnv = env; + if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) { + // Topological solver does not make any sense since the BSCC is connected. + subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType(), subEnv.solver().topological().isUnderlyingEquationSolverTypeSetFromDefault()); + } + subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); + + // Get a mapping from global state indices to local ones as well as a bitvector containing states within the BSCC. + std::unordered_map toLocalIndexMap; + storm::storage::BitVector bsccStates(this->_transitionMatrix.getRowCount(), false); + uint64_t localIndex = 0; + for (auto const& globalIndex : bscc) { + bsccStates.set(globalIndex, true); + toLocalIndexMap[globalIndex] = localIndex; + ++localIndex; + } + + // Build the auxiliary Matrix A. + auto auxMatrix = this->_transitionMatrix.getSubmatrix(false, bsccStates, bsccStates, true); // add diagonal entries! + uint64_t row = 0; + for (auto const& globalIndex : bscc) { + ValueType rateAtState = this->_exitRates ? (*this->_exitRates)[globalIndex] : storm::utility::one(); + for (auto& entry : auxMatrix.getRow(row)) { + if (entry.getColumn() == row) { + // This value is non-zero since we have a BSCC with more than one state + entry.setValue(rateAtState * (entry.getValue() - storm::utility::one())); + } else if (this->isContinuousTime()) { + entry.setValue(entry.getValue() * rateAtState); + } + } + ++row; + } + assert(row == auxMatrix.getRowCount()); + + // We need to consider A^t. This will not delete diagonal entries since they are non-zero. + auxMatrix = auxMatrix.transpose(); + + // Check whether we need the fixpoint characterization + storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; + bool isFixpointFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem; + if (isFixpointFormat) { + // Add a 1 on the diagonal + for (row = 0; row < auxMatrix.getRowCount(); ++row) { + for (auto& entry : auxMatrix.getRow(row)) { + if (entry.getColumn() == row) { + entry.setValue(storm::utility::one() + entry.getValue()); + } + } + } + } + + // We now build the equation system matrix. + // We can drop the last row of A and add ones in this row instead to assert that the variables sum up to one + // Phase 1: replace the existing entries of the last row with ones + uint64_t col = 0; + uint64_t lastRow = auxMatrix.getRowCount() - 1; + for (auto& entry : auxMatrix.getRow(lastRow)) { + entry.setColumn(col); + if (isFixpointFormat) { + if (col == lastRow) { + entry.setValue(storm::utility::zero()); + } else { + entry.setValue(-storm::utility::one()); + } + } else { + entry.setValue(storm::utility::one()); + } + ++col; + } + storm::storage::SparseMatrixBuilder builder(std::move(auxMatrix)); + for (; col <= lastRow; ++col) { + if (isFixpointFormat) { + if (col != lastRow) { + builder.addNextValue(lastRow, col, -storm::utility::one()); + } + } else { + builder.addNextValue(lastRow, col, storm::utility::one()); + } + } + + std::vector bsccEquationSystemRightSide(bscc.size(), storm::utility::zero()); + bsccEquationSystemRightSide.back() = storm::utility::one(); + + // Create a linear equation solver + auto solver = linearEquationSolverFactory.create(subEnv, builder.build()); + solver->setBounds(storm::utility::zero(), storm::utility::one()); + // Check solver requirements. + auto requirements = solver->getRequirements(subEnv); + requirements.clearLowerBounds(); + requirements.clearUpperBounds(); + STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UnmetRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); + + std::vector lraDistr(bscc.size(), storm::utility::one() / storm::utility::convertNumber(bscc.size())); + solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide); + + // Calculate final LRA Value + ValueType result = storm::utility::zero(); + auto solIt = lraDistr.begin(); + for (auto const& globalState : bscc) { + if (this->isContinuousTime()) { + result += (*solIt) * (stateValuesGetter(globalState) + (*this->_exitRates)[globalState] * actionValuesGetter(globalState)); + } else { + result += (*solIt) * (stateValuesGetter(globalState) + actionValuesGetter(globalState)); + } + ++solIt; + } + assert(solIt == lraDistr.end()); + + return std::pair>(std::move(result), std::move(lraDistr)); + } + template std::pair, std::vector> SparseDeterministicInfiniteHorizonHelper::buildSspMatrixVector(std::vector const& bsccLraValues, std::vector const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem) { @@ -203,6 +438,7 @@ namespace storm { template class SparseDeterministicInfiniteHorizonHelper; template class SparseDeterministicInfiniteHorizonHelper; + template class SparseDeterministicInfiniteHorizonHelper; } } diff --git a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h index 5c22cc1af..b5a78fad1 100644 --- a/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h +++ b/src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h @@ -46,10 +46,20 @@ namespace storm { std::pair computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); /*! - * As computeLraForMec but uses value iteration as a solution method (independent of what is set in env) + * As computeLraForComponent but uses value iteration as a solution method (independent of what is set in env) */ ValueType computeLraForBsccVi(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + /*! + * As computeLraForComponent but solves a linear equation system encoding gain and bias (independent of what is set in env) + * @see Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017), https://doi.org/10.1007/978-3-319-68167-2_25 + */ + std::pair> computeLraForBsccGainBias(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + /*! + * As computeLraForComponent but solves a linear equation system consisting encoding the long run average (steady state) distribution (independent of what is set in env) + */ + std::pair> computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc); + std::pair, std::vector> buildSspMatrixVector(std::vector const& bsccLraValues, std::vector const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem); /*!