Storm 1.10.0.1
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
DeterministicSchedsObjectiveHelper.cpp
Go to the documentation of this file.
2
21#include "storm/utility/graph.h"
23
26
28
29template<typename ModelType>
32 : model(model), objective(objective) {
33 initialize();
34}
35
36template<typename ModelType>
39 auto checkResult = mc.check(formula);
40 STORM_LOG_THROW(checkResult && checkResult->isExplicitQualitativeCheckResult(), storm::exceptions::UnexpectedException,
41 "Unexpected type of check result for subformula " << formula << ".");
42 return checkResult->asExplicitQualitativeCheckResult().getTruthValuesVector();
43}
44
45template<typename ValueType>
47 if (formula.isRewardOperatorFormula()) {
48 boost::optional<std::string> rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName();
50 rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel();
51
52 // Get a reward model where the state rewards are scaled accordingly
53 std::vector<ValueType> stateRewardWeights(model.getNumberOfStates(), storm::utility::zero<ValueType>());
54 for (auto const markovianState : model.getMarkovianStates()) {
55 stateRewardWeights[markovianState] = storm::utility::one<ValueType>() / model.getExitRate(markovianState);
56 }
57 return rewardModel.getTotalActionRewardVector(model.getTransitionMatrix(), stateRewardWeights);
58 } else {
59 assert(formula.isTimeOperatorFormula());
60 std::vector<ValueType> result(model.getNumberOfChoices(), storm::utility::zero<ValueType>());
61 for (auto const markovianState : model.getMarkovianStates()) {
62 for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(markovianState)) {
63 result[choice] = storm::utility::one<ValueType>() / model.getExitRate(markovianState);
64 }
65 }
66 return result;
67 }
68}
69
70template<typename ValueType>
71std::vector<ValueType> getTotalRewardVector(storm::models::sparse::Mdp<ValueType> const& model, storm::logic::Formula const& formula) {
72 if (formula.isRewardOperatorFormula()) {
73 boost::optional<std::string> rewardModelName = formula.asRewardOperatorFormula().getOptionalRewardModelName();
75 rewardModelName.is_initialized() ? model.getRewardModel(rewardModelName.get()) : model.getUniqueRewardModel();
76 return rewardModel.getTotalRewardVector(model.getTransitionMatrix());
77 } else {
78 assert(formula.isTimeOperatorFormula());
79 return std::vector<ValueType>(model.getNumberOfChoices(), storm::utility::one<ValueType>());
80 }
81}
82
83template<typename ModelType>
84void DeterministicSchedsObjectiveHelper<ModelType>::initialize() {
85 STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1, "Expected a single initial state.");
86 uint64_t initialState = *model.getInitialStates().begin();
87 relevantZeroRewardChoices = storm::storage::BitVector(model.getTransitionMatrix().getRowCount(), false);
88 auto const& formula = *objective.formula;
89 if (formula.isProbabilityOperatorFormula() && formula.getSubformula().isUntilFormula()) {
90 storm::storage::BitVector phiStates = evaluatePropositionalFormula(model, formula.getSubformula().asUntilFormula().getLeftSubformula());
91 storm::storage::BitVector psiStates = evaluatePropositionalFormula(model, formula.getSubformula().asUntilFormula().getRightSubformula());
92 auto backwardTransitions = model.getBackwardTransitions();
93 auto prob1States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), backwardTransitions,
94 phiStates, psiStates);
95 auto prob0States = storm::utility::graph::performProb0A(backwardTransitions, phiStates, psiStates);
96 if (prob0States.get(initialState)) {
97 constantInitialStateValue = storm::utility::zero<ValueType>();
98 } else if (prob1States.get(initialState)) {
99 constantInitialStateValue = storm::utility::one<ValueType>();
100 }
101 maybeStates = ~(prob0States | prob1States);
102 // Cut away those states that are not reachable, or only reachable via non-maybe states.
103 maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates);
104 for (auto const& state : maybeStates) {
105 for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) {
106 auto rowSum = model.getTransitionMatrix().getConstrainedRowSum(choice, prob1States);
107 if (storm::utility::isZero(rowSum)) {
108 relevantZeroRewardChoices.set(choice);
109 } else {
110 choiceRewards.emplace(choice, rowSum);
111 }
112 }
113 }
114 } else if (formula.getSubformula().isEventuallyFormula() && (formula.isRewardOperatorFormula() || formula.isTimeOperatorFormula())) {
115 storm::storage::BitVector rew0States = evaluatePropositionalFormula(model, formula.getSubformula().asEventuallyFormula().getSubformula());
116 if (formula.isRewardOperatorFormula()) {
117 auto const& baseRewardModel = formula.asRewardOperatorFormula().hasRewardModelName()
118 ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName())
119 : model.getUniqueRewardModel();
120 auto rewardModel =
121 storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), formula.getSubformula().asEventuallyFormula());
122 storm::storage::BitVector statesWithoutReward = rewardModel.get().getStatesWithZeroReward(model.getTransitionMatrix());
123 rew0States = storm::utility::graph::performProb1A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(),
124 model.getBackwardTransitions(), statesWithoutReward, rew0States);
125 }
126 if (rew0States.get(initialState)) {
127 constantInitialStateValue = storm::utility::zero<ValueType>();
128 }
129 maybeStates = ~rew0States;
130 // Cut away those states that are not reachable, or only reachable via non-maybe states.
131 maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates);
132 std::vector<ValueType> choiceBasedRewards = getTotalRewardVector(model, *objective.formula);
133 for (auto const& state : maybeStates) {
134 for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) {
135 auto const& value = choiceBasedRewards[choice];
136 if (storm::utility::isZero(value)) {
137 relevantZeroRewardChoices.set(choice);
138 } else {
139 choiceRewards.emplace(choice, value);
140 }
141 }
142 }
143 } else if (formula.isRewardOperatorFormula() && formula.getSubformula().isTotalRewardFormula()) {
144 auto const& baseRewardModel = formula.asRewardOperatorFormula().hasRewardModelName()
145 ? model.getRewardModel(formula.asRewardOperatorFormula().getRewardModelName())
146 : model.getUniqueRewardModel();
147 auto rewardModel =
148 storm::utility::createFilteredRewardModel(baseRewardModel, model.isDiscreteTimeModel(), formula.getSubformula().asTotalRewardFormula());
149 storm::storage::BitVector statesWithoutReward = rewardModel.get().getStatesWithZeroReward(model.getTransitionMatrix());
150 storm::storage::BitVector rew0States =
151 storm::utility::graph::performProbGreater0E(model.getBackwardTransitions(), statesWithoutReward, ~statesWithoutReward);
152 rew0States.complement();
153 if (rew0States.get(initialState)) {
154 constantInitialStateValue = storm::utility::zero<ValueType>();
155 }
156 maybeStates = ~rew0States;
157 // Cut away those states that are not reachable, or only reachable via non-maybe states.
158 maybeStates &= storm::utility::graph::getReachableStates(model.getTransitionMatrix(), model.getInitialStates(), maybeStates, ~maybeStates);
159 std::vector<ValueType> choiceBasedRewards = getTotalRewardVector(model, *objective.formula);
160 for (auto const& state : maybeStates) {
161 for (auto const& choice : model.getTransitionMatrix().getRowGroupIndices(state)) {
162 auto const& value = choiceBasedRewards[choice];
163 if (storm::utility::isZero(value)) {
164 relevantZeroRewardChoices.set(choice);
165 } else {
166 choiceRewards.emplace(choice, value);
167 }
168 }
169 }
170 } else {
171 STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The given formula " << formula << " is not supported.");
172 }
173
174 // negate choice rewards for minimizing objectives
175 if (storm::solver::minimize(formula.getOptimalityType())) {
176 for (auto& entry : choiceRewards) {
177 entry.second *= -storm::utility::one<ValueType>();
178 }
179 }
180
181 // Find out if we have to deal with infinite positive or negative rewards
182 storm::storage::BitVector negativeRewardChoices(model.getNumberOfChoices(), false);
183 storm::storage::BitVector positiveRewardChoices(model.getNumberOfChoices(), false);
184 for (auto const& rew : getChoiceRewards()) {
185 if (rew.second > storm::utility::zero<ValueType>()) {
186 positiveRewardChoices.set(rew.first, true);
187 } else {
188 assert(rew.second < storm::utility::zero<ValueType>());
189 negativeRewardChoices.set(rew.first, true);
190 }
191 }
192 auto backwardTransitions = model.getBackwardTransitions();
193 bool hasNegativeEC =
194 storm::utility::graph::checkIfECWithChoiceExists(model.getTransitionMatrix(), backwardTransitions, getMaybeStates(), negativeRewardChoices);
195 bool hasPositiveEc =
196 storm::utility::graph::checkIfECWithChoiceExists(model.getTransitionMatrix(), backwardTransitions, getMaybeStates(), positiveRewardChoices);
197 STORM_LOG_THROW(!(hasNegativeEC && hasPositiveEc), storm::exceptions::NotSupportedException,
198 "Objective is not convergent: Infinite positive and infinite negative reward is possible.");
199 infinityCase = hasNegativeEC ? InfinityCase::HasNegativeInfinite : (hasPositiveEc ? InfinityCase::HasPositiveInfinite : InfinityCase::AlwaysFinite);
200
201 if (infinityCase == InfinityCase::HasNegativeInfinite) {
202 storm::storage::BitVector negativeEcStates(maybeStates.size(), false);
203 storm::storage::MaximalEndComponentDecomposition<ValueType> mecs(model.getTransitionMatrix(), backwardTransitions, getMaybeStates());
204 for (auto const& mec : mecs) {
205 if (std::any_of(mec.begin(), mec.end(), [&negativeRewardChoices](auto const& sc) {
206 return std::any_of(sc.second.begin(), sc.second.end(), [&negativeRewardChoices](auto const& c) { return negativeRewardChoices.get(c); });
207 })) {
208 for (auto const& sc : mec) {
209 negativeEcStates.set(sc.first);
210 }
211 }
212 }
213 STORM_LOG_ASSERT(!negativeEcStates.empty(), "Expected some negative ec");
214 rewMinusInfEStates = storm::utility::graph::performProbGreater0E(backwardTransitions, getMaybeStates(), negativeEcStates);
215 STORM_LOG_ASSERT(model.getInitialStates().isSubsetOf(rewMinusInfEStates), "Initial state does not reach all maybestates");
216 }
217}
218
219template<typename ModelType>
223
224template<typename ModelType>
226 STORM_LOG_ASSERT(getInfinityCase() == InfinityCase::HasNegativeInfinite, "Tried to get -inf states, but there are none");
227 return rewMinusInfEStates;
228}
229
230template<typename ModelType>
232 return constantInitialStateValue.has_value();
233}
234
235template<typename ModelType>
237 assert(hasConstantInitialStateValue());
238 return *constantInitialStateValue;
239}
240
241template<typename ModelType>
242std::map<uint64_t, typename ModelType::ValueType> const& DeterministicSchedsObjectiveHelper<ModelType>::getChoiceRewards() const {
243 return choiceRewards;
244}
245
246template<typename ModelType>
250
251template<typename ModelType>
252typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper<ModelType>::getUpperValueBoundAtState(uint64_t state) const {
253 STORM_LOG_ASSERT(maybeStates.get(state), "Expected a maybestate.");
254 STORM_LOG_ASSERT(upperResultBounds.has_value(), "requested upper value bounds but they were not computed.");
255 return upperResultBounds->at(state);
256}
257
258template<typename ModelType>
259typename ModelType::ValueType const& DeterministicSchedsObjectiveHelper<ModelType>::getLowerValueBoundAtState(uint64_t state) const {
260 STORM_LOG_ASSERT(maybeStates.get(state), "Expected a maybestate.");
261 STORM_LOG_ASSERT(lowerResultBounds.has_value(), "requested lower value bounds but they were not computed.");
262 return lowerResultBounds->at(state);
263}
264
265template<typename ModelType>
269
270template<typename ModelType>
272 return objective.formula->isRewardOperatorFormula() && objective.formula->getSubformula().isTotalRewardFormula();
273}
274
275template<typename ModelType>
277 return objective.formula->hasBound();
278}
279
280template<typename ModelType>
282 STORM_LOG_ASSERT(hasThreshold(), "Trying to get a threshold but there is none");
283 STORM_LOG_THROW(!storm::logic::isStrict(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException,
284 "Objective " << *objective.originalFormula << ": Strict objective thresholds are not supported.");
285 ValueType threshold = objective.formula->template getThresholdAs<ValueType>();
286 if (storm::solver::minimize(objective.formula->getOptimalityType())) {
287 threshold *= -storm::utility::one<ValueType>(); // negate minimizing thresholds
288 STORM_LOG_THROW(!storm::logic::isLowerBound(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException,
289 "Objective " << *objective.originalFormula << ": Minimizing objective should specify an upper bound.");
290 } else {
291 STORM_LOG_THROW(storm::logic::isLowerBound(objective.formula->getBound().comparisonType), storm::exceptions::NotSupportedException,
292 "Objective " << *objective.originalFormula << ": Maximizing objective should specify a lower bound.");
293 }
294 return threshold;
295}
296
297template<typename ModelType>
299 return storm::solver::minimize(objective.formula->getOptimalityType());
300}
301
302template<typename ValueType>
304 std::vector<ValueType> const& rewards, std::vector<ValueType> const& exitProbabilities,
305 std::optional<storm::OptimizationDirection> dir, bool reqLower, bool reqUpper) {
306 if (!reqLower && !reqUpper) {
307 return; // nothing to be done!
308 }
309 STORM_LOG_ASSERT(!rewards.empty(), "empty reward vector,");
310
311 auto [minIt, maxIt] = std::minmax_element(rewards.begin(), rewards.end());
312 bool const hasNegativeValues = *minIt < storm::utility::zero<ValueType>();
313 bool const hasPositiveValues = *maxIt > storm::utility::zero<ValueType>();
314
315 std::optional<ValueType> lowerBound, upperBound;
316
317 // Get 0 as a trivial lower/upper bound if possible
318 if (!hasNegativeValues) {
319 lowerBound = storm::utility::zero<ValueType>();
320 }
321 if (!hasPositiveValues) {
322 upperBound = storm::utility::zero<ValueType>();
323 }
324
325 // Invoke the respective reward bound computers if needed
326 std::vector<ValueType> tmpRewards;
327 if (reqUpper && !upperBound.has_value()) {
328 if (hasNegativeValues) {
329 tmpRewards.resize(rewards.size());
330 storm::utility::vector::applyPointwise(rewards, tmpRewards, [](ValueType const& v) { return std::max(storm::utility::zero<ValueType>(), v); });
331 }
332 if (dir.has_value() && maximize(*dir)) {
333 solver.setUpperBound(
334 storm::modelchecker::helper::BaierUpperRewardBoundsComputer<ValueType>(matrix, hasNegativeValues ? tmpRewards : rewards, exitProbabilities)
335 .computeUpperBound());
336 } else {
337 solver.setUpperBounds(
338 storm::modelchecker::helper::DsMpiMdpUpperRewardBoundsComputer<ValueType>(matrix, hasNegativeValues ? tmpRewards : rewards, exitProbabilities)
339 .computeUpperBounds());
340 }
341 }
342 if (reqLower && !upperBound.has_value()) {
343 // For lower bounds we actually compute upper bounds for the negated rewards.
344 // We therefore need tmpRewards in any way.
345 tmpRewards.resize(rewards.size());
346 storm::utility::vector::applyPointwise(rewards, tmpRewards,
347 [](ValueType const& v) { return std::max<ValueType>(storm::utility::zero<ValueType>(), -v); });
348 if (dir.has_value() && minimize(*dir)) {
349 solver.setLowerBound(
351 } else {
352 auto lowerBounds =
354 storm::utility::vector::applyPointwise(lowerBounds, lowerBounds, [](ValueType const& v) { return -v; });
355 solver.setLowerBounds(std::move(lowerBounds));
356 }
357 }
358}
359
360template<typename ValueType>
361std::vector<ValueType> computeValuesOfReducedSystem(Environment const& env, storm::storage::SparseMatrix<ValueType> const& submatrix,
362 std::vector<ValueType> const& exitProbs, std::vector<ValueType> const& rewards,
363 storm::OptimizationDirection const& dir) {
365 minMaxSolver->setHasUniqueSolution(true);
366 minMaxSolver->setOptimizationDirection(dir);
367 auto req = minMaxSolver->getRequirements(env, dir);
368 setLowerUpperTotalRewardBoundsToSolver(*minMaxSolver, submatrix, rewards, exitProbs, dir, req.lowerBounds(), req.upperBounds());
369 req.clearBounds();
370 if (req.validInitialScheduler()) {
371 std::vector<uint64_t> initSched(submatrix.getRowGroupCount());
372 auto backwardsTransitions = submatrix.transpose(true);
373 storm::storage::BitVector statesWithChoice(initSched.size(), false);
374 storm::storage::BitVector foundStates(initSched.size(), false);
375 std::vector<uint64_t> stack;
376 for (uint64_t state = 0; state < initSched.size(); ++state) {
377 if (foundStates.get(state)) {
378 continue;
379 }
380 for (auto choice : submatrix.getRowGroupIndices(state)) {
381 if (!storm::utility::isZero(exitProbs[choice])) {
382 initSched[state] = choice - submatrix.getRowGroupIndices()[state];
383 statesWithChoice.set(state);
384 foundStates.set(state);
385 for (auto const& predecessor : backwardsTransitions.getRow(state)) {
386 if (!foundStates.get(predecessor.getColumn())) {
387 stack.push_back(predecessor.getColumn());
388 foundStates.set(predecessor.getColumn());
389 }
390 }
391 break;
392 }
393 }
394 }
395 while (!stack.empty()) {
396 auto state = stack.back();
397 stack.pop_back();
398 for (auto choice : submatrix.getRowGroupIndices(state)) {
399 auto row = submatrix.getRow(choice);
400 if (std::any_of(row.begin(), row.end(), [&statesWithChoice](auto const& entry) {
401 return !storm::utility::isZero(entry.getValue()) && statesWithChoice.get(entry.getColumn());
402 })) {
403 initSched[state] = choice - submatrix.getRowGroupIndices()[state];
404 statesWithChoice.set(state);
405 break;
406 }
407 }
408 assert(statesWithChoice.get(state));
409 for (auto const& predecessor : backwardsTransitions.getRow(state)) {
410 if (!foundStates.get(predecessor.getColumn())) {
411 stack.push_back(predecessor.getColumn());
412 foundStates.set(predecessor.getColumn());
413 }
414 }
415 }
416 assert(statesWithChoice.full());
417 minMaxSolver->setInitialScheduler(std::move(initSched));
418 req.clearValidInitialScheduler();
419 }
420 STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException,
421 "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked.");
422 minMaxSolver->setRequirementsChecked(true);
423
424 std::vector<ValueType> result(submatrix.getRowGroupCount());
425 minMaxSolver->solveEquations(env, result, rewards);
426
427 return result;
428}
429
430template<typename ValueType>
431void plusMinMaxSolverPrecision(Environment const& env, ValueType& value) {
433 auto eps = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision());
435 value += value * eps;
436 } else {
437 value += eps;
438 }
439 }
440}
441
442template<typename ValueType>
443void minusMinMaxSolverPrecision(Environment const& env, ValueType& value) {
445 auto eps = storm::utility::convertNumber<ValueType>(env.solver().minMax().getPrecision());
447 value -= value * eps;
448 } else {
449 value -= eps;
450 }
451 }
452}
453
454template<storm::OptimizationDirection Dir, typename ValueType>
455ValueType sumOfMecRewards(storm::storage::MaximalEndComponent const& mec, std::vector<ValueType> const& rewards) {
456 auto sum = storm::utility::zero<ValueType>();
457 for (auto const& stateChoices : mec) {
459 for (auto const& choice : stateChoices.second) {
460 optimalStateValue &= rewards.at(choice);
461 }
462 sum += *optimalStateValue;
463 }
464 return sum;
465}
466
467template<typename ValueType>
468ValueType getLowerBoundForNonZeroReachProb(storm::storage::SparseMatrix<ValueType> const& transitions, uint64_t init, uint64_t target) {
469 storm::storage::BitVector allStates(transitions.getRowGroupCount(), true);
470 storm::storage::BitVector initialAsBitVector(transitions.getRowGroupCount(), false);
471 initialAsBitVector.set(init, true);
472 storm::storage::BitVector targetAsBitVector(transitions.getRowGroupCount(), false);
473 targetAsBitVector.set(target, true);
474 auto relevantStates = storm::utility::graph::getReachableStates(transitions, initialAsBitVector, allStates, targetAsBitVector);
475 auto product = storm::utility::one<ValueType>();
476 for (auto state : relevantStates) {
477 if (state == target) {
478 continue;
479 }
481 for (auto const& entry : transitions.getRowGroup(state)) {
482 if (relevantStates.get(entry.getColumn())) {
483 minProb &= entry.getValue();
484 }
485 }
486 product *= *minProb;
487 }
488 return product;
489}
490
491template<typename ModelType>
493 assert(!upperResultBounds.has_value() && !lowerResultBounds.has_value());
494 auto backwardTransitions = model.getBackwardTransitions();
495 auto nonMaybeStates = ~maybeStates;
496 // Eliminate problematic mecs
497 storm::storage::MaximalEndComponentDecomposition<ValueType> problMecs(model.getTransitionMatrix(), backwardTransitions, maybeStates,
498 getRelevantZeroRewardChoices());
499 auto quotient1 = storm::transformer::EndComponentEliminator<ValueType>::transform(model.getTransitionMatrix(), problMecs, maybeStates, maybeStates);
500 auto backwardTransitions1 = quotient1.matrix.transpose(true);
501 std::vector<ValueType> rewards1(quotient1.matrix.getRowCount(), storm::utility::zero<ValueType>());
502 std::vector<ValueType> exitProbs1(quotient1.matrix.getRowCount(), storm::utility::zero<ValueType>());
503 storm::storage::BitVector subsystemChoices1(quotient1.matrix.getRowCount(), true);
504 storm::storage::BitVector allQuotient1States(quotient1.matrix.getRowGroupCount(), true);
505 for (uint64_t choice = 0; choice < quotient1.matrix.getRowCount(); ++choice) {
506 auto const& oldChoice = quotient1.newToOldRowMapping[choice];
507 if (auto findRes = choiceRewards.find(oldChoice); findRes != choiceRewards.end()) {
508 rewards1[choice] = findRes->second;
509 }
510 if (quotient1.sinkRows.get(choice)) {
511 subsystemChoices1.set(choice, false);
512 exitProbs1[choice] = storm::utility::one<ValueType>();
513 } else if (auto prob = model.getTransitionMatrix().getConstrainedRowSum(oldChoice, nonMaybeStates); !storm::utility::isZero(prob)) {
514 exitProbs1[choice] = prob;
515 subsystemChoices1.set(choice, false);
516 }
517 }
518 // Assert that every state can eventually exit the subsystem. If that is not the case, we get infinite reward as problematic (aka zero mecs) are removed.
519 auto exitStates1 = quotient1.matrix.getRowGroupFilter(~subsystemChoices1, false);
520 STORM_LOG_THROW(storm::utility::graph::performProbGreater0E(backwardTransitions1, allQuotient1States, exitStates1).full(),
521 storm::exceptions::InvalidOperationException, "Objective " << *objective.originalFormula << " does not induce finite value.");
522
523 // Compute bounds for finite cases
524 if (getInfinityCase() != InfinityCase::HasNegativeInfinite) {
525 auto result1 = computeValuesOfReducedSystem(env, quotient1.matrix, exitProbs1, rewards1, storm::OptimizationDirection::Minimize);
526 lowerResultBounds = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
527 for (auto const& state : maybeStates) {
528 ValueType val = result1.at(quotient1.oldToNewStateMapping.at(state));
530 (*lowerResultBounds)[state] = val;
531 }
532 }
533 if (getInfinityCase() != InfinityCase::HasPositiveInfinite) {
534 auto result1 = computeValuesOfReducedSystem(env, quotient1.matrix, exitProbs1, rewards1, storm::OptimizationDirection::Maximize);
535 upperResultBounds = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
536 for (auto const& state : maybeStates) {
537 ValueType val = result1.at(quotient1.oldToNewStateMapping.at(state));
539 if (infinityCase == InfinityCase::HasNegativeInfinite) { // Upper bound has to be at least 0 to allow for "trivial solution" in encoding
540 val = std::max(val, storm::utility::zero<ValueType>());
541 }
542 (*upperResultBounds)[state] = val;
543 }
544 }
545 if (getInfinityCase() != InfinityCase::AlwaysFinite) {
546 assert(getInfinityCase() == InfinityCase::HasPositiveInfinite || getInfinityCase() == InfinityCase::HasNegativeInfinite);
548 hasThreshold() || getInfinityCase() == InfinityCase::HasNegativeInfinite, storm::exceptions::NotSupportedException,
549 "The upper bound for objective " << *objective.originalFormula << " is infinity at some state. This is only supported for thresholded objectives");
550 // Eliminate remaining mecs
551 storm::storage::MaximalEndComponentDecomposition<ValueType> remainingMecs(quotient1.matrix, backwardTransitions1, allQuotient1States,
552 subsystemChoices1);
553 STORM_LOG_ASSERT(!remainingMecs.empty(), "Incorrect infinityCase: There is no MEC with rewards.");
554 auto quotient2 =
555 storm::transformer::EndComponentEliminator<ValueType>::transform(quotient1.matrix, remainingMecs, allQuotient1States, allQuotient1States);
556 std::vector<ValueType> rewards2(quotient2.matrix.getRowCount(), storm::utility::zero<ValueType>());
557 std::vector<ValueType> exitProbs2(quotient2.matrix.getRowCount(), storm::utility::zero<ValueType>());
558 for (uint64_t choice = 0; choice < quotient2.matrix.getRowCount(); ++choice) {
559 auto const& oldChoice = quotient2.newToOldRowMapping[choice];
560 if (quotient2.sinkRows.get(choice)) {
561 exitProbs2[choice] = storm::utility::one<ValueType>();
562 } else {
563 rewards2[choice] = rewards1[oldChoice];
564 exitProbs2[choice] = exitProbs1[oldChoice];
565 }
566 }
567
568 // Add rewards within ECs
569 auto initialState2 = quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(*model.getInitialStates().begin()));
570 ValueType rewardValueForPosInfCase;
571 if (getInfinityCase() == InfinityCase::HasPositiveInfinite) {
572 rewardValueForPosInfCase = getThreshold();
573 ;
574 // We need to substract a lower bound for the value at the initial state
575 std::vector<ValueType> rewards2Negative;
576 rewards2Negative.reserve(rewards2.size());
577 for (auto const& rew : rewards2) {
578 rewards2Negative.push_back(std::min(rew, storm::utility::zero<ValueType>()));
579 }
580 auto lower2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2Negative, storm::OptimizationDirection::Minimize);
581 rewardValueForPosInfCase -= lower2.at(initialState2);
582 }
583 for (auto const& mec : remainingMecs) {
584 auto mecState2 = quotient2.oldToNewStateMapping[mec.begin()->first];
585 ValueType mecReward = getInfinityCase() == InfinityCase::HasPositiveInfinite
586 ? sumOfMecRewards<storm::OptimizationDirection::Maximize>(mec, rewards1)
587 : sumOfMecRewards<storm::OptimizationDirection::Minimize>(mec, rewards1);
588 mecReward /= VisitingTimesHelper<ValueType>::computeMecTraversalLowerBound(mec, quotient1.matrix);
589 auto groupIndices = quotient2.matrix.getRowGroupIndices(mecState2);
590 for (auto const choice2 : groupIndices) {
591 rewards2[choice2] += mecReward;
592 }
593 if (getInfinityCase() == InfinityCase::HasPositiveInfinite) {
594 auto sinkChoice = quotient2.sinkRows.getNextSetIndex(*groupIndices.begin());
595 STORM_LOG_ASSERT(sinkChoice < quotient2.matrix.getRowGroupIndices()[mecState2 + 1], "EC state in quotient has no sink row.");
596 rewards2[sinkChoice] += rewardValueForPosInfCase / getLowerBoundForNonZeroReachProb(quotient2.matrix, initialState2, mecState2);
597 }
598 }
599
600 // Compute and insert missing bounds
601 if (getInfinityCase() == InfinityCase::HasNegativeInfinite) {
602 auto result2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2, storm::OptimizationDirection::Minimize);
603 lowerResultBounds = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
604 for (auto const& state : maybeStates) {
605 ValueType val = result2.at(quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(state)));
607 (*lowerResultBounds)[state] = val;
608 }
609 } else {
610 assert(getInfinityCase() == InfinityCase::HasPositiveInfinite);
611 auto result2 = computeValuesOfReducedSystem(env, quotient2.matrix, exitProbs2, rewards2, storm::OptimizationDirection::Maximize);
612 upperResultBounds = std::vector<ValueType>(model.getNumberOfStates(), storm::utility::zero<ValueType>());
613 for (auto const& state : maybeStates) {
614 ValueType val = result2.at(quotient2.oldToNewStateMapping.at(quotient1.oldToNewStateMapping.at(state)));
616 (*upperResultBounds)[state] = val;
617 }
618 }
619 }
621 std::all_of(maybeStates.begin(), maybeStates.end(), [this](auto const& s) { return this->lowerResultBounds->at(s) <= this->upperResultBounds->at(s); }),
622 storm::exceptions::UnexpectedException, "Pre-computed lower bound exceeds upper bound.");
623}
624
625template<typename ModelType>
627 Environment const& env, storm::storage::BitVector const& selectedChoices) const {
628 STORM_LOG_ASSERT(model.getInitialStates().getNumberOfSetBits() == 1u, "Expected a single initial state.");
629 STORM_LOG_ASSERT(model.getInitialStates().isSubsetOf(getMaybeStates()), "Expected initial state to be maybestate.");
630 STORM_LOG_ASSERT(selectedChoices.getNumberOfSetBits() == model.getNumberOfStates(), "invalid choice selection.");
631 storm::storage::BitVector allStates(model.getNumberOfStates(), true);
632 auto selectedMatrix = model.getTransitionMatrix().getSubmatrix(false, selectedChoices, allStates);
633 assert(selectedMatrix.getRowCount() == selectedMatrix.getRowGroupCount());
634 selectedMatrix.makeRowGroupingTrivial();
635 auto subMaybeStates =
636 getMaybeStates() & storm::utility::graph::getReachableStates(selectedMatrix, model.getInitialStates(), getMaybeStates(), ~getMaybeStates());
637 auto bsccCandidates = storm::utility::graph::performProbGreater0(selectedMatrix.transpose(true), subMaybeStates, ~subMaybeStates);
638 bsccCandidates.complement(); // i.e. states that can not reach non-maybe states
639 if (!bsccCandidates.empty()) {
641 bsccOptions.onlyBottomSccs(true).subsystem(bsccCandidates);
642 storm::storage::StronglyConnectedComponentDecomposition bsccs(selectedMatrix, bsccOptions);
643 for (auto const& bscc : bsccs) {
644 for (auto const& s : bscc) {
645 subMaybeStates.set(s, false);
646 }
647 }
648 }
649
650 if (subMaybeStates.get(*model.getInitialStates().begin())) {
653 auto eqSysMatrix = selectedMatrix.getSubmatrix(true, subMaybeStates, subMaybeStates, useEqSysFormat);
654 assert(eqSysMatrix.getRowCount() == eqSysMatrix.getRowGroupCount());
655 auto exitProbs = selectedMatrix.getConstrainedRowSumVector(subMaybeStates, ~subMaybeStates);
656 std::vector<ValueType> rewards;
657 rewards.reserve(exitProbs.size());
658 auto const& allRewards = getChoiceRewards();
659 for (auto const& state : getMaybeStates()) {
660 auto choice = selectedChoices.getNextSetIndex(model.getTransitionMatrix().getRowGroupIndices()[state]);
661 STORM_LOG_ASSERT(choice < model.getTransitionMatrix().getRowGroupIndices()[state + 1], " no choice selected at state" << state);
662 if (subMaybeStates.get(state)) {
663 if (auto findRes = allRewards.find(choice); findRes != allRewards.end()) {
664 rewards.push_back(findRes->second);
665 } else {
666 rewards.push_back(storm::utility::zero<ValueType>());
667 }
668 } else {
669 STORM_LOG_ASSERT(!bsccCandidates.get(state) || allRewards.count(choice) == 0, "Strategy selected a bscc with rewards");
670 }
671 }
672 if (useEqSysFormat) {
673 eqSysMatrix.convertToEquationSystem();
674 }
675 auto solver = factory.create(env, eqSysMatrix);
676 storm::storage::BitVector init = model.getInitialStates() % subMaybeStates;
677 assert(init.getNumberOfSetBits() == 1);
678 auto const initState = *init.begin();
679 solver->setRelevantValues(std::move(init));
680 auto req = solver->getRequirements(env);
681 if (!useEqSysFormat) {
682 setLowerUpperTotalRewardBoundsToSolver(*solver, eqSysMatrix, rewards, exitProbs, std::nullopt, req.lowerBounds(), req.upperBounds());
683 req.clearUpperBounds();
684 req.clearLowerBounds();
685 }
686 STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException,
687 "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked.");
688 std::vector<ValueType> x(rewards.size());
689 solver->solveEquations(env, x, rewards);
690
691 return x[initState];
692 } else {
693 // initial state is on a bscc
694 return storm::utility::zero<ValueType>();
695 }
696}
697
702} // namespace storm::modelchecker::multiobjective
SolverEnvironment & solver()
storm::RationalNumber const & getPrecision() const
bool const & getRelativeTerminationCriterion() const
MinMaxSolverEnvironment & minMax()
RewardOperatorFormula & asRewardOperatorFormula()
Definition Formula.cpp:461
virtual bool isRewardOperatorFormula() const
Definition Formula.cpp:176
virtual bool isTimeOperatorFormula() const
Definition Formula.cpp:140
boost::optional< std::string > const & getOptionalRewardModelName() const
Retrieves the optional representing the reward model name this property refers to.
virtual std::unique_ptr< CheckResult > check(Environment const &env, CheckTask< storm::logic::Formula, SolutionType > const &checkTask)
Checks the provided formula.
ValueType computeUpperBound()
Computes an upper bound on the expected rewards.
std::vector< ValueType > computeUpperBounds()
Computes upper bounds on the expected rewards.
DeterministicSchedsObjectiveHelper(ModelType const &model, Objective< ValueType > const &objective)
Provides helper functions for computing bounds on the expected visiting times.
This class represents a Markov automaton.
storm::storage::BitVector const & getMarkovianStates() const
Retrieves the set of Markovian states of the model.
ValueType const & getExitRate(storm::storage::sparse::state_type state) const
Retrieves the exit rate of the given state.
This class represents a (discrete-time) Markov decision process.
Definition Mdp.h:14
RewardModelType const & getUniqueRewardModel() const
Retrieves the unique reward model, if there exists exactly one.
Definition Model.cpp:303
storm::storage::SparseMatrix< ValueType > const & getTransitionMatrix() const
Retrieves the matrix representing the transitions of the model.
Definition Model.cpp:197
virtual uint_fast64_t getNumberOfStates() const override
Returns the number of states of the model.
Definition Model.cpp:162
RewardModelType const & getRewardModel(std::string const &rewardModelName) const
Retrieves the reward model with the given name, if one exists.
Definition Model.cpp:218
uint_fast64_t getNumberOfChoices(uint_fast64_t state) const
std::vector< ValueType > getTotalActionRewardVector(storm::storage::SparseMatrix< MatrixValueType > const &transitionMatrix, std::vector< MatrixValueType > const &stateRewardWeights) const
Creates a vector representing the complete action-based rewards, i.e., state-action- and transition-b...
std::vector< ValueType > getTotalRewardVector(storm::storage::SparseMatrix< MatrixValueType > const &transitionMatrix) const
Creates a vector representing the complete reward vector based on the state-, state-action- and trans...
void setUpperBound(ValueType const &value)
Sets an upper bound for the solution that can potentially be used by the solver.
void setLowerBound(ValueType const &value)
Sets a lower bound for the solution that can potentially be used by the solver.
void setUpperBounds(std::vector< ValueType > const &values)
Sets upper bounds for the solution that can potentially be used by the solver.
void setLowerBounds(std::vector< ValueType > const &values)
Sets lower bounds for the solution that can potentially be used by the solver.
virtual std::unique_ptr< LinearEquationSolver< ValueType > > create(Environment const &env) const override
Creates an equation solver with the current settings, but without a matrix.
virtual std::unique_ptr< MinMaxLinearEquationSolver< ValueType, SolutionType > > create(Environment const &env) const override
virtual LinearEquationSolverProblemFormat getEquationProblemFormat(Environment const &env) const
Retrieves the problem format that the solver expects if it was created with the current settings.
A bit vector that is internally represented as a vector of 64-bit values.
Definition BitVector.h:16
void complement()
Negates all bits in the bit vector.
uint64_t getNextSetIndex(uint64_t startingIndex) const
Retrieves the index of the bit that is the next bit set to true in the bit vector.
bool full() const
Retrieves whether all bits are set in this bit vector.
uint64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
void set(uint64_t index, bool value=true)
Sets the given truth value at the given index.
const_iterator begin() const
Returns an iterator to the indices of the set bits in the bit vector.
bool get(uint64_t index) const
Retrieves the truth value of the bit at the given index and performs a bound check.
bool empty() const
Checks if the decomposition is empty.
This class represents the decomposition of a nondeterministic model into its maximal end components.
This class represents a maximal end-component of a nondeterministic model.
A class that holds a possibly non-square matrix in the compressed row storage format.
const_rows getRow(index_type row) const
Returns an object representing the given row.
const_rows getRowGroup(index_type rowGroup) const
Returns an object representing the given row group.
index_type getRowGroupCount() const
Returns the number of row groups in the matrix.
std::vector< index_type > const & getRowGroupIndices() const
Returns the grouping of rows of this matrix.
storm::storage::SparseMatrix< value_type > transpose(bool joinGroups=false, bool keepZeros=false) const
Transposes the matrix.
This class represents the decomposition of a graph-like structure into its strongly connected compone...
static EndComponentEliminatorReturnType transform(storm::storage::SparseMatrix< ValueType > const &originalMatrix, storm::storage::MaximalEndComponentDecomposition< ValueType > ecs, storm::storage::BitVector const &subsystemStates, storm::storage::BitVector const &addSinkRowStates, bool addSelfLoopAtSinkStates=false)
Stores and manages an extremal (maximal or minimal) value.
Definition Extremum.h:15
#define STORM_LOG_ASSERT(cond, message)
Definition macros.h:11
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
bool isLowerBound(ComparisonType t)
bool isStrict(ComparisonType t)
std::vector< ValueType > getTotalRewardVector(storm::models::sparse::MarkovAutomaton< ValueType > const &model, storm::logic::Formula const &formula)
storm::storage::BitVector evaluatePropositionalFormula(ModelType const &model, storm::logic::Formula const &formula)
void minusMinMaxSolverPrecision(Environment const &env, ValueType &value)
ValueType sumOfMecRewards(storm::storage::MaximalEndComponent const &mec, std::vector< ValueType > const &rewards)
void plusMinMaxSolverPrecision(Environment const &env, ValueType &value)
ValueType getLowerBoundForNonZeroReachProb(storm::storage::SparseMatrix< ValueType > const &transitions, uint64_t init, uint64_t target)
std::vector< ValueType > computeValuesOfReducedSystem(Environment const &env, storm::storage::SparseMatrix< ValueType > const &submatrix, std::vector< ValueType > const &exitProbs, std::vector< ValueType > const &rewards, storm::OptimizationDirection const &dir)
void setLowerUpperTotalRewardBoundsToSolver(storm::solver::AbstractEquationSolver< ValueType > &solver, storm::storage::SparseMatrix< ValueType > const &matrix, std::vector< ValueType > const &rewards, std::vector< ValueType > const &exitProbabilities, std::optional< storm::OptimizationDirection > dir, bool reqLower, bool reqUpper)
bool constexpr minimize(OptimizationDirection d)
storm::storage::BitVector getReachableStates(storm::storage::SparseMatrix< T > const &transitionMatrix, storm::storage::BitVector const &initialStates, storm::storage::BitVector const &constraintStates, storm::storage::BitVector const &targetStates, bool useStepBound, uint_fast64_t maximalSteps, boost::optional< storm::storage::BitVector > const &choiceFilter)
Performs a forward depth-first search through the underlying graph structure to identify the states t...
Definition graph.cpp:47
storm::storage::BitVector performProbGreater0(storm::storage::SparseMatrix< T > const &backwardTransitions, storm::storage::BitVector const &phiStates, storm::storage::BitVector const &psiStates, bool useStepBound, uint_fast64_t maximalSteps)
Performs a backward depth-first search trough the underlying graph structure of the given model to de...
Definition graph.cpp:321
storm::storage::BitVector performProb1A(storm::models::sparse::NondeterministicModel< T, RM > const &model, storm::storage::SparseMatrix< T > const &backwardTransitions, storm::storage::BitVector const &phiStates, storm::storage::BitVector const &psiStates)
Computes the sets of states that have probability 1 of satisfying phi until psi under all possible re...
Definition graph.cpp:987
storm::storage::BitVector performProbGreater0E(storm::storage::SparseMatrix< T > const &backwardTransitions, storm::storage::BitVector const &phiStates, storm::storage::BitVector const &psiStates, bool useStepBound, uint_fast64_t maximalSteps)
Computes the sets of states that have probability greater 0 of satisfying phi until psi under at leas...
Definition graph.cpp:679
storm::storage::BitVector performProb0A(storm::storage::SparseMatrix< T > const &backwardTransitions, storm::storage::BitVector const &phiStates, storm::storage::BitVector const &psiStates)
Definition graph.cpp:739
bool checkIfECWithChoiceExists(storm::storage::SparseMatrix< T > const &transitionMatrix, storm::storage::SparseMatrix< T > const &backwardTransitions, storm::storage::BitVector const &subsystem, storm::storage::BitVector const &choices)
Checks whether there is an End Component that.
Definition graph.cpp:181
void applyPointwise(std::vector< InValueType1 > const &firstOperand, std::vector< InValueType2 > const &secondOperand, std::vector< OutValueType > &target, Operation f=Operation())
Applies the given operation pointwise on the two given vectors and writes the result to the third vec...
Definition vector.h:378
FilteredRewardModel< RewardModelType > createFilteredRewardModel(RewardModelType const &baseRewardModel, storm::logic::RewardAccumulation const &acc, bool isDiscreteTimeModel)
bool isZero(ValueType const &a)
Definition constants.cpp:41
StronglyConnectedComponentDecompositionOptions & subsystem(storm::storage::BitVector const &subsystem)
Sets a bit vector indicating which subsystem to consider for the decomposition into SCCs.
StronglyConnectedComponentDecompositionOptions & onlyBottomSccs(bool value=true)
Sets if only bottom SCCs, i.e. SCCs in which all states have no way of leaving the SCC),...