16namespace modelchecker {
17namespace multiobjective {
19template<
class SparseModelType,
typename GeometryValueType>
22 :
SparsePcaaQuery<SparseModelType, GeometryValueType>(preprocessorResult) {
24 "Invalid query Type");
26 for (uint_fast64_t objIndex = 0; objIndex < this->
objectives.size(); ++objIndex) {
27 if (!this->
objectives[objIndex].formula->hasBound()) {
28 indexOfOptimizingObjective = objIndex;
32 initializeThresholdData();
35 this->
weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(0.1));
38template<
class SparseModelType,
typename GeometryValueType>
40 thresholds.reserve(this->objectives.size());
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;
47 if (formula.hasBound()) {
48 thresholds.push_back(formula.template getThresholdAs<GeometryValueType>());
53 thresholds.back() *= -storm::utility::one<GeometryValueType>();
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());
60 thresholds.push_back(storm::utility::zero<GeometryValueType>());
67template<
class SparseModelType,
typename GeometryValueType>
69 STORM_LOG_THROW(!produceScheduler, storm::exceptions::NotImplementedException,
"Scheduler computation is not implement for achievability queries.");
72 if (this->checkAchievability(env)) {
74 GeometryValueType result = this->improveSolution(env);
77 auto const& obj = this->objectives[indexOfOptimizingObjective];
79 if (obj.considersComplementaryEvent) {
80 result = storm::utility::one<GeometryValueType>() - result;
83 if (obj.considersComplementaryEvent) {
84 result = storm::utility::one<GeometryValueType>() + result;
86 result *= -storm::utility::one<GeometryValueType>();
89 auto resultForOriginalModel = storm::utility::convertNumber<typename SparseModelType::ValueType>(result);
92 this->originalModel.getInitialStates().getNextSetIndex(0), resultForOriginalModel));
98template<
class SparseModelType,
typename GeometryValueType>
100 if (this->objectives.size() > 1) {
102 this->diracWeightVectorsToBeChecked.set(indexOfOptimizingObjective,
false);
105 WeightVector separatingVector = this->findSeparatingVector(thresholds);
106 this->updateWeightedPrecisionInAchievabilityPhase(separatingVector);
107 this->performRefinementStep(env, std::move(separatingVector),
false);
109 thresholds[indexOfOptimizingObjective] =
110 std::min(thresholds[indexOfOptimizingObjective], this->refinementSteps.back().lowerBoundPoint[indexOfOptimizingObjective]);
111 if (!checkIfThresholdsAreSatisfied(this->overApproximation)) {
114 if (checkIfThresholdsAreSatisfied(this->underApproximation)) {
122 STORM_LOG_ERROR(
"Could not check whether thresholds are achievable: Termination requested or maximum number of refinement steps exceeded.");
126template<
class SparseModelType,
typename GeometryValueType>
127void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInAchievabilityPhase(WeightVector
const& weights) {
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) {
135 STORM_LOG_ASSERT(distance >= storm::utility::zero<GeometryValueType>(),
"Negative distance between under- and over approximation was not expected");
138 distance /= GeometryValueType(2);
139 this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(distance));
145template<
class SparseModelType,
typename GeometryValueType>
146GeometryValueType SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::improveSolution(Environment
const& env) {
148 storm::exceptions::IllegalArgumentException,
"Unhandled multiobjective precision type.");
149 this->diracWeightVectorsToBeChecked.clear();
151 WeightVector directionOfOptimizingObjective(this->objectives.size(), storm::utility::zero<GeometryValueType>());
152 directionOfOptimizingObjective[indexOfOptimizingObjective] = storm::utility::one<GeometryValueType>();
158 GeometryValueType result = storm::utility::zero<GeometryValueType>();
160 if (this->refinementSteps.empty()) {
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);
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) <<
".");
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())) {
179 STORM_LOG_DEBUG(
"Solution can be improved by at most " << storm::utility::convertNumber<double>(precisionOfResult));
180 thresholds[indexOfOptimizingObjective] = optimizationRes.first[indexOfOptimizingObjective];
183 thresholds[indexOfOptimizingObjective] = result + storm::utility::one<GeometryValueType>();
185 WeightVector separatingVector = this->findSeparatingVector(thresholds);
186 this->updateWeightedPrecisionInImprovingPhase(env, separatingVector);
187 this->performRefinementStep(env, std::move(separatingVector),
false);
189 STORM_LOG_ERROR(
"Could not reach the desired precision: Termination requested or maximum number of refinement steps exceeded.");
193template<
class SparseModelType,
typename GeometryValueType>
194void SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::updateWeightedPrecisionInImprovingPhase(Environment
const& env,
195 WeightVector
const& weights) {
197 "The chosen weight-vector gives zero weight for the objective that is to be optimized.");
199 storm::exceptions::IllegalArgumentException,
"Unhandled multiobjective precision type.");
203 GeometryValueType weightedPrecision =
204 weights[this->indexOfOptimizingObjective] * storm::utility::convertNumber<GeometryValueType>(env.modelchecker().multi().getPrecision());
207 weightedPrecision *= storm::utility::convertNumber<GeometryValueType>(0.9);
209 this->weightVectorChecker->setWeightedPrecision(storm::utility::convertNumber<typename SparseModelType::ValueType>(weightedPrecision));
212template<
class SparseModelType,
typename GeometryValueType>
213bool SparsePcaaQuantitativeQuery<SparseModelType, GeometryValueType>::checkIfThresholdsAreSatisfied(
215 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = polytope->getHalfspaces();
216 for (
auto const& h : halfspaces) {
219 if (h.isPointOnBoundary(thresholds)) {
220 for (
auto strictThreshold : strictThresholds) {
221 if (h.normalVector()[strictThreshold] > storm::utility::zero<GeometryValueType>()) {
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>;
237template class SparsePcaaQuantitativeQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
238template class SparsePcaaQuantitativeQuery<storm::models::sparse::MarkovAutomaton<storm::RationalNumber>, storm::RationalNumber>;
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.
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)
#define STORM_LOG_ERROR(message)
#define STORM_LOG_ASSERT(cond, message)
#define STORM_LOG_THROW(cond, exception, message)
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.
bool isZero(ValueType const &a)
ValueType sqrt(ValueType const &number)