diff --git a/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h
index 60a5de999..768237cca 100644
--- a/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/GmmxxMdpPrctlModelChecker.h
@@ -102,7 +102,7 @@ private:
 	 * @return The solution of the system of linear equations in form of the elements of the vector
 	 * x.
 	 */
-	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<uint_fast64_t>* takenChoices = nullptr) const override {
 		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -145,6 +145,11 @@ private:
 			newX = swap;
 			++iterations;
 		}
+        
+        // If we were requested to record the taken choices, we have to construct the vector now.
+        if (takenChoices != nullptr) {
+            this->computeTakenChoices(multiplyResult, *takenChoices, nondeterministicChoiceIndices);
+        }
 
 		// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
 		// is currently stored in currentX, but x is the output vector.
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 2f7238a81..1bac8ad37 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -227,7 +227,7 @@ public:
 		return result;
 	}
 
-	/*!
+    /*!
 	 * Check the given formula that is an until formula.
 	 *
 	 * @param formula The formula to check.
@@ -235,10 +235,31 @@ public:
 	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
 	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
 	 * bounds 0 and 1.
-	 * @returns The probabilities for the given formula to hold on every state of the model associated with this model
+     * @param scheduler If <code>qualitative</code> is false and this vector is non-null and has as many elements as
+     * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability
+     * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state.
+	 * @return The probabilities for the given formula to hold on every state of the model associated with this model
 	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
 	 */
 	virtual std::vector<Type>* checkUntil(const storm::property::prctl::Until<Type>& formula, bool qualitative) const {
+        return this->checkUntil(formula, qualitative, nullptr);
+    }
+    
+	/*!
+	 * Check the given formula that is an until formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bounds 0 and 1. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bounds 0 and 1.
+     * @param scheduler If <code>qualitative</code> is false and this vector is non-null and has as many elements as
+     * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability
+     * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state.
+	 * @return The probabilities for the given formula to hold on every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact probabilities might not be computed.
+	 */
+	std::vector<Type>* checkUntil(const storm::property::prctl::Until<Type>& formula, bool qualitative, std::vector<uint_fast64_t>* scheduler) const {
 		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
 		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
 		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
@@ -378,7 +399,7 @@ public:
 		return result;
 	}
 
-	/*!
+    /*!
 	 * Checks the given formula that is a reachability reward formula.
 	 *
 	 * @param formula The formula to check.
@@ -386,10 +407,28 @@ public:
 	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
 	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
 	 * bound 0.
-	 * @returns The reward values for the given formula for every state of the model associated with this model
+	 * @return The reward values for the given formula for every state of the model associated with this model
 	 * checker. If the qualitative flag is set, exact values might not be computed.
 	 */
 	virtual std::vector<Type>* checkReachabilityReward(const storm::property::prctl::ReachabilityReward<Type>& formula, bool qualitative) const {
+        return this->checkReachabilityReward(formula, qualitative, nullptr);
+    }
+    
+	/*!
+	 * Checks the given formula that is a reachability reward formula.
+	 *
+	 * @param formula The formula to check.
+	 * @param qualitative A flag indicating whether the formula only needs to be evaluated qualitatively, i.e. if the
+	 * results are only compared against the bound 0. If set to true, this will most likely results that are only
+	 * qualitatively correct, i.e. do not represent the correct value, but only the correct relation with respect to the
+	 * bound 0.
+     * @param scheduler If <code>qualitative</code> is false and this vector is non-null and has as many elements as
+     * there are states in the MDP, this vector will represent a scheduler for the model that achieves the probability
+     * returned by model checking. To this end, the vector will hold the nondeterministic choice made for each state.
+	 * @return The reward values for the given formula for every state of the model associated with this model
+	 * checker. If the qualitative flag is set, exact values might not be computed.
+	 */
+	std::vector<Type>* checkReachabilityReward(const storm::property::prctl::ReachabilityReward<Type>& formula, bool qualitative, std::vector<uint_fast64_t>* scheduler) const {
 		// Only compute the result if the model has at least one reward model.
 		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
 			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
@@ -465,7 +504,7 @@ public:
 			}
 
 			// Solve the corresponding system of equations.
-			this->solveEquationSystem(submatrix, x, b, subNondeterministicChoiceIndices);
+			this->solveEquationSystem(submatrix, x, b, subNondeterministicChoiceIndices, scheduler);
 
 			// Set values of resulting vector according to result.
 			storm::utility::vector::setVectorValues<Type>(*result, maybeStates, x);
@@ -481,6 +520,22 @@ public:
 	}
 
 protected:
+    /*!
+     * Computes the vector of choices that need to be made to minimize/maximize the model checking result for each state.
+     *
+     * @param nondeterministicResult The model checking result for nondeterministic choices of all states.
+     * @param takenChoices The output vector that is to store the taken choices.
+     * @param nondeterministicChoiceIndices The assignment of states to their nondeterministic choices in the matrix.
+     */
+    void computeTakenChoices(std::vector<Type> const& nondeterministicResult, std::vector<uint_fast64_t>& takenChoices, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+        std::vector<Type> temporaryResult(nondeterministicChoiceIndices.size() - 1);
+        if (this->minimumOperatorStack.top()) {
+            storm::utility::vector::reduceVectorMin(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices);
+        } else {
+            storm::utility::vector::reduceVectorMax(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &takenChoices);
+        }
+    }
+    
 	/*!
 	 * A stack used for storing whether we are currently computing min or max probabilities or rewards, respectively.
 	 * The topmost element is true if and only if we are currently computing minimum probabilities or rewards.
@@ -532,9 +587,12 @@ private:
 	 * may be ignored.
 	 * @param b The right-hand side of the equation system.
      * @param nondeterministicChoiceIndices The assignment of states to their rows in the matrix.
+     * @param takenChoices If not null, this vector will be filled with the nondeterministic choices taken by the states
+     * to achieve the solution of the equation system. This assumes that the given vector has at least as many elements
+     * as there are states in the MDP.
 	 * @returns The solution vector x of the system of linear equations as the content of the parameter x.
 	 */
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<uint_fast64_t>* takenChoices = nullptr) const {
 		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -575,6 +633,11 @@ private:
 			newX = swap;
 			++iterations;
 		}
+        
+        // If we were requested to record the taken choices, we have to construct the vector now.
+        if (takenChoices != nullptr) {
+            this->computeTakenChoices(multiplyResult, *takenChoices, nondeterministicChoiceIndices);
+        }
 
 		// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
 		// is currently stored in currentX, but x is the output vector.
@@ -592,7 +655,7 @@ private:
 			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 		}
 	}
-
+    
 	/*!
 	 * Computes the nondeterministic choice indices vector resulting from reducing the full system to the states given
 	 * by the parameter constraint.
diff --git a/src/models/AbstractDeterministicModel.h b/src/models/AbstractDeterministicModel.h
index 9cb6990e9..c09f42f90 100644
--- a/src/models/AbstractDeterministicModel.h
+++ b/src/models/AbstractDeterministicModel.h
@@ -71,7 +71,9 @@ class AbstractDeterministicModel: public AbstractModel<T> {
             
             for (auto const& transition : *this->transitionMatrix) {
                 if (transition.value() != storm::utility::constGetZero<T>()) {
-                    outStream << "\t" << transition.row() << " -> " << transition.column() << " [ label= \"" << transition.value() << "\" ];" << std::endl;
+                    if (subsystem == nullptr || subsystem->get(transition.column())) {
+                        outStream << "\t" << transition.row() << " -> " << transition.column() << " [ label= \"" << transition.value() << "\" ];" << std::endl;
+                    }
                 }
             }
                         
diff --git a/src/models/AbstractModel.h b/src/models/AbstractModel.h
index b1233f1c0..31c019e66 100644
--- a/src/models/AbstractModel.h
+++ b/src/models/AbstractModel.h
@@ -347,59 +347,61 @@ protected:
             outStream << "digraph deterministicModel {" << std::endl;
         
             for (uint_fast64_t state = 0, highestStateIndex = this->getNumberOfStates() - 1; state <= highestStateIndex; ++state) {
-                outStream << "\t" << state;
-                if (includeLabeling || firstValue != nullptr || secondValue != nullptr || stateColoring != nullptr) {
-                    outStream << " [ ";
-                    if (includeLabeling || firstValue != nullptr || secondValue != nullptr) {
-                        outStream << "label = \"" << state << ": ";
+                if (subsystem == nullptr || subsystem->get(state)) {
+                    outStream << "\t" << state;
+                    if (includeLabeling || firstValue != nullptr || secondValue != nullptr || stateColoring != nullptr) {
+                        outStream << " [ ";
+                        if (includeLabeling || firstValue != nullptr || secondValue != nullptr) {
+                            outStream << "label = \"" << state << ": ";
                     
-                        // Now print the state labeling to the stream if requested.
-                        if (includeLabeling) {
-                            outStream << "{";
-                            bool includeComma = false;
-                            for (std::string const& label : this->getLabelsForState(state)) {
-                                if (includeComma) {
-                                    outStream << ", ";
-                                } else {
-                                    includeComma = true;
+                            // Now print the state labeling to the stream if requested.
+                            if (includeLabeling) {
+                                outStream << "{";
+                                bool includeComma = false;
+                                for (std::string const& label : this->getLabelsForState(state)) {
+                                    if (includeComma) {
+                                        outStream << ", ";
+                                    } else {
+                                        includeComma = true;
+                                    }
+                                    outStream << label;
                                 }
-                                outStream << label;
+                                outStream << "}";
                             }
-                            outStream << "}";
-                        }
                     
-                        // If we are to include some values for the state as well, we do so now.
-                        if (firstValue != nullptr || secondValue != nullptr) {
-                            outStream << "[";
-                            if (firstValue != nullptr) {
-                                outStream << (*firstValue)[state];
+                            // If we are to include some values for the state as well, we do so now.
+                            if (firstValue != nullptr || secondValue != nullptr) {
+                                outStream << "[";
+                                if (firstValue != nullptr) {
+                                    outStream << (*firstValue)[state];
+                                    if (secondValue != nullptr) {
+                                        outStream << ", ";
+                                    }
+                                }
                                 if (secondValue != nullptr) {
-                                    outStream << ", ";
+                                    outStream << (*secondValue)[state];
                                 }
+                                outStream << "]";
                             }
-                            if (secondValue != nullptr) {
-                                outStream << (*secondValue)[state];
-                            }
-                            outStream << "]";
-                        }
-                        outStream << "\"";
+                            outStream << "\"";
                     
-                        // Now, we color the states if there were colors given.
-                        if (stateColoring != nullptr && colors != nullptr) {
-                            outStream << ", ";
-                            outStream << " fillcolor = \"" << (*colors)[(*stateColoring)[state]] << "\"";
+                            // Now, we color the states if there were colors given.
+                            if (stateColoring != nullptr && colors != nullptr) {
+                                outStream << ", ";
+                                outStream << " fillcolor = \"" << (*colors)[(*stateColoring)[state]] << "\"";
+                            }
                         }
+                        outStream << " ]";
                     }
-                    outStream << " ]";
+                    outStream << ";" << std::endl;
                 }
-                outStream << ";" << std::endl;
             }
             
             if (finalizeOutput) {
                 outStream << "}" << std::endl;
             }
         }
-
+        
 		/*! A matrix representing the likelihoods of moving between states. */
 		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix;
 
diff --git a/src/models/AbstractNondeterministicModel.h b/src/models/AbstractNondeterministicModel.h
index 39b4a6ed2..e9f101315 100644
--- a/src/models/AbstractNondeterministicModel.h
+++ b/src/models/AbstractNondeterministicModel.h
@@ -144,17 +144,19 @@ class AbstractNondeterministicModel: public AbstractModel<T> {
                     // Now draw all probabilitic arcs that belong to this nondeterminstic choice.
                     transitionIte.moveToNextRow();
                     for (; transitionIt != transitionIte; ++transitionIt) {
-                        outStream << "\t\"" << state << "c" << row << "\" -> " << transitionIt.column() << " [ label= \"" << transitionIt.value() << "\" ]";
+                        if (subsystem == nullptr || subsystem->get(transitionIt.column())) {
+                            outStream << "\t\"" << state << "c" << row << "\" -> " << transitionIt.column() << " [ label= \"" << transitionIt.value() << "\" ]";
                         
-                        // If we were given a scheduler to highlight, we do so now.
-                        if (scheduler != nullptr) {
-                            if (highlightChoice) {
-                                outStream << " [color=\"red\", penwidth = 2]";
-                            } else {
-                                outStream << " [style = \"dotted\"]";
+                            // If we were given a scheduler to highlight, we do so now.
+                            if (scheduler != nullptr) {
+                                if (highlightChoice) {
+                                    outStream << " [color=\"red\", penwidth = 2]";
+                                } else {
+                                    outStream << " [style = \"dotted\"]";
+                                }
                             }
+                            outStream << ";" << std::endl;
                         }
-                        outStream << ";" << std::endl;
                     }
                 }
             }
diff --git a/src/utility/vector.h b/src/utility/vector.h
index fff25fd52..9b973de04 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -155,16 +155,22 @@ void addVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
  * @param rowGrouping A vector that specifies the begin and end of each group of elements in the values vector.
  * @param filter A function that compares two elements v1 and v2 according to some filter criterion. This function must
  * return true iff v1 is supposed to be taken instead of v2.
+ * @param choices If non-null, this vector is used to store the choices made during the selection.
  */
 template<class T>
-void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::function<bool (T const&, T const&)> filter) {
+void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::function<bool (T const&, T const&)> filter, std::vector<uint_fast64_t>* choices = nullptr) {
     uint_fast64_t currentSourceRow = 0;
     uint_fast64_t currentTargetRow = -1;
-    for (auto it = source.cbegin(), ite = source.cend(); it != ite; ++it, ++currentSourceRow) {
+    uint_fast64_t currentLocalRow = 0;
+    for (auto it = source.cbegin(), ite = source.cend(); it != ite; ++it, ++currentSourceRow, ++currentLocalRow) {
         // Check whether we have considered all source rows for the current target row.
         if (rowGrouping[currentTargetRow + 1] <= currentSourceRow || currentSourceRow == 0) {
+            currentLocalRow = 0;
             ++currentTargetRow;
             target[currentTargetRow] = source[currentSourceRow];
+            if (choices != nullptr) {
+                (*choices)[currentTargetRow] = 0;
+            }
             continue;
         }
             
@@ -172,6 +178,9 @@ void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vec
         // value is actually lower.
         if (filter(*it, target[currentTargetRow])) {
             target[currentTargetRow] = *it;
+            if (choices != nullptr) {
+                (*choices)[currentTargetRow] = currentLocalRow;
+            }
         }
     }
 }
@@ -181,10 +190,11 @@ void reduceVector(std::vector<T> const& source, std::vector<T>& target, std::vec
  *
  * @param source The source vector which is to be reduced.
  * @param target The target vector into which a single element from each row group is written.
- * @param rowGrouping A vector that specifies the begin and end of each group of elements in the values vector.
+ * @param rowGrouping A vector that specifies the begin and end of each group of elements in the source vector.
+ * @param choices If non-null, this vector is used to store the choices made during the selection.
  */
 template<class T>
-void reduceVectorMin(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping) {
+void reduceVectorMin(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
 	reduceVector<T>(source, target, rowGrouping, std::less<T>());
 }
 
@@ -193,13 +203,14 @@ void reduceVectorMin(std::vector<T> const& source, std::vector<T>& target, std::
  *
  * @param source The source vector which is to be reduced.
  * @param target The target vector into which a single element from each row group is written.
- * @param rowGrouping A vector that specifies the begin and end of each group of elements in the values vector.
+ * @param rowGrouping A vector that specifies the begin and end of each group of elements in the source vector.
+ * @param choices If non-null, this vector is used to store the choices made during the selection.
  */
 template<class T>
-void reduceVectorMax(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping) {
+void reduceVectorMax(std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
     reduceVector<T>(source, target, rowGrouping, std::greater<T>());
 }
-
+    
 /*!
  * Compares the given elements and determines whether they are equal modulo the given precision. The provided flag
  * additionaly specifies whether the error is computed in relative or absolute terms.