Storm 1.10.0.1
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
SparsePcaaQuantitativeQuery.cpp
Go to the documentation of this file.
2
12
14
15namespace storm {
16namespace modelchecker {
17namespace multiobjective {
18
19template<class SparseModelType, typename GeometryValueType>
22 : SparsePcaaQuery<SparseModelType, GeometryValueType>(preprocessorResult) {
24 "Invalid query Type");
25
26 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
27 if (!this->objectives[objIndex].formula->hasBound()) {
28 indexOfOptimizingObjective = objIndex;
29 break;
30 }
31 }
32 initializeThresholdData();
33
34 // Set the maximum distance between lower and upper bound of the weightVectorChecker result.
35 this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(0.1));
36}
37
38template<class SparseModelType, typename GeometryValueType>
40 thresholds.reserve(this->objectives.size());
41 strictThresholds = storm::storage::BitVector(this->objectives.size(), false);
42 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> thresholdConstraints;
43 thresholdConstraints.reserve(this->objectives.size() - 1);
44 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
45 auto const& formula = *this->objectives[objIndex].formula;
46
47 if (formula.hasBound()) {
48 thresholds.push_back(formula.template getThresholdAs<GeometryValueType>());
49 if (storm::solver::minimize(formula.getOptimalityType())) {
50 STORM_LOG_ASSERT(!storm::logic::isLowerBound(formula.getBound().comparisonType), "Minimizing objective should not specify an upper bound.");
51 // Values for minimizing objectives will be negated in order to convert them to maximizing objectives.
52 // Hence, we also negate the threshold
53 thresholds.back() *= -storm::utility::one<GeometryValueType>();
54 }
55 strictThresholds.set(objIndex, storm::logic::isStrict(formula.getBound().comparisonType));
56 WeightVector normalVector(this->objectives.size(), storm::utility::zero<GeometryValueType>());
57 normalVector[objIndex] = -storm::utility::one<GeometryValueType>();
58 thresholdConstraints.emplace_back(std::move(normalVector), -thresholds.back());
59 } else {
60 thresholds.push_back(storm::utility::zero<GeometryValueType>());
61 }
62 }
63 // Note: If we have a single objective (i.e., no objectives with thresholds), thresholdsAsPolytope gets no constraints
64 thresholdsAsPolytope = storm::storage::geometry::Polytope<GeometryValueType>::create(thresholdConstraints);
65}
66
67template<class SparseModelType, typename GeometryValueType>
68std::unique_ptr<CheckResult> SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::check(Environment const& env, bool produceScheduler) {
69 STORM_LOG_THROW(!produceScheduler, storm::exceptions::NotImplementedException, "Scheduler computation is not implement for achievability queries.");
70
71 // First find one solution that achieves the given thresholds ...
72 if (this->checkAchievability(env)) {
73 // ... then improve it
74 GeometryValueType result = this->improveSolution(env);
75
76 // transform the obtained result for the preprocessed model to a result w.r.t. the original model and return the checkresult
77 auto const& obj = this->objectives[indexOfOptimizingObjective];
78 if (storm::solver::maximize(obj.formula->getOptimalityType())) {
79 if (obj.considersComplementaryEvent) {
80 result = storm::utility::one<GeometryValueType>() - result;
81 }
82 } else {
83 if (obj.considersComplementaryEvent) {
84 result = storm::utility::one<GeometryValueType>() + result; // 1 - (-result)
85 } else {
86 result *= -storm::utility::one<GeometryValueType>();
87 }
88 }
89 auto resultForOriginalModel = storm::utility::convertNumber<typename SparseModelType::ValueType>(result);
90
92 this->originalModel.getInitialStates().getNextSetIndex(0), resultForOriginalModel));
93 } else {
94 return std::unique_ptr<CheckResult>(new ExplicitQualitativeCheckResult(this->originalModel.getInitialStates().getNextSetIndex(0), false));
95 }
96}
97
98template<class SparseModelType, typename GeometryValueType>
100 if (this->objectives.size() > 1) {
101 // We don't care for the optimizing objective at this point
102 this->diracWeightVectorsToBeChecked.set(indexOfOptimizingObjective, false);
103
104 while (!this->maxStepsPerformed(env) && !storm::utility::resources::isTerminate()) {
105 WeightVector separatingVector = this->findSeparatingVector(thresholds);
106 this->updateWeightedPrecisionInAchievabilityPhase(separatingVector);
107 this->performRefinementStep(env, std::move(separatingVector), false); // scheduler computation currently not supported
108 // Pick the threshold for the optimizing objective low enough so valid solutions are not excluded
109 thresholds[indexOfOptimizingObjective] =
110 std::min(thresholds[indexOfOptimizingObjective], this->refinementSteps.back().lowerBoundPoint[indexOfOptimizingObjective]);
111 if (!checkIfThresholdsAreSatisfied(this->overApproximation)) {
112 return false;
113 }
114 if (checkIfThresholdsAreSatisfied(this->underApproximation)) {
115 return true;
116 }
117 }
118 } else {
119 // If there is only one objective than its the optimizing one. Thus the query has to be achievable.
120 return true;
121 }
122 STORM_LOG_ERROR("Could not check whether thresholds are achievable: Termination requested or maximum number of refinement steps exceeded.");
123 return false;
124}
125
126template<class SparseModelType, typename GeometryValueType>
127void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInAchievabilityPhase(WeightVector const& weights) {
128 // Our heuristic considers the distance between the under- and the over approximation w.r.t. the given direction
129 std::pair<Point, bool> optimizationResOverApprox = this->overApproximation->optimize(weights);
130 if (optimizationResOverApprox.second) {
131 std::pair<Point, bool> optimizationResUnderApprox = this->underApproximation->optimize(weights);
132 if (optimizationResUnderApprox.second) {
133 GeometryValueType distance = storm::utility::vector::dotProduct(optimizationResOverApprox.first, weights) -
134 storm::utility::vector::dotProduct(optimizationResUnderApprox.first, weights);
135 STORM_LOG_ASSERT(distance >= storm::utility::zero<GeometryValueType>(), "Negative distance between under- and over approximation was not expected");
136 // Normalize the distance by dividing it with the Euclidean Norm of the weight-vector
137 distance /= storm::utility::sqrt(storm::utility::vector::dotProduct(weights, weights));
138 distance /= GeometryValueType(2);
139 this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(distance));
140 }
141 }
142 // do not update the precision if one of the approximations is unbounded in the provided direction
143}
144
145template<class SparseModelType, typename GeometryValueType>
146GeometryValueType SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::improveSolution(Environment const& env) {
147 STORM_LOG_THROW(env.modelchecker().multi().getPrecisionType() == MultiObjectiveModelCheckerEnvironment::PrecisionType::Absolute,
148 storm::exceptions::IllegalArgumentException, "Unhandled multiobjective precision type.");
149 this->diracWeightVectorsToBeChecked.clear(); // Only check weight vectors that can actually improve the solution
150
151 WeightVector directionOfOptimizingObjective(this->objectives.size(), storm::utility::zero<GeometryValueType>());
152 directionOfOptimizingObjective[indexOfOptimizingObjective] = storm::utility::one<GeometryValueType>();
153
154 // Improve the found solution.
155 // Note that we do not care whether a threshold is strict anymore, because the resulting optimum should be
156 // the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict
157 // thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds.
158 GeometryValueType result = storm::utility::zero<GeometryValueType>();
159 while (!this->maxStepsPerformed(env) && !storm::utility::resources::isTerminate()) {
160 if (this->refinementSteps.empty()) {
161 // We did not make any refinement steps during the checkAchievability phase (e.g., because there is only one objective).
162 this->weightVectorChecker->setWeightedPrecision(
163 storm::utility::convertNumber<typename SparseModelType::ValueType>(env.modelchecker().multi().getPrecision()));
164 WeightVector separatingVector = directionOfOptimizingObjective;
165 this->performRefinementStep(env, std::move(separatingVector), false); // scheduler computation currently not supported
166 }
167 std::pair<Point, bool> optimizationRes = this->underApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective);
168 STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty.");
169 result = optimizationRes.first[indexOfOptimizingObjective];
170 STORM_LOG_DEBUG("Best solution found so far is ~" << storm::utility::convertNumber<double>(result) << ".");
171 // Compute an upper bound for the optimum and check for convergence
172 optimizationRes = this->overApproximation->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective);
173 if (optimizationRes.second) {
174 GeometryValueType precisionOfResult = optimizationRes.first[indexOfOptimizingObjective] - result;
175 if (precisionOfResult < storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision())) {
176 // Goal precision reached!
177 return result;
178 } else {
179 STORM_LOG_DEBUG("Solution can be improved by at most " << storm::utility::convertNumber<double>(precisionOfResult));
180 thresholds[indexOfOptimizingObjective] = optimizationRes.first[indexOfOptimizingObjective];
181 }
182 } else {
183 thresholds[indexOfOptimizingObjective] = result + storm::utility::one<GeometryValueType>();
184 }
185 WeightVector separatingVector = this->findSeparatingVector(thresholds);
186 this->updateWeightedPrecisionInImprovingPhase(env, separatingVector);
187 this->performRefinementStep(env, std::move(separatingVector), false); // scheduler computation currently not supported
188 }
189 STORM_LOG_ERROR("Could not reach the desired precision: Termination requested or maximum number of refinement steps exceeded.");
190 return result;
191}
192
193template<class SparseModelType, typename GeometryValueType>
194void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInImprovingPhase(Environment const& env,
195 WeightVector const& weights) {
196 STORM_LOG_THROW(!storm::utility::isZero(weights[this->indexOfOptimizingObjective]), exceptions::UnexpectedException,
197 "The chosen weight-vector gives zero weight for the objective that is to be optimized.");
198 STORM_LOG_THROW(env.modelchecker().multi().getPrecisionType() == MultiObjectiveModelCheckerEnvironment::PrecisionType::Absolute,
199 storm::exceptions::IllegalArgumentException, "Unhandled multiobjective precision type.");
200 // If weighs[indexOfOptimizingObjective] is low, the computation of the weightVectorChecker needs to be more precise.
201 // Our heuristic ensures that if p is the new vertex of the under-approximation, then max{ eps | p' = p + (0..0 eps 0..0) is in the over-approximation } <=
202 // multiobjective_precision/0.9
203 GeometryValueType weightedPrecision =
204 weights[this->indexOfOptimizingObjective] * storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
205 // Normalize by division with the Euclidean Norm of the weight-vector
206 weightedPrecision /= storm::utility::sqrt(storm::utility::vector::dotProduct(weights, weights));
207 weightedPrecision *= storm::utility::convertNumber<GeometryValueType>(0.9);
208
209 this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(weightedPrecision));
210}
211
212template<class SparseModelType, typename GeometryValueType>
213bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::checkIfThresholdsAreSatisfied(
214 std::shared_ptr<storm::storage::geometry::Polytope<GeometryValueType>> const& polytope) {
215 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = polytope->getHalfspaces();
216 for (auto const& h : halfspaces) {
217 if (storm::utility::isZero(h.distance(thresholds))) {
218 // Check if the threshold point is on the boundary of the halfspace and whether this is violates strict thresholds
219 if (h.isPointOnBoundary(thresholds)) {
220 for (auto strictThreshold : strictThresholds) {
221 if (h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
222 return false;
223 }
224 }
225 }
226 } else {
227 return false;
228 }
229 }
230 return true;
231}
232
233#ifdef STORM_HAVE_CARL
234template class SparsePcaaQuantitativeQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>;
235template class SparsePcaaQuantitativeQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
236
237template class SparsePcaaQuantitativeQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
238template class SparsePcaaQuantitativeQuery<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
239#endif
240} // namespace multiobjective
241} // namespace modelchecker
242} // namespace storm
SparsePcaaQuantitativeQuery(preprocessing::SparseMultiObjectivePreprocessorResult< SparseModelType > &preprocessorResult)
virtual std::unique_ptr< CheckResult > check(Environment const &env, bool produceScheduler) override
std::unique_ptr< PcaaWeightVectorChecker< SparseModelType > > weightVectorChecker
std::vector< Objective< typename SparseModelType::ValueType > > objectives
A bit vector that is internally represented as a vector of 64-bit values.
Definition BitVector.h:16
static std::shared_ptr< Polytope< ValueType > > create(std::vector< Halfspace< ValueType > > const &halfspaces)
Creates a polytope from the given halfspaces.
#define STORM_LOG_DEBUG(message)
Definition logging.h:18
#define STORM_LOG_ERROR(message)
Definition logging.h:26
#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)
bool constexpr maximize(OptimizationDirection d)
bool constexpr minimize(OptimizationDirection d)
bool isTerminate()
Check whether the program should terminate (due to some abort signal).
T dotProduct(std::vector< T > const &firstOperand, std::vector< T > const &secondOperand)
Computes the dot product (aka scalar product) and returns the result.
Definition vector.h:477
bool isZero(ValueType const &a)
Definition constants.cpp:41
ValueType sqrt(ValueType const &number)
LabParser.cpp.