sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Deterministic infinite horizon: Added gain-bias and lra-distribution based solution methods

main
Tim Quatmann 5 years ago
parent
0cc2b1c749
commit
c674de5893
2 changed files with 259 additions and 13 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 260
      src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp
  2. 12
      src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h

260
src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.cpp
View File

@ -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 <typename ValueType>
ValueType SparseDeterministicInfiniteHorizonHelper<ValueType>::computeLraForComponent(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) {
ValueType SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>(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<ValueType>(env, bscc, rateMatrix, valueGetter, exitRateVector).first;*/
return computeLraForBsccGainBias(env, stateValueGetter, actionValueGetter, component).first;
}
template <typename ValueType>
std::pair<bool, ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::computeLraForTrivialBscc(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& component) {
std::pair<bool, ValueType> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>()};
@ -98,8 +100,12 @@ namespace storm {
}
template <>
storm::RationalFunction SparseDeterministicInfiniteHorizonHelper<storm::RationalFunction>::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 <typename ValueType>
ValueType SparseDeterministicInfiniteHorizonHelper<ValueType>::computeLraForBsccVi(Environment const& env, ValueGetter const& stateRewardsGetter, ValueGetter const& actionRewardsGetter, storm::storage::StronglyConnectedComponent const& bscc) {
ValueType SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>(env.solver().lra().getAperiodicFactor());
@ -108,12 +114,241 @@ namespace storm {
if (this->isContinuousTime()) {
// We assume a CTMC (with deterministic timed states and no instant states)
storm::modelchecker::helper::internal::LraViHelper<ValueType, storm::storage::StronglyConnectedComponent, storm::modelchecker::helper::internal::LraViTransitionsType::DetTsNoIs> 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<ValueType, storm::storage::StronglyConnectedComponent, storm::modelchecker::helper::internal::LraViTransitionsType::DetTsNoIs> viHelper(bscc, this->_transitionMatrix, aperiodicFactor);
return viHelper.performValueIteration(env, stateRewardsGetter, actionRewardsGetter);
return viHelper.performValueIteration(env, stateValueGetter, actionValueGetter);
}
}
template <typename ValueType>
std::pair<ValueType, std::vector<ValueType>> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<uint64_t, uint64_t> 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<ValueType> linearEquationSolverFactory;
bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
storm::storage::SparseMatrixBuilder<ValueType> builder(bscc.size(), bscc.size());
std::vector<ValueType> 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<ValueType>();
// Coefficient for the gain variable
if (isEquationSystemFormat) {
// '1-0' in row 0 and -(-1) in other rows
builder.addNextValue(row, 0, storm::utility::one<ValueType>());
} else if (row > 0) {
// No coeficient in row 0, othwerise substract the gain
builder.addNextValue(row, 0, -storm::utility::one<ValueType>());
}
// 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<ValueType>();
if (row > 0) {
if (isEquationSystemFormat) {
diagonalValue = rateAtState;
} else {
diagonalValue = storm::utility::one<ValueType>() - 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<ValueType> eqSysSol(bscc.size(), storm::utility::zero<ValueType>());
// Take the mean of the rewards as an initial guess for the gain
//eqSysSol.front() = std::accumulate(eqSysVector.begin(), eqSysVector.end(), storm::utility::zero<ValueType>()) / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size());
solver->solveEquations(subEnv, eqSysSol, eqSysVector);
ValueType gain = eqSysSol.front();
// insert bias value for state 0
eqSysSol.front() = storm::utility::zero<ValueType>();
// Return the gain and the bias values
return std::pair<ValueType, std::vector<ValueType>>(std::move(gain), std::move(eqSysSol));
}
template <typename ValueType>
std::pair<ValueType, std::vector<ValueType>> SparseDeterministicInfiniteHorizonHelper<ValueType>::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<ValueType>()} };
}
// 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<uint64_t, uint64_t> 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<ValueType>();
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<ValueType>()));
} 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<ValueType> 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<ValueType>() + 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<ValueType>());
} else {
entry.setValue(-storm::utility::one<ValueType>());
}
} else {
entry.setValue(storm::utility::one<ValueType>());
}
++col;
}
storm::storage::SparseMatrixBuilder<ValueType> builder(std::move(auxMatrix));
for (; col <= lastRow; ++col) {
if (isFixpointFormat) {
if (col != lastRow) {
builder.addNextValue(lastRow, col, -storm::utility::one<ValueType>());
}
} else {
builder.addNextValue(lastRow, col, storm::utility::one<ValueType>());
}
}
std::vector<ValueType> bsccEquationSystemRightSide(bscc.size(), storm::utility::zero<ValueType>());
bsccEquationSystemRightSide.back() = storm::utility::one<ValueType>();
// Create a linear equation solver
auto solver = linearEquationSolverFactory.create(subEnv, builder.build());
solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
// 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<ValueType> lraDistr(bscc.size(), storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size()));
solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide);
// Calculate final LRA Value
ValueType result = storm::utility::zero<ValueType>();
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<ValueType, std::vector<ValueType>>(std::move(result), std::move(lraDistr));
}
template <typename ValueType>
@ -203,6 +438,7 @@ namespace storm {
template class SparseDeterministicInfiniteHorizonHelper<double>;
template class SparseDeterministicInfiniteHorizonHelper<storm::RationalNumber>;
template class SparseDeterministicInfiniteHorizonHelper<storm::RationalFunction>;
}
}

12
src/storm/modelchecker/helper/infinitehorizon/SparseDeterministicInfiniteHorizonHelper.h
View File

@ -46,10 +46,20 @@ namespace storm {
std::pair<bool, ValueType> 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<ValueType, std::vector<ValueType>> 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<ValueType, std::vector<ValueType>> computeLraForBsccLraDistr(Environment const& env, ValueGetter const& stateValuesGetter, ValueGetter const& actionValuesGetter, storm::storage::StronglyConnectedComponent const& bscc);
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> buildSspMatrixVector(std::vector<ValueType> const& bsccLraValues, std::vector<uint64_t> const& inputStateToBsccIndexMap, storm::storage::BitVector const& statesNotInComponent, bool asEquationSystem);
/*!

Write Preview
Loading…
Cancel
Save