diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd9d84115..83f581951 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -32,6 +32,7 @@ option(USE_INTELTBB "Sets whether the Intel TBB Extensions should be used." OFF)
 option(STORM_USE_COTIRE "Sets whether Cotire should be used (for building precompiled headers)." OFF)
 option(LINK_LIBCXXABI "Sets whether libc++abi should be linked." OFF)
 set(GUROBI_ROOT "" CACHE STRING "The root directory of Gurobi (if available).")
+set(Z3_ROOT "" CACHE STRING "The root directory of Z3 (if available).")
 set(ADDITIONAL_INCLUDE_DIRS "" CACHE STRING "Additional directories added to the include directories.")
 set(ADDITIONAL_LINK_DIRS "" CACHE STRING "Additional directories added to the link directories.")
 set(CUSTOM_BOOST_ROOT "" CACHE STRING "A custom path to the Boost root directory.")
@@ -84,6 +85,12 @@ else()
     set(ENABLE_GUROBI ON)
 endif()
 
+if ("${Z3_ROOT}" STREQUAL "")
+    set(ENABLE_Z3 OFF)
+else()
+    set(ENABLE_Z3 ON)
+endif()
+
 message(STATUS "StoRM - CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}")
 message(STATUS "StoRM - CMAKE_BUILD_TYPE (ENV): $ENV{CMAKE_BUILD_TYPE}")
 
@@ -171,6 +178,13 @@ else()
 	set(STORM_CPP_GUROBI_DEF "undef")
 endif()
 
+# Z3 Defines
+if (ENABLE_Z3)
+	set(STORM_CPP_Z3_DEF "define")
+else()
+	set(STORM_CPP_Z3_DEF "undef")
+endif()
+
 # Intel TBB Defines
 if (TBB_FOUND AND USE_INTELTBB)
 	set(STORM_CPP_INTELTBB_DEF "define")
diff --git a/src/adapters/ExplicitModelAdapter.h b/src/adapters/ExplicitModelAdapter.h
index 20d2c1133..a38b7c27b 100644
--- a/src/adapters/ExplicitModelAdapter.h
+++ b/src/adapters/ExplicitModelAdapter.h
@@ -63,7 +63,7 @@ namespace storm {
         };
         
         /*!
-         * A helper class that provides the functionality to compare states of the model.
+         * A helper class that provides the functionality to compare states of the model wrt. equality.
          */
         class StateCompare {
         public:
@@ -72,8 +72,30 @@ namespace storm {
             }
         };
         
+        /*!
+         * A helper class that provides the functionality to compare states of the model wrt. the relation "<".
+         */
+        class StateLess {
+        public:
+            bool operator()(StateType* state1, StateType* state2) const {
+                // Compare boolean variables.
+                for (uint_fast64_t i = 0; i < state1->first.size(); ++i) {
+                    if (!state1->first.at(i) && state2->first.at(i)) {
+                        return true;
+                    }
+                }
+                // Then compare integer variables.
+                for (uint_fast64_t i = 0; i < state1->second.size(); ++i) {
+                    if (!state1->second.at(i) && state2->second.at(i)) {
+                        return true;
+                    }
+                }
+                return false;
+            }
+        };
+        
         // A structure holding information about a particular choice.
-        template<typename ValueType>
+        template<typename ValueType, typename KeyType=uint_fast64_t, typename Compare=std::less<uint_fast64_t>>
         struct Choice {
         public:
             Choice(std::string const& actionLabel) : distribution(), actionLabel(actionLabel), choiceLabels() {
@@ -85,7 +107,7 @@ namespace storm {
              *
              * @return An iterator to the first element of this choice.
              */
-            typename std::map<uint_fast64_t, ValueType>::iterator begin() {
+            typename std::map<KeyType, ValueType>::iterator begin() {
                 return distribution.begin();
             }
             
@@ -94,7 +116,7 @@ namespace storm {
              *
              * @return An iterator to the first element of this choice.
              */
-            typename std::map<uint_fast64_t, ValueType>::const_iterator begin() const {
+            typename std::map<KeyType, ValueType>::const_iterator begin() const {
                 return distribution.cbegin();
             }
             
@@ -103,7 +125,7 @@ namespace storm {
              *
              * @return An iterator that points past the elements of this choice.
              */
-            typename std::map<uint_fast64_t, ValueType>::iterator end() {
+            typename std::map<KeyType, ValueType>::iterator end() {
                 return distribution.end();
             }
             
@@ -112,7 +134,7 @@ namespace storm {
              *
              * @return An iterator that points past the elements of this choice.
              */
-            typename std::map<uint_fast64_t, ValueType>::const_iterator end() const {
+            typename std::map<KeyType, ValueType>::const_iterator end() const {
                 return distribution.cend();
             }
             
@@ -123,7 +145,7 @@ namespace storm {
              * @param value The value to find.
              * @return An iterator to the element with the given key, if there is one.
              */
-            typename std::map<uint_fast64_t, ValueType>::iterator find(uint_fast64_t value) {
+            typename std::map<KeyType, ValueType>::iterator find(uint_fast64_t value) {
                 return distribution.find(value);
             }
             
@@ -213,7 +235,7 @@ namespace storm {
             
         private:
             // The distribution that is associated with the choice.
-            std::map<uint_fast64_t, ValueType> distribution;
+            std::map<KeyType, ValueType, Compare> distribution;
             
             // The label of the choice.
             std::string actionLabel;
@@ -257,7 +279,7 @@ namespace storm {
         public:
             // A structure holding information about the reachable state space.
             struct StateInformation {
-                StateInformation() : reachableStates(), stateToIndexMap(), numberOfChoices(), numberOfTransitions(0) {
+                StateInformation() : reachableStates(), stateToIndexMap() {
                     // Intentionally left empty.
                 }
                 
@@ -266,12 +288,6 @@ namespace storm {
                 
                 // A mapping from states to indices in the list of reachable states.
                 std::unordered_map<StateType*, uint_fast64_t, StateHash, StateCompare> stateToIndexMap;
-                
-                // A vector storing the number of choices for each state.
-                std::vector<uint_fast64_t> numberOfChoices;
-                
-                // The total number of transitions discovered.
-                uint_fast64_t numberOfTransitions;
             };
             
             // A structure holding the individual components of a model.
@@ -339,16 +355,16 @@ namespace storm {
                 std::shared_ptr<storm::models::AbstractModel<ValueType>> result;
                 switch (program.getModelType()) {
                     case storm::ir::Program::DTMC:
-                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Dtmc<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.stateRewards), std::move(modelComponents.transitionRewardMatrix), std::move(modelComponents.choiceLabeling)));
+                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Dtmc<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), rewardModelName != "" ? std::move(modelComponents.stateRewards) : boost::optional<std::vector<ValueType>>(), rewardModelName != "" ? std::move(modelComponents.transitionRewardMatrix) : boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(modelComponents.choiceLabeling)));
                         break;
                     case storm::ir::Program::CTMC:
-                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Ctmc<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.stateRewards), std::move(modelComponents.transitionRewardMatrix), std::move(modelComponents.choiceLabeling)));
+                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Ctmc<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), rewardModelName != "" ? std::move(modelComponents.stateRewards) : boost::optional<std::vector<ValueType>>(), rewardModelName != "" ? std::move(modelComponents.transitionRewardMatrix) : boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(modelComponents.choiceLabeling)));
                         break;
                     case storm::ir::Program::MDP:
-                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Mdp<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.nondeterministicChoiceIndices), std::move(modelComponents.stateRewards), std::move(modelComponents.transitionRewardMatrix), std::move(modelComponents.choiceLabeling)));
+                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Mdp<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.nondeterministicChoiceIndices), rewardModelName != "" ? std::move(modelComponents.stateRewards) : boost::optional<std::vector<ValueType>>(), rewardModelName != "" ? std::move(modelComponents.transitionRewardMatrix) : boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(modelComponents.choiceLabeling)));
                         break;
                     case storm::ir::Program::CTMDP:
-                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Ctmdp<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.nondeterministicChoiceIndices), std::move(modelComponents.stateRewards), std::move(modelComponents.transitionRewardMatrix), std::move(modelComponents.choiceLabeling)));
+                        result = std::shared_ptr<storm::models::AbstractModel<ValueType>>(new storm::models::Ctmdp<ValueType>(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.nondeterministicChoiceIndices), rewardModelName != "" ? std::move(modelComponents.stateRewards) : boost::optional<std::vector<ValueType>>(), rewardModelName != "" ? std::move(modelComponents.transitionRewardMatrix) : boost::optional<storm::storage::SparseMatrix<ValueType>>(), std::move(modelComponents.choiceLabeling)));
                         break;
                     default:
                         LOG4CPLUS_ERROR(logger, "Error while creating model from probabilistic program: cannot handle this model type.");
@@ -572,7 +588,6 @@ namespace storm {
                 if (indexIt == stateInformation.stateToIndexMap.end()) {
                     // The state has not been seen, yet, so add it to the list of all reachable states.
                     stateInformation.stateToIndexMap[state] = stateInformation.reachableStates.size();
-                    stateInformation.numberOfChoices.push_back(0);
                     stateInformation.reachableStates.push_back(state);
                     return std::make_pair(false, stateInformation.stateToIndexMap[state]);
                 } else {
@@ -581,59 +596,7 @@ namespace storm {
                     return std::make_pair(true, indexIt->second);
                 }
             }
-            
-            /*!
-             * Expands all unlabeled (i.e. independent) transitions of the given state and adds newly discovered states
-             * to the given queue.
-             *
-             * @param program The program that defines the transitions.
-             * @param stateInformation Provides information about the discovered state space.
-             * @param variableInformation A structure with information about the variables in the program.
-             * @param stateIndex The index of the state to expand.
-             * @param stateQueue The queue of states into which newly discovered states are inserted.
-             */
-            static void exploreUnlabeledTransitions(storm::ir::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue<uint_fast64_t>& stateQueue) {
-                StateType const* currentState = stateInformation.reachableStates[stateIndex];
-                
-                // Iterate over all modules.
-                for (uint_fast64_t i = 0; i < program.getNumberOfModules(); ++i) {
-                    storm::ir::Module const& module = program.getModule(i);
-
-                    // Iterate over all commands.
-                    for (uint_fast64_t j = 0; j < module.getNumberOfCommands(); ++j) {
-                        storm::ir::Command const& command = module.getCommand(j);
-                        
-                        // Only consider unlabeled commands.
-                        if (command.getActionName() != "") continue;
-                        // Skip the command, if it is not enabled.
-                        if (!command.getGuard()->getValueAsBool(currentState)) continue;
-                        
-                        ++stateInformation.numberOfChoices[stateIndex];
-                        
-                        std::set<uint_fast64_t> targetStates;
-                        
-                        // Iterate over all updates of the current command.
-                        for (uint_fast64_t k = 0; k < command.getNumberOfUpdates(); ++k) {
-                            storm::ir::Update const& update = command.getUpdate(k);
-                            
-                            // Obtain target state index.
-                            std::pair<bool, uint_fast64_t> flagAndTargetStateIndexPair = getOrAddStateIndex(applyUpdate(variableInformation, currentState, update), stateInformation);
-                            
-                            // If the state has not yet been discovered, add it to the queue.
-                            if (!flagAndTargetStateIndexPair.first) {
-                                stateQueue.push(flagAndTargetStateIndexPair.second);
-                            }
-                            
-                            // In any case, record that there is a transition to the newly found state to count the number
-                            // of transitions correctly.
-                            targetStates.insert(flagAndTargetStateIndexPair.second);
-                        }
-                        
-                        stateInformation.numberOfTransitions += targetStates.size();
-                    }
-                }
-            }
-
+    
             /*!
              * Retrieves all commands that are labeled with the given label and enabled in the given state, grouped by
              * modules.
@@ -684,149 +647,6 @@ namespace storm {
                 return result;
             }
             
-            /*!
-             * Expands all labeled (i.e. synchronizing) transitions of the given state and adds newly discovered states
-             * to the given queue.
-             *
-             * @param program The program that defines the transitions.
-             * @param stateInformation Provides information about the discovered state space.
-             * @param variableInformation A structure with information about the variables in the program.
-             * @param stateIndex The index of the state to expand.
-             * @param stateQueue The queue of states into which newly discovered states are inserted.
-             */
-            static void exploreLabeledTransitions(storm::ir::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue<uint_fast64_t>& stateQueue) {
-                for (std::string const& action : program.getActions()) {
-                    StateType const* currentState = stateInformation.reachableStates[stateIndex];
-                    boost::optional<std::vector<std::list<storm::ir::Command>>> optionalActiveCommandLists = getActiveCommandsByAction(program, currentState, action);
-                    
-                    // Only process this action label, if there is at least one feasible solution.
-                    if (optionalActiveCommandLists) {
-                        std::vector<std::list<storm::ir::Command>> const& activeCommandList = optionalActiveCommandLists.get();
-                        std::vector<std::list<storm::ir::Command>::const_iterator> iteratorList(activeCommandList.size());
-                        
-                        // Initialize the list of iterators.
-                        for (size_t i = 0; i < activeCommandList.size(); ++i) {
-                            iteratorList[i] = activeCommandList[i].cbegin();
-                        }
-                        
-                        // As long as there is one feasible combination of commands, keep on expanding it.
-                        bool done = false;
-                        while (!done) {
-                            ++stateInformation.numberOfChoices[stateIndex];
-                            
-                            std::unordered_set<StateType*, StateHash, StateCompare>* currentTargetStates = new std::unordered_set<StateType*, StateHash, StateCompare>();
-                            std::unordered_set<StateType*, StateHash, StateCompare>* newTargetStates = new std::unordered_set<StateType*, StateHash, StateCompare>();
-                            currentTargetStates->insert(new StateType(*currentState));
-                            
-                            // FIXME: This does not check whether a global variable is written multiple times. While the
-                            // behaviour for this is undefined anyway, a warning should be issued in that case.
-                            for (uint_fast64_t i = 0; i < iteratorList.size(); ++i) {
-                                storm::ir::Command const& command = *iteratorList[i];
-                                
-                                for (uint_fast64_t j = 0; j < command.getNumberOfUpdates(); ++j) {
-                                    storm::ir::Update const& update = command.getUpdate(j);
-                                    
-                                    for (auto const& targetState : *currentTargetStates) {
-                                        StateType* newTargetState = applyUpdate(variableInformation, targetState, currentState, update);
-                                        
-                                        auto existingState = newTargetStates->find(newTargetState);
-                                        if (existingState == newTargetStates->end()) {
-                                            newTargetStates->insert(newTargetState);
-                                        } else {
-                                            // If the state was already seen in one of the other updates, we need to delete this copy.
-                                            delete newTargetState;
-                                        }
-                                    }
-                                }
-                                
-                                // If there is one more command to come, shift the target states one time step back.
-                                if (i < iteratorList.size() - 1) {
-                                    for (auto const& state : *currentTargetStates) {
-                                        delete state;
-                                    }
-                                    delete currentTargetStates;
-                                    currentTargetStates = newTargetStates;
-                                    newTargetStates = new std::unordered_set<StateType*, StateHash, StateCompare>();
-                                }
-                            }
-                            
-                            // If there are states that have not yet been discovered, add them to the queue.
-                            for (auto const& state : *newTargetStates) {
-                                std::pair<bool, uint_fast64_t> flagAndNewStateIndexPair = getOrAddStateIndex(state, stateInformation);
-                                
-                                if (!flagAndNewStateIndexPair.first) {
-                                    stateQueue.push(flagAndNewStateIndexPair.second);
-                                }
-                            }
-                            
-                            stateInformation.numberOfTransitions += newTargetStates->size();
-                            
-                            // Dispose of the temporary maps.
-                            delete currentTargetStates;
-                            delete newTargetStates;
-                            
-                            // Now, check whether there is one more command combination to consider.
-                            bool movedIterator = false;
-                            for (int_fast64_t j = iteratorList.size() - 1; j >= 0; --j) {
-                                ++iteratorList[j];
-                                if (iteratorList[j] != activeCommandList[j].end()) {
-                                    movedIterator = true;
-                                } else {
-                                    // Reset the iterator to the beginning of the list.
-                                    iteratorList[j] = activeCommandList[j].begin();
-                                }
-                            }
-                            
-                            done = !movedIterator;
-                        }
-                    }
-                }
-            }
-            
-            /*!
-             * Explores the reachable state space given by the program and returns a list of all states as well as a
-             * mapping from states to their indices.
-             *
-             * @param program The program whose state space to explore.
-             * @param variableInformation A structure with information about the variables in the program.
-             * @return A structure containing all reachable states and a mapping from states to their indices.
-             */
-            static StateInformation exploreReachableStateSpace(storm::ir::Program const& program, VariableInformation const& variableInformation) {
-                StateInformation stateInformation;
-                
-                // Initialize a queue and insert the initial state.
-                std::queue<uint_fast64_t> stateQueue;
-                StateType* initialState = getInitialState(program, variableInformation);
-                getOrAddStateIndex(initialState, stateInformation);
-                stateQueue.push(stateInformation.stateToIndexMap[initialState]);
-                
-                // Now explore the current state until there is no more reachable state.
-                while (!stateQueue.empty()) {
-                    uint_fast64_t currentStateIndex = stateQueue.front();
-                    
-                    exploreUnlabeledTransitions(program, stateInformation, variableInformation, currentStateIndex, stateQueue);
-                    exploreLabeledTransitions(program, stateInformation, variableInformation, currentStateIndex, stateQueue);
-                    
-                    // If the state does not have any outgoing transitions, we fix this by adding a self-loop if the
-                    // option was set.
-                    if (stateInformation.numberOfChoices[currentStateIndex] == 0) {
-                        if (storm::settings::Settings::getInstance()->isSet("fixDeadlocks")) {
-                            ++stateInformation.numberOfTransitions;
-                            ++stateInformation.numberOfChoices[currentStateIndex] = 1;
-                        } else {
-                            LOG4CPLUS_ERROR(logger, "Error while creating sparse matrix from probabilistic program: found deadlock state. For fixing these, please provide the appropriate option.");
-                            throw storm::exceptions::WrongFormatException() << "Error while creating sparse matrix from probabilistic program: found deadlock state. For fixing these, please provide the appropriate option.";
-                        }
-                    }
-
-                    
-                    // After exploring a state, remove it from the queue.
-                    stateQueue.pop();
-                }
-                
-                return stateInformation;
-            }
-            
             /*!
              * Aggregates information about the variables in the program.
              *
@@ -879,7 +699,7 @@ namespace storm {
                 return result;
             }
             
-            static std::list<Choice<ValueType>> getUnlabeledTransitions(storm::ir::Program const& program, StateInformation const& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex) {
+            static std::list<Choice<ValueType>> getUnlabeledTransitions(storm::ir::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue<uint_fast64_t>& stateQueue) {
                 std::list<Choice<ValueType>> result;
                 
                 StateType const* currentState = stateInformation.reachableStates[stateIndex];
@@ -907,12 +727,17 @@ namespace storm {
                             storm::ir::Update const& update = command.getUpdate(k);
                             
                             // Obtain target state index.
-                            uint_fast64_t targetStateIndex = stateInformation.stateToIndexMap.at(applyUpdate(variableInformation, currentState, update));
+                            std::pair<bool, uint_fast64_t> flagTargetStateIndexPair = getOrAddStateIndex(applyUpdate(variableInformation, currentState, update), stateInformation);
+                            
+                            // If the state has not been discovered yet, add it to the queue of states to be explored.
+                            if (!flagTargetStateIndexPair.first) {
+                                stateQueue.push(flagTargetStateIndexPair.second);
+                            }
                             
                             // Update the choice by adding the probability/target state to it.
                             double probabilityToAdd = update.getLikelihoodExpression()->getValueAsDouble(currentState);
                             probabilitySum += probabilityToAdd;
-                            addProbabilityToChoice(choice, targetStateIndex, probabilityToAdd, {update.getGlobalIndex()});
+                            addProbabilityToChoice(choice, flagTargetStateIndexPair.second, probabilityToAdd, {update.getGlobalIndex()});
                         }
                         
                         // Check that the resulting distribution is in fact a distribution.
@@ -926,7 +751,7 @@ namespace storm {
                 return result;
             }
             
-            static std::list<Choice<ValueType>> getLabeledTransitions(storm::ir::Program const& program, StateInformation const& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex) {
+            static std::list<Choice<ValueType>> getLabeledTransitions(storm::ir::Program const& program, StateInformation& stateInformation, VariableInformation const& variableInformation, uint_fast64_t stateIndex, std::queue<uint_fast64_t>& stateQueue) {
                 std::list<Choice<ValueType>> result;
                 
                 for (std::string const& action : program.getActions()) {
@@ -946,9 +771,9 @@ namespace storm {
                         // As long as there is one feasible combination of commands, keep on expanding it.
                         bool done = false;
                         while (!done) {
-                            std::unordered_map<StateType*, double, StateHash, StateCompare>* currentTargetStates = new std::unordered_map<StateType*, double, StateHash, StateCompare>();
-                            std::unordered_map<StateType*, double, StateHash, StateCompare>* newTargetStates = new std::unordered_map<StateType*, double, StateHash, StateCompare>();
-                            (*currentTargetStates)[new StateType(*currentState)] = 1.0;
+                            std::unordered_map<StateType*, storm::storage::LabeledValues<double>, StateHash, StateCompare>* currentTargetStates = new std::unordered_map<StateType*, storm::storage::LabeledValues<double>, StateHash, StateCompare>();
+                            std::unordered_map<StateType*, storm::storage::LabeledValues<double>, StateHash, StateCompare>* newTargetStates = new std::unordered_map<StateType*, storm::storage::LabeledValues<double>, StateHash, StateCompare>();
+                            (*currentTargetStates)[new StateType(*currentState)] = storm::storage::LabeledValues<double>(1.0);
                             
                             // FIXME: This does not check whether a global variable is written multiple times. While the
                             // behaviour for this is undefined anyway, a warning should be issued in that case.
@@ -960,12 +785,24 @@ namespace storm {
                                     
                                     for (auto const& stateProbabilityPair : *currentTargetStates) {
                                         StateType* newTargetState = applyUpdate(variableInformation, stateProbabilityPair.first, currentState, update);
+
+                                        storm::storage::LabeledValues<double> newProbability;
+                                        
+                                        double updateProbability = update.getLikelihoodExpression()->getValueAsDouble(currentState);
+                                        for (auto const& valueLabelSetPair : stateProbabilityPair.second) {
+                                            // Copy the label set, so we can modify it.
+                                            std::set<uint_fast64_t> newLabelSet = valueLabelSetPair.second;
+                                            newLabelSet.insert(update.getGlobalIndex());
+                                            
+                                            newProbability.addValue(valueLabelSetPair.first * updateProbability, newLabelSet);
+                                        }
                                         
                                         auto existingStateProbabilityPair = newTargetStates->find(newTargetState);
                                         if (existingStateProbabilityPair == newTargetStates->end()) {
-                                            (*newTargetStates)[newTargetState] = stateProbabilityPair.second * update.getLikelihoodExpression()->getValueAsDouble(currentState);
-                                        } else {
-                                            existingStateProbabilityPair->second += stateProbabilityPair.second * update.getLikelihoodExpression()->getValueAsDouble(currentState);
+                                            (*newTargetStates)[newTargetState] = newProbability;
+                                        } else {                                            
+                                            existingStateProbabilityPair->second += newProbability;
+                                            
                                             // If the state was already seen in one of the other updates, we need to delete this copy.
                                             delete newTargetState;
                                         }
@@ -977,9 +814,10 @@ namespace storm {
                                     for (auto const& stateProbabilityPair : *currentTargetStates) {
                                         delete stateProbabilityPair.first;
                                     }
+                                    
                                     delete currentTargetStates;
                                     currentTargetStates = newTargetStates;
-                                    newTargetStates = new std::unordered_map<StateType*, double, StateHash, StateCompare>();
+                                    newTargetStates = new std::unordered_map<StateType*, storm::storage::LabeledValues<double>, StateHash, StateCompare>();
                                 }
                             }
                             
@@ -997,11 +835,18 @@ namespace storm {
                             }
                             
                             double probabilitySum = 0;
-                            std::set<uint_fast64_t> labels;
                             for (auto const& stateProbabilityPair : *newTargetStates) {
-                                uint_fast64_t newStateIndex = stateInformation.stateToIndexMap.at(stateProbabilityPair.first);
-                                addProbabilityToChoice(choice, newStateIndex, stateProbabilityPair.second, labels);
-                                probabilitySum += stateProbabilityPair.second;
+                                std::pair<bool, uint_fast64_t> flagTargetStateIndexPair = getOrAddStateIndex(stateProbabilityPair.first, stateInformation);
+
+                                // If the state has not been discovered yet, add it to the queue of states to be explored.
+                                if (!flagTargetStateIndexPair.first) {
+                                    stateQueue.push(flagTargetStateIndexPair.second);
+                                }
+                                
+                                for (auto const& probabilityLabelPair : stateProbabilityPair.second) {
+                                    addProbabilityToChoice(choice, flagTargetStateIndexPair.second, probabilityLabelPair.first, probabilityLabelPair.second);
+                                    probabilitySum += probabilityLabelPair.first;
+                                }
                             }
                             
                             // Check that the resulting distribution is in fact a distribution.
@@ -1035,18 +880,39 @@ namespace storm {
             }
             
             /*!
+             * Builds the transition matrix and the transition reward matrix based for the given program.
              *
+             * @param program The program for which to build the matrices.
+             * @param variableInformation A structure containing information about the variables in the program.
+             * @param transitionRewards A list of transition rewards that are to be considered in the transition reward
+             * matrix.
+             * @param stateInformation A structure containing information about the states of the program.
+             * @param deterministicModel A flag indicating whether the model is supposed to be deterministic or not.
+             * @param transitionMatrix A reference to an initialized matrix which is filled with all transitions by this
+             * function.
+             * @param transitionRewardMatrix A reference to an initialized matrix which is filled with all transition
+             * rewards by this function.
+             * @return A tuple containing a vector with all rows at which the nondeterministic choices of each state begin
+             * and a vector containing the labels associated with each choice.
              */
-            static std::pair<std::vector<uint_fast64_t>, std::vector<std::set<uint_fast64_t>>> buildMatrices(storm::ir::Program const& program, VariableInformation const& variableInformation, std::vector<storm::ir::TransitionReward> const& transitionRewards, StateInformation const& stateInformation, bool deterministicModel, storm::storage::SparseMatrix<ValueType>& transitionMatrix, storm::storage::SparseMatrix<ValueType>& transitionRewardMatrix) {
-                std::vector<uint_fast64_t> nondeterministicChoiceIndices(stateInformation.reachableStates.size() + 1);
+            static std::pair<std::vector<uint_fast64_t>, std::vector<std::set<uint_fast64_t>>> buildMatrices(storm::ir::Program const& program, VariableInformation const& variableInformation, std::vector<storm::ir::TransitionReward> const& transitionRewards, StateInformation& stateInformation, bool deterministicModel, storm::storage::SparseMatrix<ValueType>& transitionMatrix, storm::storage::SparseMatrix<ValueType>& transitionRewardMatrix) {
+                std::vector<uint_fast64_t> nondeterministicChoiceIndices;
                 std::vector<std::set<uint_fast64_t>> choiceLabels;
                 
-                // Add transitions and rewards for all states.
+                // Initialize a queue and insert the initial state.
+                std::queue<uint_fast64_t> stateQueue;
+                StateType* initialState = getInitialState(program, variableInformation);
+                getOrAddStateIndex(initialState, stateInformation);
+                stateQueue.push(stateInformation.stateToIndexMap[initialState]);
+                
+                // Now explore the current state until there is no more reachable state.
                 uint_fast64_t currentRow = 0;
-                for (uint_fast64_t currentState = 0; currentState < stateInformation.reachableStates.size(); ++currentState) {
+                while (!stateQueue.empty()) {
+                    uint_fast64_t currentState = stateQueue.front();
+            
                     // Retrieve all choices for the current state.
-                    std::list<Choice<ValueType>> allUnlabeledChoices = getUnlabeledTransitions(program, stateInformation, variableInformation, currentState);
-                    std::list<Choice<ValueType>> allLabeledChoices = getLabeledTransitions(program, stateInformation, variableInformation, currentState);
+                    std::list<Choice<ValueType>> allUnlabeledChoices = getUnlabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue);
+                    std::list<Choice<ValueType>> allLabeledChoices = getLabeledTransitions(program, stateInformation, variableInformation, currentState, stateQueue);
                     
                     uint_fast64_t totalNumberOfChoices = allUnlabeledChoices.size() + allLabeledChoices.size();
                     
@@ -1054,7 +920,10 @@ namespace storm {
                     // requested and issue an error otherwise.
                     if (totalNumberOfChoices == 0) {
                         if (storm::settings::Settings::getInstance()->isSet("fixDeadlocks")) {
-                            transitionMatrix.insertNextValue(currentRow, currentState, storm::utility::constGetOne<ValueType>());
+                            transitionMatrix.insertNextValue(currentRow, currentState, storm::utility::constGetOne<ValueType>(), true);
+                            if (transitionRewards.size() > 0) {
+                                transitionRewardMatrix.insertEmptyRow(true);
+                            }
                             ++currentRow;
                         } else {
                             LOG4CPLUS_ERROR(logger, "Error while creating sparse matrix from probabilistic program: found deadlock state. For fixing these, please provide the appropriate option.");
@@ -1065,6 +934,7 @@ namespace storm {
                         // or compose them to one choice.
                         if (deterministicModel) {
                             Choice<ValueType> globalChoice("");
+                            std::map<uint_fast64_t, ValueType> stateToRewardMap;
                             std::set<uint_fast64_t> allChoiceLabels;
                             
                             // Combine all the choices and scale them with the total number of choices of the current state.
@@ -1072,28 +942,51 @@ namespace storm {
                                 globalChoice.addChoiceLabels(choice.getChoiceLabels());
                                 for (auto const& stateProbabilityPair : choice) {
                                     globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                    
+                                    // Now add all rewards that match this choice.
+                                    for (auto const& transitionReward : transitionRewards) {
+                                        if (transitionReward.getActionName() == "" && transitionReward.getStatePredicate()->getValueAsBool(stateInformation.reachableStates.at(currentState))) {
+                                            stateToRewardMap[stateProbabilityPair.first] += ValueType(transitionReward.getRewardValue()->getValueAsDouble(stateInformation.reachableStates.at(currentState)));
+                                        }
+                                    }
                                 }
                             }
                             for (auto const& choice : allLabeledChoices) {
                                 globalChoice.addChoiceLabels(choice.getChoiceLabels());
                                 for (auto const& stateProbabilityPair : choice) {
                                     globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                                
+                                    // Now add all rewards that match this choice.
+                                    for (auto const& transitionReward : transitionRewards) {
+                                        if (transitionReward.getActionName() == choice.getActionLabel() && transitionReward.getStatePredicate()->getValueAsBool(stateInformation.reachableStates.at(currentState))) {
+                                            stateToRewardMap[stateProbabilityPair.first] += ValueType(transitionReward.getRewardValue()->getValueAsDouble(stateInformation.reachableStates.at(currentState)));
+                                        }
+                                    }
                                 }
                             }
 
                             
                             // Now add the resulting distribution as the only choice of the current state.
-                            nondeterministicChoiceIndices[currentState] = currentRow;
+                            nondeterministicChoiceIndices.push_back(currentRow);
                             choiceLabels.push_back(globalChoice.getChoiceLabels());
                             
                             for (auto const& stateProbabilityPair : globalChoice) {
-                                transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second);
+                                transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second, true);
+                            }
+                            
+                            // Add all transition rewards to the matrix and add dummy entry if there is none.
+                            if (stateToRewardMap.size() > 0) {
+                                for (auto const& stateRewardPair : stateToRewardMap) {
+                                    transitionRewardMatrix.insertNextValue(currentRow, stateRewardPair.first, stateRewardPair.second, true);
+                                }
+                            } else if (transitionRewards.size() > 0) {
+                                transitionRewardMatrix.insertEmptyRow(true);
                             }
                             
                             ++currentRow;
                         } else {
                             // If the model is nondeterministic, we add all choices individually.
-                            nondeterministicChoiceIndices[currentState] = currentRow;
+                            nondeterministicChoiceIndices.push_back(currentRow);
                             
                             // First, process all unlabeled choices.
                             for (auto const& choice : allUnlabeledChoices) {
@@ -1101,7 +994,7 @@ namespace storm {
                                 choiceLabels.emplace_back(std::move(choice.getChoiceLabels()));
                                 
                                 for (auto const& stateProbabilityPair : choice) {
-                                    transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second);
+                                    transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second, true);
                                     
                                     // Now add all rewards that match this choice.
                                     for (auto const& transitionReward : transitionRewards) {
@@ -1112,9 +1005,13 @@ namespace storm {
 
                                 }
                                 
-                                // Add all transition rewards to the matrix.
-                                for (auto const& stateRewardPair : stateToRewardMap) {
-                                    transitionRewardMatrix.insertNextValue(currentRow, stateRewardPair.first, stateRewardPair.second);
+                                // Add all transition rewards to the matrix and add dummy entry if there is none.
+                                if (stateToRewardMap.size() > 0) {
+                                    for (auto const& stateRewardPair : stateToRewardMap) {
+                                        transitionRewardMatrix.insertNextValue(currentRow, stateRewardPair.first, stateRewardPair.second, true);
+                                    }
+                                } else if (transitionRewards.size() > 0) {
+                                    transitionRewardMatrix.insertEmptyRow(true);
                                 }
                                 
                                 ++currentRow;
@@ -1126,29 +1023,35 @@ namespace storm {
                                 choiceLabels.emplace_back(std::move(choice.getChoiceLabels()));
 
                                 for (auto const& stateProbabilityPair : choice) {
-                                    transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second);
+                                    transitionMatrix.insertNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second, true);
                                     
                                     // Now add all rewards that match this choice.
                                     for (auto const& transitionReward : transitionRewards) {
-                                        if (transitionReward.getActionName() == "" && transitionReward.getStatePredicate()->getValueAsBool(stateInformation.reachableStates.at(currentState))) {
+                                        if (transitionReward.getActionName() == choice.getActionLabel() && transitionReward.getStatePredicate()->getValueAsBool(stateInformation.reachableStates.at(currentState))) {
                                             stateToRewardMap[stateProbabilityPair.first] += ValueType(transitionReward.getRewardValue()->getValueAsDouble(stateInformation.reachableStates.at(currentState)));
                                         }
                                     }
 
                                 }
                                 
-                                // Add all transition rewards to the matrix.
-                                for (auto const& stateRewardPair : stateToRewardMap) {
-                                    transitionRewardMatrix.insertNextValue(currentRow, stateRewardPair.first, stateRewardPair.second);
+                                // Add all transition rewards to the matrix and add dummy entry if there is none.
+                                if (stateToRewardMap.size() > 0) {
+                                    for (auto const& stateRewardPair : stateToRewardMap) {
+                                        transitionRewardMatrix.insertNextValue(currentRow, stateRewardPair.first, stateRewardPair.second, true);
+                                    }
+                                } else if (transitionRewards.size() > 0) {
+                                    transitionRewardMatrix.insertEmptyRow(true);
                                 }
 
                                 ++currentRow;
                             }
                         }
                     }
+                    
+                    stateQueue.pop();
                 }
                 
-                nondeterministicChoiceIndices[stateInformation.reachableStates.size()] = currentRow;
+                nondeterministicChoiceIndices.push_back(currentRow);
                 
                 return std::make_pair(nondeterministicChoiceIndices, choiceLabels);
             }
@@ -1166,38 +1069,28 @@ namespace storm {
                 
                 VariableInformation variableInformation = createVariableInformation(program);
                 
-                // Start by exploring the state space.
-                StateInformation stateInformation = exploreReachableStateSpace(program, variableInformation);
+                // Initialize the matrices.
+                modelComponents.transitionMatrix.initialize();
+                modelComponents.transitionRewardMatrix.initialize();
                 
-                // Then build the transition and reward matrices for the reachable states.
-                bool deterministicModel = program.getModelType() == storm::ir::Program::DTMC || program.getModelType() == storm::ir::Program::CTMC;
-                if (deterministicModel) {
-                    modelComponents.transitionMatrix = storm::storage::SparseMatrix<ValueType>(stateInformation.reachableStates.size());
-                    modelComponents.transitionMatrix.initialize();
-                    modelComponents.transitionRewardMatrix = storm::storage::SparseMatrix<ValueType>(stateInformation.reachableStates.size());
-                    modelComponents.transitionRewardMatrix.initialize();
-                } else {
-                    uint_fast64_t totalNumberOfChoices = 0;
-                    for (auto numberOfChoices : stateInformation.numberOfChoices) {
-                        totalNumberOfChoices += numberOfChoices;
-                    }
-                    modelComponents.transitionMatrix = storm::storage::SparseMatrix<ValueType>(totalNumberOfChoices, stateInformation.reachableStates.size());
-                    modelComponents.transitionMatrix.initialize();
-                    modelComponents.transitionRewardMatrix = storm::storage::SparseMatrix<ValueType>(totalNumberOfChoices);
-                    modelComponents.transitionRewardMatrix.initialize();
-                }
+                // Create the structure for storing the reachable state space.
+                StateInformation stateInformation;
                 
                 // Get the selected reward model or create an empty one if none is selected.
                 storm::ir::RewardModel const& rewardModel = rewardModelName != "" ? program.getRewardModel(rewardModelName) : storm::ir::RewardModel();
                 
+                // Determine whether we have to combine different choices to one or whether this model can have more than
+                // one choice per state.
+                bool deterministicModel = program.getModelType() == storm::ir::Program::DTMC || program.getModelType() == storm::ir::Program::CTMC;
+
                 // Build the transition and reward matrices.
                 std::pair<std::vector<uint_fast64_t>, std::vector<std::set<uint_fast64_t>>> nondeterministicChoiceIndicesAndChoiceLabelsPair = buildMatrices(program, variableInformation, rewardModel.getTransitionRewards(), stateInformation, deterministicModel, modelComponents.transitionMatrix, modelComponents.transitionRewardMatrix);
                 modelComponents.nondeterministicChoiceIndices = std::move(nondeterministicChoiceIndicesAndChoiceLabelsPair.first);
                 modelComponents.choiceLabeling = std::move(nondeterministicChoiceIndicesAndChoiceLabelsPair.second);
                 
                 // Finalize the resulting matrices.
-                modelComponents.transitionMatrix.finalize();
-                modelComponents.transitionRewardMatrix.finalize();
+                modelComponents.transitionMatrix.finalize(true);
+                modelComponents.transitionRewardMatrix.finalize(true);
                 
                 // Now build the state labeling.
                 modelComponents.stateLabeling = buildStateLabeling(program, variableInformation, stateInformation);
@@ -1213,6 +1106,14 @@ namespace storm {
                 return modelComponents;
             }
             
+            /*!
+             * Builds the state labeling for the given program.
+             *
+             * @param program The program for which to build the state labeling.
+             * @param variableInformation Information about the variables in the program.
+             * @param stateInformation Information about the state space of the program.
+             * @return The state labeling of the given program.
+             */
             static storm::models::AtomicPropositionsLabeling buildStateLabeling(storm::ir::Program const& program, VariableInformation const& variableInformation, StateInformation const& stateInformation) {
                 std::map<std::string, std::shared_ptr<storm::ir::expressions::BaseExpression>> const& labels = program.getLabels();
                 
@@ -1241,6 +1142,13 @@ namespace storm {
                 return result;
             }
 
+            /*!
+             * Builds the state rewards for the given state space.
+             *
+             * @param rewards A vector of state rewards to consider.
+             * @param stateInformation Information about the state space.
+             * @return A vector containing the state rewards for the state space.
+             */
             static std::vector<ValueType> buildStateRewards(std::vector<storm::ir::StateReward> const& rewards, StateInformation const& stateInformation) {
                 std::vector<ValueType> result(stateInformation.reachableStates.size());
                 for (uint_fast64_t index = 0; index < stateInformation.reachableStates.size(); index++) {
diff --git a/src/counterexamples/MinimalLabelSetGenerator.h b/src/counterexamples/MILPMinimalLabelSetGenerator.h
similarity index 99%
rename from src/counterexamples/MinimalLabelSetGenerator.h
rename to src/counterexamples/MILPMinimalLabelSetGenerator.h
index 9f77a605d..611da0f7e 100644
--- a/src/counterexamples/MinimalLabelSetGenerator.h
+++ b/src/counterexamples/MILPMinimalLabelSetGenerator.h
@@ -1,14 +1,14 @@
 /*
- * MinimalLabelSetGenerator.h
+ * MILPMinimalLabelSetGenerator.h
  *
  *  Created on: 15.09.2013
  *      Author: Christian Dehnert
  */
 
-#ifndef STORM_COUNTEREXAMPLES_MINIMALCOMMANDSETGENERATOR_MDP_H_
-#define STORM_COUNTEREXAMPLES_MINIMALCOMMANDSETGENERATOR_MDP_H_
+#ifndef STORM_COUNTEREXAMPLES_MILPMINIMALLABELSETGENERATOR_MDP_H_
+#define STORM_COUNTEREXAMPLES_MILPMINIMALLABELSETGENERATOR_MDP_H_
 
-// To detect whether the usage of Gurobi is possible, this include is neccessary
+// To detect whether the usage of Gurobi is possible, this include is neccessary.
 #include "storm-config.h"
 
 #ifdef STORM_HAVE_GUROBI
@@ -386,7 +386,7 @@ namespace storm {
                     variableNameBuffer.str("");
                     variableNameBuffer.clear();
                     variableNameBuffer << "p" << state;
-                    error = GRBaddvar(model, 0, nullptr, nullptr, maximizeProbability ? (labeledMdp.getInitialStates().get(state) ? -0.5 : 0) : 0, 0.0, 1.0, GRB_CONTINUOUS, variableNameBuffer.str().c_str());
+                    error = GRBaddvar(model, 0, nullptr, nullptr, 0.0, 0.0, 1.0, GRB_CONTINUOUS, variableNameBuffer.str().c_str());
                     if (error) {
                         LOG4CPLUS_ERROR(logger, "Could not create Gurobi variable (" << GRBgeterrormsg(env) << ").");
                         throw storm::exceptions::InvalidStateException() << "Could not create Gurobi variable (" << GRBgeterrormsg(env) << ").";
@@ -1258,7 +1258,7 @@ namespace storm {
             
         };
         
-    }
-}
+    } // namespace counterexamples
+} // namespace storm
 
-#endif /* STORM_COUNTEREXAMPLES_MINIMALLABELSETGENERATOR_MDP_H_ */
+#endif /* STORM_COUNTEREXAMPLES_MILPMINIMALLABELSETGENERATOR_MDP_H_ */
diff --git a/src/counterexamples/SMTMinimalCommandSetGenerator.h b/src/counterexamples/SMTMinimalCommandSetGenerator.h
new file mode 100644
index 000000000..4608f4e8e
--- /dev/null
+++ b/src/counterexamples/SMTMinimalCommandSetGenerator.h
@@ -0,0 +1,50 @@
+/*
+ * SMTMinimalCommandSetGenerator.h
+ *
+ *  Created on: 01.10.2013
+ *      Author: Christian Dehnert
+ */
+
+#ifndef STORM_COUNTEREXAMPLES_SMTMINIMALCOMMANDSETGENERATOR_MDP_H_
+#define STORM_COUNTEREXAMPLES_SMTMINIMALCOMMANDSETGENERATOR_MDP_H_
+
+// To detect whether the usage of Z3 is possible, this include is neccessary.
+#include "storm-config.h"
+
+namespace storm {
+    namespace counterexamples {
+        
+        /*!
+         * This class provides functionality to generate a minimal counterexample to a probabilistic reachability
+         * property in terms of used labels.
+         */
+        template <class T>
+        class MinimalLabelSetGenerator {
+#ifdef STORM_HAVE_Z3
+        private:
+            
+#endif
+            
+        public:
+            static std::unordered_set<uint_fast64_t> getMinimalCommandSet(storm::ir::Program const& program, storm::models::Mdp<T> const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, double probabilityThreshold, bool checkThresholdFeasible = false) {
+#ifdef STORM_HAVE_Z3
+                // (0) Check whether the MDP is indeed labeled.
+                if (!labeledMdp.hasChoiceLabels()) {
+                    throw storm::exceptions::InvalidArgumentException() << "Minimal command set generation is impossible for unlabeled model.";
+                }
+                
+                
+                
+                return std::unordered_set<uint_fast64_t>();
+                
+#else
+                throw storm::exceptions::NotImplementedException() << "This functionality is unavailable since StoRM has been compiled without support for Z3.";
+#endif
+            }
+            
+        }
+        
+    } // namespace counterexamples
+} // namespace storm
+
+#endif /* STORM_COUNTEREXAMPLES_SMTMINIMALCOMMANDSETGENERATOR_MDP_H_ */
diff --git a/src/models/AbstractNondeterministicModel.h b/src/models/AbstractNondeterministicModel.h
index bb6b0fd37..c4569473f 100644
--- a/src/models/AbstractNondeterministicModel.h
+++ b/src/models/AbstractNondeterministicModel.h
@@ -6,246 +6,241 @@
 #include <memory>
 
 namespace storm {
-
-namespace models {
-
-/*!
- *	@brief	Base class for all non-deterministic model classes.
- *
- *	This is base class defines a common interface for all non-deterministic models.
- */
-template<class T>
-class AbstractNondeterministicModel: public AbstractModel<T> {
-
-	public:
-		/*! Constructs an abstract non-determinstic model from the given parameters.
-		 * All values are copied.
-		 * @param transitionMatrix The matrix representing the transitions in the model.
-		 * @param stateLabeling The labeling that assigns a set of atomic
-		 * propositions to each state.
-		 * @param choiceIndices A mapping from states to rows in the transition matrix.
-		 * @param stateRewardVector The reward values associated with the states.
-		 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
-         * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
-		 */
-		AbstractNondeterministicModel(
-			storm::storage::SparseMatrix<T> const& transitionMatrix, 
-			storm::models::AtomicPropositionsLabeling const& stateLabeling,
-			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices,
-			boost::optional<std::vector<T>> const& optionalStateRewardVector, 
-			boost::optional<storm::storage::SparseMatrix<T>> const& optionalTransitionRewardMatrix,
-            boost::optional<std::vector<std::set<uint_fast64_t>>> const& optionalChoiceLabeling)
+    namespace models {
+        
+        /*!
+         *	@brief	Base class for all non-deterministic model classes.
+         *
+         *	This is base class defines a common interface for all non-deterministic models.
+         */
+        template<class T>
+        class AbstractNondeterministicModel: public AbstractModel<T> {
+            
+        public:
+            /*! Constructs an abstract non-determinstic model from the given parameters.
+             * All values are copied.
+             * @param transitionMatrix The matrix representing the transitions in the model.
+             * @param stateLabeling The labeling that assigns a set of atomic
+             * propositions to each state.
+             * @param choiceIndices A mapping from states to rows in the transition matrix.
+             * @param stateRewardVector The reward values associated with the states.
+             * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+             * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state.
+             */
+            AbstractNondeterministicModel(storm::storage::SparseMatrix<T> const& transitionMatrix,
+                                          storm::models::AtomicPropositionsLabeling const& stateLabeling,
+                                          std::vector<uint_fast64_t> const& nondeterministicChoiceIndices,
+                                          boost::optional<std::vector<T>> const& optionalStateRewardVector,
+                                          boost::optional<storm::storage::SparseMatrix<T>> const& optionalTransitionRewardMatrix,
+                                          boost::optional<std::vector<std::set<uint_fast64_t>>> const& optionalChoiceLabeling)
 			: AbstractModel<T>(transitionMatrix, stateLabeling, optionalStateRewardVector, optionalTransitionRewardMatrix, optionalChoiceLabeling) {
 				this->nondeterministicChoiceIndices = nondeterministicChoiceIndices;
-		}
-
-		/*! Constructs an abstract non-determinstic model from the given parameters.
-		 * All values are moved.
-		 * @param transitionMatrix The matrix representing the transitions in the model.
-		 * @param stateLabeling The labeling that assigns a set of atomic
-		 * propositions to each state.
-		 * @param choiceIndices A mapping from states to rows in the transition matrix.
-		 * @param stateRewardVector The reward values associated with the states.
-		 * @param transitionRewardMatrix The reward values associated with the transitions of the model.
-		 */
-		AbstractNondeterministicModel(
-			storm::storage::SparseMatrix<T>&& transitionMatrix, 
-			storm::models::AtomicPropositionsLabeling&& stateLabeling,
-			std::vector<uint_fast64_t>&& nondeterministicChoiceIndices,
-			boost::optional<std::vector<T>>&& optionalStateRewardVector, 
-			boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,
-            boost::optional<std::vector<std::set<uint_fast64_t>>>&& optionalChoiceLabeling)
-			// The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class
-			: AbstractModel<T>(std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix),
+            }
+            
+            /*! Constructs an abstract non-determinstic model from the given parameters.
+             * All values are moved.
+             * @param transitionMatrix The matrix representing the transitions in the model.
+             * @param stateLabeling The labeling that assigns a set of atomic
+             * propositions to each state.
+             * @param choiceIndices A mapping from states to rows in the transition matrix.
+             * @param stateRewardVector The reward values associated with the states.
+             * @param transitionRewardMatrix The reward values associated with the transitions of the model.
+             */
+            AbstractNondeterministicModel(storm::storage::SparseMatrix<T>&& transitionMatrix,
+                                          storm::models::AtomicPropositionsLabeling&& stateLabeling,
+                                          std::vector<uint_fast64_t>&& nondeterministicChoiceIndices,
+                                          boost::optional<std::vector<T>>&& optionalStateRewardVector,
+                                          boost::optional<storm::storage::SparseMatrix<T>>&& optionalTransitionRewardMatrix,
+                                          boost::optional<std::vector<std::set<uint_fast64_t>>>&& optionalChoiceLabeling)
+                // The std::move call must be repeated here because otherwise this calls the copy constructor of the Base Class
+                : AbstractModel<T>(std::move(transitionMatrix), std::move(stateLabeling), std::move(optionalStateRewardVector), std::move(optionalTransitionRewardMatrix),
                                std::move(optionalChoiceLabeling)), nondeterministicChoiceIndices(std::move(nondeterministicChoiceIndices)) {
-			// Intentionally left empty.	
-		}
-
-		/*!
-		 * Destructor.
-		 */
-		virtual ~AbstractNondeterministicModel() {
-			// Intentionally left empty.
-		}
-
-		/*!
-		 * Copy Constructor.
-		 */
-		AbstractNondeterministicModel(AbstractNondeterministicModel const& other) : AbstractModel<T>(other),
-				nondeterministicChoiceIndices(other.nondeterministicChoiceIndices) {
-			// Intentionally left empty.
-		}
-
-		/*!
-		 * Move Constructor.
-		 */
-		AbstractNondeterministicModel(AbstractNondeterministicModel&& other) : AbstractModel<T>(std::move(other)),
-				nondeterministicChoiceIndices(std::move(other.nondeterministicChoiceIndices)) {
-			// Intentionally left empty.
-		}
-
-		/*!
-		 * Returns the number of choices for all states of the MDP.
-		 * @return The number of choices for all states of the MDP.
-		 */
-		uint_fast64_t getNumberOfChoices() const {
-			return this->transitionMatrix.getRowCount();
-		}
-    
-		/*!
-		 * Retrieves the size of the internal representation of the model in memory.
-		 * @return the size of the internal representation of the model in memory
-		 * measured in bytes.
-		 */
-		virtual uint_fast64_t getSizeInMemory() const {
-			return AbstractModel<T>::getSizeInMemory() + nondeterministicChoiceIndices.size() * sizeof(uint_fast64_t);
-		}
-
-		/*!
-		 * Retrieves the vector indicating which matrix rows represent non-deterministic choices
-		 * of a certain state.
-		 * @param the vector indicating which matrix rows represent non-deterministic choices
-		 * of a certain state.
-		 */
-		std::vector<uint_fast64_t> const& getNondeterministicChoiceIndices() const {
-			return nondeterministicChoiceIndices;
-		}
-    
-        virtual typename storm::storage::SparseMatrix<T>::Rows getRows(uint_fast64_t state) const override {
-            return this->transitionMatrix.getRows(nondeterministicChoiceIndices[state], nondeterministicChoiceIndices[state + 1] - 1);
-        }
-    
-        virtual typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIteratorBegin(uint_fast64_t state) const override {
-            return this->transitionMatrix.begin(nondeterministicChoiceIndices[state]);
-        }
+                // Intentionally left empty.
+            }
+            
+            /*!
+             * Destructor.
+             */
+            virtual ~AbstractNondeterministicModel() {
+                // Intentionally left empty.
+            }
+            
+            /*!
+             * Copy Constructor.
+             */
+            AbstractNondeterministicModel(AbstractNondeterministicModel const& other) : AbstractModel<T>(other),
+                nondeterministicChoiceIndices(other.nondeterministicChoiceIndices) {
+                // Intentionally left empty.
+            }
+            
+            /*!
+             * Move Constructor.
+             */
+            AbstractNondeterministicModel(AbstractNondeterministicModel&& other) : AbstractModel<T>(std::move(other)),
+                nondeterministicChoiceIndices(std::move(other.nondeterministicChoiceIndices)) {
+                // Intentionally left empty.
+            }
+            
+            /*!
+             * Returns the number of choices for all states of the MDP.
+             * @return The number of choices for all states of the MDP.
+             */
+            uint_fast64_t getNumberOfChoices() const {
+                return this->transitionMatrix.getRowCount();
+            }
+            
+            /*!
+             * Retrieves the size of the internal representation of the model in memory.
+             *
+             * @return the size of the internal representation of the model in memory
+             * measured in bytes.
+             */
+            virtual uint_fast64_t getSizeInMemory() const {
+                return AbstractModel<T>::getSizeInMemory() + nondeterministicChoiceIndices.size() * sizeof(uint_fast64_t);
+            }
+            
+            /*!
+             * Retrieves the vector indicating which matrix rows represent non-deterministic choices
+             * of a certain state.
+             * @param the vector indicating which matrix rows represent non-deterministic choices
+             * of a certain state.
+             */
+            std::vector<uint_fast64_t> const& getNondeterministicChoiceIndices() const {
+                return nondeterministicChoiceIndices;
+            }
+            
+            virtual typename storm::storage::SparseMatrix<T>::Rows getRows(uint_fast64_t state) const override {
+                return this->transitionMatrix.getRows(nondeterministicChoiceIndices[state], nondeterministicChoiceIndices[state + 1] - 1);
+            }
+        
+            virtual typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIteratorBegin(uint_fast64_t state) const override {
+                return this->transitionMatrix.begin(nondeterministicChoiceIndices[state]);
+            }
     
-        virtual typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIteratorEnd(uint_fast64_t state) const override {
+            virtual typename storm::storage::SparseMatrix<T>::ConstRowIterator rowIteratorEnd(uint_fast64_t state) const override {
             return this->transitionMatrix.end(nondeterministicChoiceIndices[state + 1] - 1);
-        }
+            }
 
-		/*!
-		 * Calculates a hash over all values contained in this Model.
-		 * @return size_t A Hash Value
-		 */
-		virtual size_t getHash() const override {
-			std::size_t result = 0;
-			std::size_t hashTmp = storm::utility::Hash<uint_fast64_t>::getHash(nondeterministicChoiceIndices);
-			std::cout << "nondeterministicChoiceIndices Hash: " << hashTmp << std::endl;
+            /*!
+             * Calculates a hash over all values contained in this Model.
+             * @return size_t A Hash Value
+             */
+            virtual size_t getHash() const override {
+                std::size_t result = 0;
+                std::size_t hashTmp = storm::utility::Hash<uint_fast64_t>::getHash(nondeterministicChoiceIndices);
+                boost::hash_combine(result, AbstractModel<T>::getHash());
+                boost::hash_combine(result, hashTmp);
+                return result;
+            }
 
-			boost::hash_combine(result, AbstractModel<T>::getHash());
-			boost::hash_combine(result, hashTmp);
-			return result;
-		}
+            /*!
+            * Prints information about the model to the specified stream.
+            * @param out The stream the information is to be printed to.
+            */
+            virtual void printModelInformationToStream(std::ostream& out) const override {
+                out << "-------------------------------------------------------------- " << std::endl;
+                out << "Model type: \t\t" << this->getType() << std::endl;
+                out << "States: \t\t" << this->getNumberOfStates() << std::endl;
+                out << "Transitions: \t\t" << this->getNumberOfTransitions() << std::endl;
+                out << "Choices: \t\t" << this->getNumberOfChoices() << std::endl;
+                this->getStateLabeling().printAtomicPropositionsInformationToStream(out);
+                out << "Size in memory: \t" << (this->getSizeInMemory())/1024 << " kbytes" << std::endl;
+                out << "-------------------------------------------------------------- " << std::endl;
+            }
 
-        /*!
-         * Prints information about the model to the specified stream.
-         * @param out The stream the information is to be printed to.
-         */
-        virtual void printModelInformationToStream(std::ostream& out) const override {
-            out << "-------------------------------------------------------------- " << std::endl;
-            out << "Model type: \t\t" << this->getType() << std::endl;
-            out << "States: \t\t" << this->getNumberOfStates() << std::endl;
-            out << "Transitions: \t\t" << this->getNumberOfTransitions() << std::endl;
-            out << "Choices: \t\t" << this->getNumberOfChoices() << std::endl;
-            this->getStateLabeling().printAtomicPropositionsInformationToStream(out);
-            out << "Size in memory: \t" << (this->getSizeInMemory())/1024 << " kbytes" << std::endl;
-            out << "-------------------------------------------------------------- " << std::endl;
-        }
+            virtual void writeDotToStream(std::ostream& outStream, bool includeLabeling = true, storm::storage::BitVector const* subsystem = nullptr, std::vector<T> const* firstValue = nullptr, std::vector<T> const* secondValue = nullptr, std::vector<uint_fast64_t> const* stateColoring = nullptr, std::vector<std::string> const* colors = nullptr, std::vector<uint_fast64_t>* scheduler = nullptr, bool finalizeOutput = true) const override {
+                AbstractModel<T>::writeDotToStream(outStream, includeLabeling, subsystem, firstValue, secondValue, stateColoring, colors, scheduler, false);
     
-        virtual void writeDotToStream(std::ostream& outStream, bool includeLabeling = true, storm::storage::BitVector const* subsystem = nullptr, std::vector<T> const* firstValue = nullptr, std::vector<T> const* secondValue = nullptr, std::vector<uint_fast64_t> const* stateColoring = nullptr, std::vector<std::string> const* colors = nullptr, std::vector<uint_fast64_t>* scheduler = nullptr, bool finalizeOutput = true) const override {
-            AbstractModel<T>::writeDotToStream(outStream, includeLabeling, subsystem, firstValue, secondValue, stateColoring, colors, scheduler, false);
-
-            // Write the probability distributions for all the states.
-            auto rowIt = this->transitionMatrix.begin();
-            for (uint_fast64_t state = 0, highestStateIndex = this->getNumberOfStates() - 1; state <= highestStateIndex; ++state) {
-                uint_fast64_t rowCount = nondeterministicChoiceIndices[state + 1] - nondeterministicChoiceIndices[state];
-                bool highlightChoice = true;
-                
-                // For this, we need to iterate over all available nondeterministic choices in the current state.
-                for (uint_fast64_t row = 0; row < rowCount; ++row, ++rowIt) {
-                    if (scheduler != nullptr) {
-                        // If the scheduler picked the current choice, we will not make it dotted, but highlight it.
-                        if ((*scheduler)[state] == row) {
-                            highlightChoice = true;
-                        } else {
-                            highlightChoice = false;
-                        }
-                    }
-                    
-                    // For each nondeterministic choice, we draw an arrow to an intermediate node to better display
-                    // the grouping of transitions.
-                    outStream << "\t\"" << state << "c" << row << "\" [shape = \"point\"";
-                    
-                    // If we were given a scheduler to highlight, we do so now.
-                    if (scheduler != nullptr) {
-                        if (highlightChoice) {
-                            outStream << ", fillcolor=\"red\"";
+                // Write the probability distributions for all the states.
+                auto rowIt = this->transitionMatrix.begin();
+                for (uint_fast64_t state = 0, highestStateIndex = this->getNumberOfStates() - 1; state <= highestStateIndex; ++state) {
+                    uint_fast64_t rowCount = nondeterministicChoiceIndices[state + 1] - nondeterministicChoiceIndices[state];
+                    bool highlightChoice = true;
+        
+                    // For this, we need to iterate over all available nondeterministic choices in the current state.
+                    for (uint_fast64_t row = 0; row < rowCount; ++row, ++rowIt) {
+                        if (scheduler != nullptr) {
+                            // If the scheduler picked the current choice, we will not make it dotted, but highlight it.
+                            if ((*scheduler)[state] == row) {
+                                highlightChoice = true;
+                            } else {
+                                highlightChoice = false;
+                            }
                         }
-                    }
-                    outStream << "];" << std::endl;
-                    
-                    outStream << "\t" << state << " -> \"" << state << "c" << row << "\"";
-                    
-                    // 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\"]";
+            
+                        // For each nondeterministic choice, we draw an arrow to an intermediate node to better display
+                        // the grouping of transitions.
+                        outStream << "\t\"" << state << "c" << row << "\" [shape = \"point\"";
+            
+                        // If we were given a scheduler to highlight, we do so now.
+                        if (scheduler != nullptr) {
+                            if (highlightChoice) {
+                                outStream << ", fillcolor=\"red\"";
+                            }
                         }
-                    }
-                    outStream << ";" << std::endl;
-                    
-                    // Now draw all probabilitic arcs that belong to this nondeterminstic choice.
-                    for (auto transitionIt = rowIt.begin(), transitionIte = rowIt.end(); transitionIt != transitionIte; ++transitionIt) {
-                        if (subsystem == nullptr || subsystem->get(transitionIt.column())) {
-                            outStream << "\t\"" << state << "c" << row << "\" -> " << transitionIt.column() << " [ label= \"" << transitionIt.value() << "\" ]";
+                        outStream << "];" << std::endl;
                         
-                            // 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 << "\t" << state << " -> \"" << state << "c" << row << "\"";
+            
+                        // 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;
+            
+                        // Now draw all probabilitic arcs that belong to this nondeterminstic choice.
+                        for (auto transitionIt = rowIt.begin(), transitionIte = rowIt.end(); transitionIt != transitionIte; ++transitionIt) {
+                            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\"]";
+                                    }
                                 }
+                                outStream << ";" << std::endl;
                             }
-                            outStream << ";" << std::endl;
                         }
                     }
                 }
+    
+                if (finalizeOutput) {
+                    outStream << "}" << std::endl;
+                }
             }
-        
-            if (finalizeOutput) {
-                outStream << "}" << std::endl;
-            }
-        }
-		
-		/*!
-		 * Assigns this model a new set of choiceLabels, giving each choice a label with the stateId
-		 * @return void
-		 */
-		virtual void setStateIdBasedChoiceLabeling() override {
-			std::vector<std::set<uint_fast64_t>> newChoiceLabeling;
-
-			size_t stateCount = this->getNumberOfStates();
-			size_t choiceCount = this->getNumberOfChoices();
-			newChoiceLabeling.resize(choiceCount);
-
-			for (size_t state = 0; state < stateCount; ++state) {
-				for (size_t choice = this->nondeterministicChoiceIndices.at(state); choice < this->nondeterministicChoiceIndices.at(state + 1); ++choice) {
-					newChoiceLabeling.at(choice).insert(state);
-				}
-			}
-
-			this->choiceLabeling.reset(newChoiceLabeling);
-		}
-	private:
-		/*! A vector of indices mapping states to the choices (rows) in the transition matrix. */
-		std::vector<uint_fast64_t> nondeterministicChoiceIndices;
-};
-
-} // namespace models
 
+            /*!
+             * Assigns this model a new set of choiceLabels, giving each choice a label with the stateId
+             * @return void
+             */
+            virtual void setStateIdBasedChoiceLabeling() override {
+                std::vector<std::set<uint_fast64_t>> newChoiceLabeling;
+    
+                size_t stateCount = this->getNumberOfStates();
+                size_t choiceCount = this->getNumberOfChoices();
+                newChoiceLabeling.resize(choiceCount);
+    
+                for (size_t state = 0; state < stateCount; ++state) {
+                    for (size_t choice = this->nondeterministicChoiceIndices.at(state); choice < this->nondeterministicChoiceIndices.at(state + 1); ++choice) {
+                        newChoiceLabeling.at(choice).insert(state);
+                    }
+                }
+    
+                this->choiceLabeling.reset(newChoiceLabeling);
+            }
+                        
+        private:
+            /*! A vector of indices mapping states to the choices (rows) in the transition matrix. */
+            std::vector<uint_fast64_t> nondeterministicChoiceIndices;
+        };
+    } // namespace models
 } // namespace storm
 
 #endif /* STORM_MODELS_ABSTRACTDETERMINISTICMODEL_H_ */
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index 849a42344..7b012803d 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -50,7 +50,7 @@ public:
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
 		if (this->hasTransitionRewards()) {
-			if (!this->getTransitionRewardMatrix().containsAllPositionsOf(this->getTransitionMatrix())) {
+			if (!this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix())) {
 				LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
 				throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
 			}
@@ -78,7 +78,7 @@ public:
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
 		if (this->hasTransitionRewards()) {
-			if (!this->getTransitionRewardMatrix().containsAllPositionsOf(this->getTransitionMatrix())) {
+			if (!this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix())) {
 				LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
 				throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
 			}
diff --git a/src/models/Mdp.h b/src/models/Mdp.h
index cb988c90b..0f97e304a 100644
--- a/src/models/Mdp.h
+++ b/src/models/Mdp.h
@@ -17,6 +17,7 @@
 #include "src/storage/SparseMatrix.h"
 #include "src/settings/Settings.h"
 #include "src/models/AbstractNondeterministicModel.h"
+#include "src/utility/set.h"
 
 namespace storm {
 
@@ -54,6 +55,12 @@ public:
 			LOG4CPLUS_ERROR(logger, "Probability matrix is invalid.");
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
+        if (this->hasTransitionRewards()) {
+            if (!this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix())) {
+                LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
+                throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
+            }
+        }
 	}
 
 	/*!
@@ -78,6 +85,12 @@ public:
 			LOG4CPLUS_ERROR(logger, "Probability matrix is invalid.");
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
+        if (this->hasTransitionRewards()) {
+            if (!this->getTransitionRewardMatrix().isSubmatrixOf(this->getTransitionMatrix())) {
+                LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
+                throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
+            }
+        }
 	}
 
 	/*!
@@ -113,6 +126,47 @@ public:
 		return MDP;
 	}
 
+    /*!
+     * Constructs an MDP by copying the given MDP and restricting the choices of each state to the ones whose label set
+     * is contained in the given label set.
+     *
+     * @param originalModel The model to restrict.
+     * @param enabledChoiceLabels A set of labels that determines which choices of the original model can be taken
+     * and which ones need to be ignored.
+     */
+    Mdp<T> restrictChoiceLabels(Mdp<T> const& originalModel, std::set<uint_fast64_t> const& enabledChoiceLabels) {
+        // Only perform this operation if the given model has choice labels.
+        if (!originalModel.hasChoiceLabels()) {
+            throw storm::exceptions::InvalidArgumentException() << "Restriction to label set is impossible for unlabeled model.";
+        }
+        
+        std::vector<std::set<uint_fast64_t>> const& choiceLabeling = this->getChoiceLabeling();
+        
+        storm::storage::SparseMatrix<T> transitionMatrix;
+        transitionMatrix.initialize();
+        
+        // Check for each choice of each state, whether the choice labels are fully contained in the given label set.
+        for(uint_fast64_t state = 0; state < this->getNumberOfStates(); ++state) {
+            for (uint_fast64_t choice = this->getNondeterministicChoiceIndices()[state]; choice < this->getNondeterministicChoiceIndices()[state + 1]; ++choice) {
+                bool choiceValid = storm::utility::set::isSubsetOf(choiceLabeling[state], enabledChoiceLabels);
+                
+                // If the choice is valid, copy over all its elements.
+                if (choiceValid) {
+                    typename storm::storage::SparseMatrix<T>::Rows row = this->getTransitionMatrix().getRows(choice, choice);
+                    for (typename storm::storage::SparseMatrix<T>::ConstIterator rowIt = row.begin(), rowIte = row.end(); rowIt != rowIte; ++rowIt) {
+                        transitionMatrix.insertNextValue(choice, rowIt.column(), rowIt.value(), true);
+                    }
+                } else {
+                    // If the choice may not be taken, we insert a self-loop to the state instead.
+                    transitionMatrix.insertNextValue(choice, state, storm::utility::constGetOne<T>(), true);
+                }
+            }
+        }
+        
+        Mdp<T> restrictedMdp(std::move(transitionMatrix), storm::models::AtomicPropositionsLabeling(this->getStateLabeling()), std::vector<uint_fast64_t>(this->getNondeterministicChoiceIndices()), this->hasStateRewards() ? boost::optional<std::vector<T>>(this->getStateRewardVector()) : boost::optional<std::vector<T>>(), this->hasTransitionRewards() ? boost::optional<storm::storage::SparseMatrix<T>>(this->getTransitionRewardMatrix()) : boost::optional<storm::storage::SparseMatrix<T>>(), boost::optional<std::vector<std::set<uint_fast64_t>>>(this->getChoiceLabeling()));
+        return restrictedMdp;
+    }
+    
 	/*!
 	 * Calculates a hash over all values contained in this Model.
 	 * @return size_t A Hash Value
@@ -120,6 +174,7 @@ public:
 	virtual std::size_t getHash() const override {
 		return AbstractNondeterministicModel<T>::getHash();
 	}
+    
 private:
 
 	/*!
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 728ac1265..0623090ac 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -401,10 +401,6 @@ public:
 			triggerErrorState();
 			LOG4CPLUS_ERROR(logger, "Trying to initialize matrix that is not uninitialized.");
 			throw storm::exceptions::InvalidStateException("Trying to initialize matrix that is not uninitialized.");
-		} else if ((rowCount == 0) || (colCount == 0)) {
-			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Trying to create initialize a matrix with 0 rows or 0 columns.");
-			throw storm::exceptions::InvalidArgumentException("Trying to create initialize a matrix with 0 rows or 0 columns.");
 		} else if ((rowCount * colCount) < nonZeroEntries) {
 			triggerErrorState();
 			LOG4CPLUS_ERROR(logger, "Trying to initialize a matrix with more non-zero entries than there can be.");
@@ -426,7 +422,7 @@ public:
     
 	/*!
 	 * Sets the matrix element at the given row and column to the given value. After all elements have been added,
-     * a call to finalize() is mandatory.
+     * a call to finalize(false) is mandatory.
 	 * NOTE: This is a linear setter. It must be called consecutively for each element,
 	 * row by row *and* column by column.
      * NOTE: This method is different from insertNextValue(...) in that the number of nonzero elements must be known
@@ -464,7 +460,7 @@ public:
     
     /*!
      * Inserts a value at the given row and column with the given value. After all elements have been inserted,
-     * a call to finalize() is mandatory.
+     * a call to finalize(true) is mandatory.
      * NOTE: This is a linear inserter. It must be called consecutively for each element, row by row *and* column by
      * column.
      * NOTE: This method is different from addNextValue(...) in that the number of nonzero elements need not be known
@@ -473,22 +469,29 @@ public:
 	 * @param row The row in which the matrix element is to be set.
 	 * @param col The column in which the matrix element is to be set.
 	 * @param value The value that is to be set.
+     * @param pushRowIndication If set to true, the next row indication value is pushed, otherwise it is added. If the
+     * number of rows was not set in the beginning, then this needs to be true and false otherwise.
      */
-    void insertNextValue(const uint_fast64_t row, const uint_fast64_t col,	T const& value) {
+    void insertNextValue(const uint_fast64_t row, const uint_fast64_t col,	T const& value, bool pushRowIndication = false) {
 		// Check whether the given row and column positions are valid and throw
 		// error otherwise.
-		if ((row > rowCount) || (col > colCount)) {
+		if (row < lastRow) {
 			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Trying to insert a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ").");
-			throw storm::exceptions::OutOfRangeException() << "Trying to insert a value at illegal position (" << row << ", " << col << ") in matrix of size (" << rowCount << ", " << colCount << ").";
+			LOG4CPLUS_ERROR(logger, "Trying to insert a value at illegal position (" << row << ", " << col << ").");
+			throw storm::exceptions::OutOfRangeException() << "Trying to insert a value at illegal position (" << row << ", " << col << ").";
 		}
         
 		// If we switched to another row, we have to adjust the missing entries in the rowIndications array.
 		if (row != lastRow) {
 			for (uint_fast64_t i = lastRow + 1; i <= row; ++i) {
-				rowIndications[i] = currentSize;
+                if (pushRowIndication) {
+                    rowIndications.push_back(currentSize);
+                } else {
+                    rowIndications[i] = currentSize;
+                }
 			}
-			lastRow = row;
+            rowCount = row + 1;
+            lastRow = row;
 		}
         
 		// Finally, set the element and increase the current size.
@@ -496,14 +499,33 @@ public:
 		columnIndications.push_back(col);
         ++nonZeroEntryCount;
 		++currentSize;
+        
+        // Check that we also have the correct number of columns.
+        colCount = std::max(colCount, col + 1);
 	}
-
-
+    
+    /*!
+     * Inserts an empty row in the matrix.
+     *
+     * @param pushRowIndication If set to true, the next row indication value is pushed, otherwise it is added. If the
+     * number of rows was not set in the beginning, then this needs to be true and false otherwise.
+     */
+    void insertEmptyRow(bool pushRowIndication = false) {
+        if (pushRowIndication) {
+            rowIndications.push_back(currentSize);
+        } else {
+            rowIndications[lastRow + 1] = currentSize;
+        }
+        
+        ++rowCount;
+        ++lastRow;
+    }
+    
 	/*
 	 * Finalizes the sparse matrix to indicate that initialization has been completed and the matrix may now be used.
      *
      * @param pushSentinelElement A boolean flag that indicates whether the sentinel element is to be pushed or inserted
-     * at a fixed location. If the elements have been added to the matrix via insertNextElement, this needs to be true
+     * at a fixed location. If the elements have been added to the matrix via insertNextValue, this needs to be true
      * and false otherwise.
 	 */
 	void finalize(bool pushSentinelElement = false) {
@@ -519,10 +541,8 @@ public:
 		} else {
 			// Fill in the missing entries in the row_indications array.
 			// (Can happen because of empty rows at the end.)
-			if (lastRow != rowCount) {
-				for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) {
-					rowIndications[i] = currentSize;
-				}
+			for (uint_fast64_t i = lastRow + 1; i < rowCount; ++i) {
+				rowIndications[i] = currentSize;
 			}
 
 			// Set a sentinel element at the last position of the row_indications array. This eases iteration work, as
@@ -1308,18 +1328,18 @@ public:
 	}
 
 	/*!
-	 * Checks if the given matrix is a submatrix of the current matrix, where A matrix A is called a
+	 * Checks if the current matrix is a submatrix of the given matrix, where a matrix A is called a
 	 * submatrix of B if a value in A is only nonzero, if the value in B at the same position is
-	 * also nonzero. Furthermore, A and B have to have the same size.
+	 * also nonzero. Additionally, the matrices must be of equal size.
 	 *
-	 * @param matrix The matrix that is possibly a submatrix of the current matrix.
-	 * @returns True iff the given matrix is a submatrix of the current matrix.
+	 * @param matrix The matrix that possibly is a "supermatrix" of the current matrix.
+	 * @returns True iff the current matrix is a submatrix of the given matrix.
 	 */
-	bool containsAllPositionsOf(SparseMatrix<T> const& matrix) const {
-		// Check for mismatching sizes.
-		if (this->getRowCount() != matrix.getRowCount()) return false;
-		if (this->getColumnCount() != matrix.getColumnCount()) return false;
-
+	bool isSubmatrixOf(SparseMatrix<T> const& matrix) const {
+        // Check for matching sizes.
+        if (this->getRowCount() != matrix.getRowCount()) return false;
+        if (this->getColumnCount() != matrix.getColumnCount()) return false;
+        
 		// Check the subset property for all rows individually.
 		for (uint_fast64_t row = 0; row < this->getRowCount(); ++row) {
 			for (uint_fast64_t elem = rowIndications[row], elem2 = matrix.rowIndications[row]; elem < rowIndications[row + 1] && elem < matrix.rowIndications[row + 1]; ++elem) {
@@ -1328,7 +1348,9 @@ public:
 				while (elem2 < matrix.rowIndications[row + 1] && matrix.columnIndications[elem2] < columnIndications[elem]) {
 					++elem2;
 				}
-				if (!(elem2 < matrix.rowIndications[row + 1]) || columnIndications[elem] != matrix.columnIndications[elem2]) return false;
+				if (!(elem2 < matrix.rowIndications[row + 1]) || columnIndications[elem] != matrix.columnIndications[elem2]) {
+                    return false;
+                }
 			}
 		}
 		return true;
@@ -1501,8 +1523,9 @@ private:
 	}
 
 	/*!
-	 * Prepares the internal storage. For this, it requires the number of non-zero entries and the
-	 * amount of rows to be set correctly.
+	 * Prepares the internal storage. It relies on the number of non-zero entries and the
+	 * amount of rows to be set correctly. They may, however, be zero, but then insertNextValue needs to be used rather
+     * than addNextElement for filling the matrix.
 	 *
 	 * @param initializeElements If set to true, all entries are initialized.
 	 * @return True on success, false otherwise (allocation failed).
diff --git a/src/storm.cpp b/src/storm.cpp
index 19ec8baee..39497409a 100644
--- a/src/storm.cpp
+++ b/src/storm.cpp
@@ -26,7 +26,7 @@
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
-#include "src/counterexamples/MinimalLabelSetGenerator.h"
+#include "src/counterexamples/MILPMinimalLabelSetGenerator.h"
 #include "src/parser/AutoParser.h"
 #include "src/parser/PrctlParser.h"
 #include "src/utility/ErrorHandling.h"
diff --git a/src/utility/set.h b/src/utility/set.h
new file mode 100644
index 000000000..d264034eb
--- /dev/null
+++ b/src/utility/set.h
@@ -0,0 +1,52 @@
+/*
+ * set.h
+ *
+ *  Created on: 06.12.2012
+ *      Author: Christian Dehnert
+ */
+
+#ifndef STORM_UTILITY_SET_H_
+#define STORM_UTILITY_SET_H_
+
+#include <set>
+
+#include "log4cplus/logger.h"
+#include "log4cplus/loggingmacros.h"
+
+extern log4cplus::Logger logger;
+
+namespace storm {
+    namespace utility {
+        namespace set {
+            
+            template<typename T, typename Compare>
+            bool isSubsetOf(std::set<T, Compare> const& set1, std::set<T, Compare> const& set2) {
+                // First, get a comparator object.
+                typename std::set<T, Compare>::key_compare comparator = set1.key_comp();
+                
+                for (typename std::set<T, Compare>::const_iterator it1 = set1.begin(), it2 = set2.begin(); it1 != set1.end() && it2 != set2.end(); ++it1) {
+                    // If the value in set1 is smaller than the value in set2, set1 is not a subset of set2.
+                    if (comparator(*it1, *it2)) {
+                        return false;
+                    }
+                    
+                    // If the value in the second set is smaller, we need to move the iterator until the comparison is false.
+                    while(comparator(*it2, *it1) && it2 != set2.end()) {
+                        ++it2;
+                    }
+                        
+                    // If we have reached the end of set2 or the element we found is actually larger than the one in set1
+                    // we know that the subset property is violated.
+                    if (it2 == set2.end() || comparator(*it1, *it2)) {
+                        return false;
+                    }
+                        
+                    // Otherwise, we have found an equivalent element and can continue with the next one.
+                }
+            }
+            
+        } // namespace set
+    } // namespace utility
+} // namespace storm
+
+#endif /* STORM_UTILITY_SET_H_ */
diff --git a/test/functional/parser/ParsePrismTest.cpp b/test/functional/parser/ParsePrismTest.cpp
index 63eca2d31..0a8abb420 100644
--- a/test/functional/parser/ParsePrismTest.cpp
+++ b/test/functional/parser/ParsePrismTest.cpp
@@ -24,7 +24,7 @@ TEST(ParsePrismTest, parseTwoDice) {
 	storm::ir::Program program;
 	ASSERT_NO_THROW(program = storm::parser::PrismParserFromFile(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.nm"));
 
-	std::shared_ptr<storm::models::AbstractModel<double>> model = storm::adapters::ExplicitModelAdapter<double>::translateProgram(program);
+	std::shared_ptr<storm::models::AbstractModel<double>> model = storm::adapters::ExplicitModelAdapter<double>::translateProgram(program, "", "coinflips");
 	
 	ASSERT_EQ(model->getNumberOfStates(), 169ull);
 	ASSERT_EQ(model->as<storm::models::AbstractNondeterministicModel<double>>()->getNumberOfChoices(), 254ull);