some minor tweaks plus polling example

Former-commit-id: eebe4ca6d6

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

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);

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) {

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);