diff --git a/examples/multiobjective/ma/polling/polling.ma b/examples/multiobjective/ma/polling/polling.ma
new file mode 100644
index 000000000..ef9c4958e
--- /dev/null
+++ b/examples/multiobjective/ma/polling/polling.ma
@@ -0,0 +1,105 @@
+// Translation of the MAPA Specification of a polling system into PRISM code
+// http://wwwhome.cs.utwente.nl/~timmer/scoop/papers/qest13/index.html
+
+ma
+
+const int J; // number of job types (should be in [1..6])
+const int Q; // Maximum queue size in each station
+const bool samerate=false; // true if all stations should have the same rate
+
+// Formulae to control the LIFO queue of the stations.
+// The queue is represented by some integer whose base J representation has at most Q digits, each representing one of the job types 0, 1, ..., J-1.
+// In addition, we store the current size of the queue which is needed to distinguish an empty queue from a queue holding job of type 0
+formula queue1_empty = q1Size=0;
+formula queue1_full = q1Size=Q;
+formula queue1_pop = floor(q1/J);
+formula queue1_head = q1 - (queue1_pop * J); // i.e. q1 modulo J
+formula queue1_push = q1*J;
+formula queue2_empty = q2Size=0;
+formula queue2_full = q2Size=Q;
+formula queue2_pop = floor(q2/J);
+formula queue2_head = q2 - (queue2_pop * J); // i.e. q2 modulo J
+formula queue2_push = q2*J;
+
+const int queue_maxValue = (J^Q)-1;
+
+
+
+const double inRate1 = 3; // = (2*1) + 1;
+const double inRate2 = 5; // = (2*2) + 1;
+const double inRate3 = 7; // = (2*3) + 1;
+
+module pollingsys
+ // The queues for the stations
+ q1 : [0..queue_maxValue];
+ q1Size : [0..Q];
+ q2 : [0..queue_maxValue];
+ q2Size : [0..Q];
+
+ // Store the job that is currently processed by the server. j=J means that no job is processed.
+ j : [0..J] init J;
+
+ // Flag indicating whether a new job arrived
+ newJob1 : bool init false;
+ newJob2 : bool init false;
+
+ //<> !newJob1 & !newJob2 & !queue1_full & queue2_full & j=J -> inRate1 : (newJob1'=true);
+ //<> !newJob1 & !newJob2 & queue1_full & !queue2_full & j=J -> inRate2 : (newJob2'=true);
+ <> !newJob1 & !newJob2 & !queue1_full & !queue2_full & j=J -> inRate1 : (newJob1'=true) + inRate2 : (newJob2'=true);
+ <> !newJob1 & !newJob2 & queue1_full & queue2_full & j<J -> 2*(j+1) : (j'=J);
+ <> !newJob1 & !newJob2 & !queue1_full & queue2_full & j<J -> inRate1 : (newJob1'=true) + 2*(j+1) : (j'=J);
+ <> !newJob1 & !newJob2 & queue1_full & !queue2_full & j<J -> inRate2 : (newJob2'=true) + 2*(j+1) : (j'=J);
+ <> !newJob1 & !newJob2 & !queue1_full & !queue2_full & j<J -> inRate1 : (newJob1'=true) + inRate2 : (newJob2'=true) + 2*(j+1) : (j'=J);
+
+ [] newJob1 & J>=1 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+0) & (newJob1'=false);
+ [] newJob1 & J>=2 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+1) & (newJob1'=false);
+ [] newJob1 & J>=3 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+2) & (newJob1'=false);
+ [] newJob1 & J>=4 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+3) & (newJob1'=false);
+ [] newJob1 & J>=5 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+4) & (newJob1'=false);
+ [] newJob1 & J>=6 -> 1 : (q1Size'=q1Size+1) & (q1'=queue1_push+5) & (newJob1'=false);
+
+ [] newJob2 & J>=1 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+0) & (newJob2'=false);
+ [] newJob2 & J>=2 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+1) & (newJob2'=false);
+ [] newJob2 & J>=3 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+2) & (newJob2'=false);
+ [] newJob2 & J>=4 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+3) & (newJob2'=false);
+ [] newJob2 & J>=5 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+4) & (newJob2'=false);
+ [] newJob2 & J>=6 -> 1 : (q2Size'=q2Size+1) & (q2'=queue2_push+5) & (newJob2'=false);
+
+ [copy1] !newJob1 & !newJob2 & !queue1_empty & j=J -> 0.9 : (j'=queue1_head) & (q1Size'=q1Size-1) & (q1'=queue1_pop) + 0.1 : (j'=queue1_head);
+ [copy2] !newJob1 & !newJob2 & !queue2_empty & j=J -> 0.9 : (j'=queue2_head) & (q2Size'=q2Size-1) & (q2'=queue2_pop) + 0.1 : (j'=queue2_head);
+
+endmodule
+
+
+
+label "q1full" = q1Size=Q;
+label "q2full" = q2Size=Q;
+label "allqueuesfull" = q1Size=Q & q2Size=Q;
+
+
+// Rewards adapted from Guck et al.: Modelling and Analysis of Markov Reward Automata
+
+rewards "processedjobs1"
+ [copy1] true : 1;
+endrewards
+
+rewards "processedjobs2"
+ [copy1] true : 1;
+endrewards
+
+rewards "processedjobs"
+ [copy1] true : 1;
+ [copy2] true : 1;
+endrewards
+
+rewards "waiting1"
+ true : (q1Size);
+endrewards
+
+rewards "waiting2"
+ true : (q2Size);
+endrewards
+
+rewards "waiting"
+ true : (q1Size + q2Size);
+endrewards
\ No newline at end of file
diff --git a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
index af76131ba..a891f1b1c 100644
--- a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
+++ b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp
@@ -78,10 +78,10 @@ namespace storm {
// Only perform such a step if there is time left.
if(currentEpoch>0) {
performMSStep(MS, PS, consideredObjectives);
- if(currentEpoch % 1000 == 0) {
- STORM_LOG_DEBUG(currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0]);
- std::cout << std::endl << currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0] << std::endl;
- }
+ // if(currentEpoch % 1000 == 0) {
+ // STORM_LOG_DEBUG(currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0]);
+ // std::cout << std::endl << currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0] << std::endl;
+ // }
--currentEpoch;
} else {
break;
@@ -261,17 +261,19 @@ namespace storm {
VT const maxRate = this->data.preprocessedModel.getMaximalExitRate();
for(uint_fast64_t objIndex = 0; objIndex < this->data.objectives.size(); ++objIndex) {
auto const& obj = this->data.objectives[objIndex];
+ VT errorTowardsZero;
+ VT errorAwayFromZero;
if(obj.lowerTimeBound) {
uint_fast64_t digitizedBound = storm::utility::convertNumber<uint_fast64_t>((*obj.lowerTimeBound)/digitizationConstant);
auto timeBoundIt = lowerTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->data.objectives.size(), false))).first;
timeBoundIt->second.set(objIndex);
VT digitizationError = storm::utility::one<VT>();
digitizationError -= std::exp(-maxRate * (*obj.lowerTimeBound)) * storm::utility::pow(storm::utility::one<VT>() + maxRate * digitizationConstant, digitizedBound);
- this->offsetsToLowerBound[objIndex] = -digitizationError;
- this->offsetsToUpperBound[objIndex] = storm::utility::one<VT>() - std::exp(-maxRate * digitizationConstant);;
+ errorTowardsZero = -digitizationError;
+ errorAwayFromZero = storm::utility::one<VT>() - std::exp(-maxRate * digitizationConstant);;
} else {
- this->offsetsToLowerBound[objIndex] = storm::utility::zero<VT>();
- this->offsetsToUpperBound[objIndex] = storm::utility::zero<VT>();
+ errorTowardsZero = storm::utility::zero<VT>();
+ errorAwayFromZero = storm::utility::zero<VT>();
}
if(obj.upperTimeBound) {
uint_fast64_t digitizedBound = storm::utility::convertNumber<uint_fast64_t>((*obj.upperTimeBound)/digitizationConstant);
@@ -279,9 +281,16 @@ namespace storm {
timeBoundIt->second.set(objIndex);
VT digitizationError = storm::utility::one<VT>();
digitizationError -= std::exp(-maxRate * (*obj.upperTimeBound)) * storm::utility::pow(storm::utility::one<VT>() + maxRate * digitizationConstant, digitizedBound);
- this->offsetsToUpperBound[objIndex] += digitizationError;
+ errorAwayFromZero += digitizationError;
+ }
+ STORM_LOG_ASSERT(errorTowardsZero + errorAwayFromZero <= this->maximumLowerUpperBoundGap, "Precision not sufficient.");
+ if (obj.rewardsArePositive) {
+ this->offsetsToLowerBound[objIndex] = -errorTowardsZero;
+ this->offsetsToUpperBound[objIndex] = errorAwayFromZero;
+ } else {
+ this->offsetsToLowerBound[objIndex] = -errorAwayFromZero;
+ this->offsetsToUpperBound[objIndex] = errorTowardsZero;
}
- STORM_LOG_ASSERT(this->offsetsToUpperBound[objIndex] - this->offsetsToLowerBound[objIndex] <= this->maximumLowerUpperBoundGap, "Precision not sufficient.");
}
}
@@ -364,14 +373,14 @@ namespace storm {
// check whether the linEqSolver needs to be updated, i.e., whether the scheduler has changed
if(newOptimalChoices != optimalChoicesAtCurrentEpoch) {
- std::cout << "X";
+ // std::cout << "X";
optimalChoicesAtCurrentEpoch.swap(newOptimalChoices);
linEq.solver = nullptr;
storm::storage::SparseMatrix<ValueType> linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true);
linEqMatrix.convertToEquationSystem();
linEq.solver = linEq.factory.create(std::move(linEqMatrix));
} else {
- std::cout << " ";
+ // std::cout << " ";
}
for(auto objIndex : consideredObjectives) {
PS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], PS.auxChoiceValues);
diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
index 02a05bd3d..627db4eba 100644
--- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
+++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveHelper.cpp
@@ -57,8 +57,8 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::achievabilityQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
- //Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type.
- weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1)); // TODO try other values?
+ //Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type. Note that we could also try other values
+ weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1));
// Get a point that represents the thresholds
Point thresholds;
thresholds.reserve(preprocessorData.objectives.size());
@@ -84,8 +84,8 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::numericalQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
STORM_LOG_ASSERT(preprocessorData.indexOfOptimizingObjective, "Detected numerical query but index of optimizing objective is not set.");
- // Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type.
- weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1)); // TODO try other values?
+ // Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type. Note that we could also try other values
+ weightVectorChecker->setMaximumLowerUpperBoundGap(storm::utility::convertNumber<SparseModelValueType>(0.1));
// initialize some data
uint_fast64_t optimizingObjIndex = *preprocessorData.indexOfOptimizingObjective;
Point thresholds;
@@ -131,11 +131,12 @@ namespace storm {
// Note that we do not have to care whether a threshold is strict anymore, because the resulting optimum should be
// the supremum over all strategies. Hence, one could combine a scheduler inducing the optimum value (but possibly violating strict
// thresholds) and (with very low probability) a scheduler that satisfies all (possibly strict) thresholds.
- // The euclidean distance between the lower and upper bounds of the results of the weightVectorChecker should be less than the precision given in the settings
+
+ // Pick the maximum gap between lower and upper bound of the weightVectorChecker result.
SparseModelValueType gap = storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
- gap -= storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
- gap /= storm::utility::sqrt(static_cast<SparseModelValueType>(preprocessorData.objectives.size()));
- weightVectorChecker->setMaximumLowerUpperBoundGap(gap); // TODO try other values?
+ gap /= (storm::utility::one<SparseModelValueType>() + storm::utility::one<SparseModelValueType>());
+ weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
+
while(true) {
std::pair<Point, bool> optimizationRes = resultData.underApproximation()->intersection(thresholdsAsPolytope)->optimize(directionOfOptimizingObjective);
STORM_LOG_THROW(optimizationRes.second, storm::exceptions::UnexpectedException, "The underapproximation is either unbounded or empty.");
@@ -160,12 +161,13 @@ namespace storm {
template <class SparseModelType, typename RationalNumberType>
void SparseMultiObjectiveHelper<SparseModelType, RationalNumberType>::paretoQuery(PreprocessorData const& preprocessorData, WeightVectorCheckerType weightVectorChecker, ResultData& resultData) {
- // Set the maximum gap between lower and upper bound of the weightVectorChecker result we initially allow for this query type.
- // The euclidean distance between the lower and upper bounds of the results of the weightVectorChecker should be less than the precision given in the settings
+ // Set the maximum gap between lower and upper bound of the weightVectorChecker result.
+ // This is the maximal edge length of the box we have to consider around each computed point
+ // We pick the gap such that the maximal distance between two points within this box is less than the given precision divided by two.
SparseModelValueType gap = storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::MultiObjectiveSettings>().getPrecision());
- gap -= storm::utility::convertNumber<SparseModelValueType>(storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision());
+ gap /= (storm::utility::one<SparseModelValueType>() + storm::utility::one<SparseModelValueType>());
gap /= storm::utility::sqrt(static_cast<SparseModelValueType>(preprocessorData.objectives.size()));
- weightVectorChecker->setMaximumLowerUpperBoundGap(gap); // TODO try other values?
+ weightVectorChecker->setMaximumLowerUpperBoundGap(gap);
//First consider the objectives individually
for(uint_fast64_t objIndex = 0; objIndex<preprocessorData.objectives.size() && !maxStepsPerformed(resultData); ++objIndex) {
diff --git a/src/utility/storm.h b/src/utility/storm.h
index 502cdd215..c78906434 100644
--- a/src/utility/storm.h
+++ b/src/utility/storm.h
@@ -327,11 +327,11 @@ namespace storm {
template<typename ValueType>
std::unique_ptr<storm::modelchecker::CheckResult> verifySparseMarkovAutomaton(std::shared_ptr<storm::models::sparse::MarkovAutomaton<ValueType>> ma, storm::modelchecker::CheckTask<storm::logic::Formula> const& task) {
std::unique_ptr<storm::modelchecker::CheckResult> result;
+ // Close the MA, if it is not already closed.
+ if (!ma->isClosed()) {
+ ma->close();
+ }
if(task.getFormula().isMultiObjectiveFormula()) {
- // Close the MA, if it is not already closed.
- if (!ma->isClosed()) {
- ma->close();
- }
storm::modelchecker::SparseMaMultiObjectiveModelChecker<storm::models::sparse::MarkovAutomaton<ValueType>> modelchecker(*ma);
if (modelchecker.canHandle(task)) {
result = modelchecker.check(task);