diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp
index 0ecbdfb24..d6d8cf96f 100644
--- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp
+++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveModelCheckerHelper.cpp
@@ -93,7 +93,46 @@ namespace storm {
return h;
}
}
- return underApproximation->findSeparatingHalfspace(pointToBeSeparated);
+
+ // We now use LP solving to find a seperating halfspace.
+ // Note that StoRM's LPSover does not support rational numbers. Hence, we optimize a polytope instead.
+ // The variables of the LP are the normalVector entries w_0 ... w_{n-1} and the minimal distance d between the halfspace and the paretoOptimalPoints
+ uint_fast64_t numVariables = pointToBeSeparated.size() + 1;
+ std::vector<storm::storage::geometry::Halfspace<RationalNumberType>> constraints;
+ constraints.reserve((numVariables-1) + 2 + paretoOptimalPoints.size());
+ // w_i >= 0 and d>= 0 <==> -w_i <= 0 and -d <= 0
+ for(uint_fast64_t i = 0; i<numVariables; ++i) {
+ WeightVector constraintCoefficients(numVariables, storm::utility::zero<RationalNumberType>());
+ constraintCoefficients[i] = -storm::utility::one<RationalNumberType>();
+ RationalNumberType constraintOffset = storm::utility::zero<RationalNumberType>();
+ constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+ }
+ // sum_i w_i == 1 <==> sum_i w_i <=1 and sum_i -w_i <= -1
+ {
+ WeightVector constraintCoefficients(numVariables, storm::utility::one<RationalNumberType>());
+ constraintCoefficients.back() = storm::utility::zero<RationalNumberType>();
+ RationalNumberType constraintOffset = storm::utility::one<RationalNumberType>();
+ constraints.emplace_back(constraintCoefficients, constraintOffset);
+ storm::utility::vector::scaleVectorInPlace(constraintCoefficients, -storm::utility::one<RationalNumberType>());
+ constraintOffset = -storm::utility::one<RationalNumberType>();
+ constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+ }
+ // let q=pointToBeSeparated. For each paretoPoint x: q*w - x*w >= d <==> (x-q) * w + d <= 0
+ for(auto const& paretoPoint : paretoOptimalPoints){
+ WeightVector constraintCoefficients(numVariables-1);
+ storm::utility::vector::subtractVectors(paretoPoint.first, pointToBeSeparated, constraintCoefficients);
+ constraintCoefficients.push_back(storm::utility::one<RationalNumberType>());
+ RationalNumberType constraintOffset = storm::utility::zero<RationalNumberType>();
+ constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset));
+ }
+
+ WeightVector optimizationVector(numVariables);
+ optimizationVector.back() = storm::utility::one<RationalNumberType>(); // maximize d
+ std::pair<WeightVector, bool> optimizeResult = storm::storage::geometry::Polytope<RationalNumberType>::create(constraints)->optimize(optimizationVector);
+ STORM_LOG_THROW(optimizeResult.second, storm::exceptions::UnexpectedException, "There seems to be no seperating halfspace.");
+ optimizeResult.first.pop_back(); //drop the distance
+ RationalNumberType offset = storm::utility::vector::dotProduct(optimizeResult.first, pointToBeSeparated);
+ return storm::storage::geometry::Halfspace<RationalNumberType>(std::move(optimizeResult.first), std::move(offset));
}
template <class SparseModelType, typename RationalNumberType>
diff --git a/src/storage/geometry/HyproPolytope.cpp b/src/storage/geometry/HyproPolytope.cpp
index 0def7ae64..4efb965c9 100644
--- a/src/storage/geometry/HyproPolytope.cpp
+++ b/src/storage/geometry/HyproPolytope.cpp
@@ -148,7 +148,7 @@ namespace storm {
std::shared_ptr<Polytope<ValueType>> HyproPolytope<ValueType>::downwardClosure(boost::optional<Point> const& upperBounds) const {
if(this->isUniversal() || this->isEmpty()) {
// In these cases, the polytope does not change, i.e., we return a copy of this
- return std::make_shared<HyproPolytope<ValueType>>(internPolytope);
+ return std::make_shared<HyproPolytope<ValueType>>(*this);
}
// Only keep the halfspaces where all entries of the normal vector are non-negative
std::vector<hypro::Halfspace<ValueType>> halfspaces;
@@ -184,6 +184,22 @@ namespace storm {
return std::make_shared<HyproPolytope<ValueType>>(HyproPolytopeType(std::move(halfspaces)));
}
+
+ template <typename ValueType>
+ std::pair<typename HyproPolytope<ValueType>::Point, bool> HyproPolytope<ValueType>::optimize(Point const& direction) const {
+ hypro::EvaluationResult<ValueType> evalRes = internPolytope.evaluate(storm::adapters::toHypro(direction));
+ switch (evalRes.errorCode) {
+ case hypro::SOLUTION::FEAS:
+ return std::make_pair(storm::adapters::fromHypro(evalRes.optimumValue), true);
+ case hypro::SOLUTION::UNKNOWN:
+ STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unexpected eror code for Polytope evaluation");
+ return std::make_pair(Point(), false);
+ default:
+ //solution is infinity or infeasible
+ return std::make_pair(Point(), false);
+ }
+ }
+
template <typename ValueType>
bool HyproPolytope<ValueType>::isHyproPolytope() const {
return true;
diff --git a/src/storage/geometry/HyproPolytope.h b/src/storage/geometry/HyproPolytope.h
index c17dfd775..2f523c769 100644
--- a/src/storage/geometry/HyproPolytope.h
+++ b/src/storage/geometry/HyproPolytope.h
@@ -107,6 +107,15 @@ namespace storm {
* @param upperBounds If given, this vector is considered for y (hence, max{x_i | x i \in P does not need to be computed)
*/
virtual std::shared_ptr<Polytope<ValueType>> downwardClosure(boost::optional<Point> const& upperBounds = boost::none) const override;
+
+ /*!
+ * Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,
+ * a point that maximizes dotPorduct(point, direction).
+ * If such a point does not exist, the returned bool is false. There are two reasons for this:
+ * - The polytope is empty
+ * - The polytope is not bounded in the given direction
+ */
+ virtual std::pair<Point, bool> optimize(Point const& direction) const override;
virtual bool isHyproPolytope() const override;
diff --git a/src/storage/geometry/Polytope.cpp b/src/storage/geometry/Polytope.cpp
index 0b7c75d69..f2fefa2bd 100644
--- a/src/storage/geometry/Polytope.cpp
+++ b/src/storage/geometry/Polytope.cpp
@@ -129,9 +129,9 @@ namespace storm {
}
template <typename ValueType>
- Halfspace<ValueType> Polytope<ValueType>::findSeparatingHalfspace(Point const& pointToBeSeparated) const {
+ std::pair<typename Polytope<ValueType>::Point, bool> Polytope<ValueType>::optimize(Point const& direction) const {
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented.");
- return Halfspace<ValueType>(Point(), storm::utility::zero<ValueType>());
+ return std::make_pair(Point(), false);
}
template <typename ValueType>
diff --git a/src/storage/geometry/Polytope.h b/src/storage/geometry/Polytope.h
index e95c894c2..8226a70e2 100644
--- a/src/storage/geometry/Polytope.h
+++ b/src/storage/geometry/Polytope.h
@@ -99,15 +99,13 @@ namespace storm {
virtual std::shared_ptr<Polytope<ValueType>> downwardClosure(boost::optional<Point> const& upperBounds = boost::none) const;
/*!
- * Returns a halfspace h that separates this polytope from the given point p, i.e.,
- * - p lies on the border of h and
- * - for each x in this polytope, it holds that h contains x
- *
- * A halfspace with maximal distance to the polytope is preferred.
- *
- * @param pointToBeSeparated the point that is to be seperated
+ * Finds an optimal point inside this polytope w.r.t. the given direction, i.e.,
+ * a point that maximizes dotPorduct(point, direction).
+ * If such a point does not exist, the returned bool is false. There are two reasons for this:
+ * - The polytope is empty
+ * - The polytope is not bounded in the given direction
*/
- virtual Halfspace<ValueType> findSeparatingHalfspace(Point const& pointToBeSeparated) const;
+ virtual std::pair<Point, bool> optimize(Point const& direction) const;
/*
* Returns a string representation of this polytope.