From fb1fa2f23c63ee8bf6d6c5700d5fbb67664723c8 Mon Sep 17 00:00:00 2001 From: TimQu Date: Sun, 5 Jun 2016 12:24:53 +0200 Subject: [PATCH] implemented the LP solving to find a separating halfspace Former-commit-id: d88558db0b027e81498e22cbce8862a9d06c7ae8 --- ...SparseMultiObjectiveModelCheckerHelper.cpp | 41 ++++++++++++++++++- src/storage/geometry/HyproPolytope.cpp | 18 +++++++- src/storage/geometry/HyproPolytope.h | 9 ++++ src/storage/geometry/Polytope.cpp | 4 +- src/storage/geometry/Polytope.h | 14 +++---- 5 files changed, 74 insertions(+), 12 deletions(-) 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> 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()); + constraintCoefficients[i] = -storm::utility::one(); + RationalNumberType constraintOffset = storm::utility::zero(); + 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()); + constraintCoefficients.back() = storm::utility::zero(); + RationalNumberType constraintOffset = storm::utility::one(); + constraints.emplace_back(constraintCoefficients, constraintOffset); + storm::utility::vector::scaleVectorInPlace(constraintCoefficients, -storm::utility::one()); + constraintOffset = -storm::utility::one(); + 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 constraintOffset = storm::utility::zero(); + constraints.emplace_back(std::move(constraintCoefficients), std::move(constraintOffset)); + } + + WeightVector optimizationVector(numVariables); + optimizationVector.back() = storm::utility::one(); // maximize d + std::pair optimizeResult = storm::storage::geometry::Polytope::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(std::move(optimizeResult.first), std::move(offset)); } template 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> HyproPolytope::downwardClosure(boost::optional 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>(internPolytope); + return std::make_shared>(*this); } // Only keep the halfspaces where all entries of the normal vector are non-negative std::vector> halfspaces; @@ -184,6 +184,22 @@ namespace storm { return std::make_shared>(HyproPolytopeType(std::move(halfspaces))); } + + template + std::pair::Point, bool> HyproPolytope::optimize(Point const& direction) const { + hypro::EvaluationResult 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 bool HyproPolytope::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> downwardClosure(boost::optional 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 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 - Halfspace Polytope::findSeparatingHalfspace(Point const& pointToBeSeparated) const { + std::pair::Point, bool> Polytope::optimize(Point const& direction) const { STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Functionality not implemented."); - return Halfspace(Point(), storm::utility::zero()); + return std::make_pair(Point(), false); } template 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> downwardClosure(boost::optional 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 findSeparatingHalfspace(Point const& pointToBeSeparated) const; + virtual std::pair optimize(Point const& direction) const; /* * Returns a string representation of this polytope.