diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp index 18c4f268e..cccc2bbe2 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp @@ -1,6 +1,7 @@ #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h" #include +#include #include "src/utility/constants.h" #include "src/utility/macros.h" @@ -10,10 +11,12 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/solver/LpSolver.h" + #include "src/exceptions/InvalidPropertyException.h" +#include "src/storage/expressions/Expressions.h" #include "src/storage/MaximalEndComponentDecomposition.h" -#include "src/storage/MaximalEndComponent.h" namespace storm { namespace modelchecker { @@ -320,28 +323,28 @@ namespace storm { template std::vector SparseMdpPrctlModelChecker::computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const { // If there are no goal states, we avoid the computation and directly return zero. + auto numOfStates = this->getModel().getNumberOfStates(); if (psiStates.empty()) { - return std::vector(model.getNumberOfStates(), storm::utility::zero()); + return std::vector(numOfStates, storm::utility::zero()); } // Likewise, if all bits are set, we can avoid the computation and set. if ((~psiStates).empty()) { - return std::vector(model.getNumberOfStates(), storm::utility::one()); + return std::vector(numOfStates, storm::utility::one()); } // Start by decomposing the MDP into its MECs. - storm::storage::MaximalEndComponentDecomposition mecDecomposition(model); + storm::storage::MaximalEndComponentDecomposition mecDecomposition(this->getModelAs>()); // Get some data members for convenience. - typename storm::storage::SparseMatrix const& transitionMatrix = model.getTransitionMatrix(); - std::vector const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices(); + typename storm::storage::SparseMatrix const& transitionMatrix = this->getModel().getTransitionMatrix(); // Now start with compute the long-run average for all end components in isolation. std::vector lraValuesForEndComponents; // While doing so, we already gather some information for the following steps. - std::vector stateToMecIndexMap(model.getNumberOfStates()); - storm::storage::BitVector statesInMecs(model.getNumberOfStates()); + std::vector stateToMecIndexMap(numOfStates); + storm::storage::BitVector statesInMecs(numOfStates); for (uint_fast64_t currentMecIndex = 0; currentMecIndex < mecDecomposition.size(); ++currentMecIndex) { storm::storage::MaximalEndComponent const& mec = mecDecomposition[currentMecIndex]; @@ -362,7 +365,7 @@ namespace storm { storm::storage::BitVector statesNotContainedInAnyMec = ~statesInMecs; // Prepare result vector. - std::vector result(model.getNumberOfStates()); + std::vector result(numOfStates); //Set the values for all states in MECs. for (auto state : statesInMecs) { @@ -374,10 +377,10 @@ namespace storm { for (auto state : statesNotContainedInAnyMec){ //calculate what result values the reachable states in MECs have - storm::storage::BitVector currentState(model.getNumberOfStates); + storm::storage::BitVector currentState(numOfStates); currentState.set(state); storm::storage::BitVector reachableStates = storm::utility::graph::getReachableStates( - transitionMatrix, currentState, storm::storage::BitVector(model.getNumberOfStates, true), statesInMecs + transitionMatrix, currentState, storm::storage::BitVector(numOfStates, true), statesInMecs ); storm::storage::BitVector reachableMecStates = statesInMecs & reachableStates; @@ -406,62 +409,47 @@ namespace storm { } template - ValueType SparseMarkovAutomatonCslModelChecker::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { + ValueType SparseMdpPrctlModelChecker::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec) { std::shared_ptr solver = storm::utility::solver::getLpSolver("LRA for MEC"); solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize); //// First, we need to create the variables for the problem. - //std::map stateToVariableMap; - //for (auto const& stateChoicesPair : mec) { - // std::string variableName = "x" + std::to_string(stateChoicesPair.first); - // stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName); - //} - //storm::expressions::Variable k = solver->addUnboundedContinuousVariable("k", 1); - //solver->update(); - - //// Now we encode the problem as constraints. - //for (auto const& stateChoicesPair : mec) { - // uint_fast64_t state = stateChoicesPair.first; - - // // Now, based on the type of the state, create a suitable constraint. - // if (markovianStates.get(state)) { - // storm::expressions::Expression constraint = stateToVariableMap.at(state); - - // for (auto element : transitionMatrix.getRow(nondeterministicChoiceIndices[state])) { - // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); - // } - - // constraint = constraint + solver->getConstant(storm::utility::one() / exitRates[state]) * k; - // storm::expressions::Expression rightHandSide = psiStates.get(state) ? solver->getConstant(storm::utility::one() / exitRates[state]) : solver->getConstant(0); - // if (minimize) { - // constraint = constraint <= rightHandSide; - // } else { - // constraint = constraint >= rightHandSide; - // } - // solver->addConstraint("state" + std::to_string(state), constraint); - // } else { - // // For probabilistic states, we want to add the constraint x_s <= sum P(s, a, s') * x_s' where a is the current action - // // and the sum ranges over all states s'. - // for (auto choice : stateChoicesPair.second) { - // storm::expressions::Expression constraint = stateToVariableMap.at(state); - - // for (auto element : transitionMatrix.getRow(choice)) { - // constraint = constraint - stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); - // } - - // storm::expressions::Expression rightHandSide = solver->getConstant(storm::utility::zero()); - // if (minimize) { - // constraint = constraint <= rightHandSide; - // } else { - // constraint = constraint >= rightHandSide; - // } - // solver->addConstraint("state" + std::to_string(state), constraint); - // } - // } - //} - - //solver->optimize(); - //return solver->getContinuousValue(k); + std::map stateToVariableMap; + for (auto const& stateChoicesPair : mec) { + std::string variableName = "h" + std::to_string(stateChoicesPair.first); + stateToVariableMap[stateChoicesPair.first] = solver->addUnboundedContinuousVariable(variableName); + } + storm::expressions::Variable lambda = solver->addUnboundedContinuousVariable("L", 1); + solver->update(); + + // Now we encode the problem as constraints. + for (auto const& stateChoicesPair : mec) { + uint_fast64_t state = stateChoicesPair.first; + + // Now, based on the type of the state, create a suitable constraint. + for (auto choice : stateChoicesPair.second) { + storm::expressions::Expression constraint = solver->getConstant(1); + ValueType w = 0; + + for (auto element : transitionMatrix.getRow(choice)) { + constraint = constraint + stateToVariableMap.at(element.getColumn()) * solver->getConstant(element.getValue()); + if (psiStates.get(element.getColumn())) { + w += element.getValue(); + } + } + constraint = constraint - solver->getConstant(w) * lambda; + + if (minimize) { + constraint = stateToVariableMap.at(state) <= constraint; + } else { + constraint = stateToVariableMap.at(state) >= constraint; + } + solver->addConstraint("state" + std::to_string(state) + "," + std::to_string(choice), constraint); + } + } + + solver->optimize(); + return solver->getContinuousValue(lambda); } template diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h index 3dad8e619..9718457cd 100644 --- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h @@ -5,6 +5,8 @@ #include "src/models/Mdp.h" #include "src/utility/solver.h" #include "src/solver/NondeterministicLinearEquationSolver.h" +#include "src/storage/MaximalEndComponent.h" + namespace storm { namespace counterexamples {