diff --git a/examples/dft/seq1.dft b/examples/dft/seq1.dft
new file mode 100644
index 000000000..8f5459fd2
--- /dev/null
+++ b/examples/dft/seq1.dft
@@ -0,0 +1,5 @@
+toplevel "A";
+"A" and "B" "C";
+"X" seq "B" "C";
+"B" lambda=0.5 dorm=0.3;
+"C" lambda=0.5 dorm=0.3;
diff --git a/examples/dft/spare_cold.dft b/examples/dft/spare_cold.dft
new file mode 100644
index 000000000..d01174bd5
--- /dev/null
+++ b/examples/dft/spare_cold.dft
@@ -0,0 +1,5 @@
+toplevel "A";
+"A" wsp "I" "M";
+"I" lambda=0.5 dorm=0.0;
+"M" lambda=0.5 dorm=0.0;
+
diff --git a/examples/dft/symmetry3.dft b/examples/dft/symmetry3.dft
new file mode 100644
index 000000000..edda963f9
--- /dev/null
+++ b/examples/dft/symmetry3.dft
@@ -0,0 +1,12 @@
+toplevel "A";
+"A" and "B" "B'" "B''";
+"B" and "C" "D";
+"B'" and "C'" "D'";
+"B''" and "C''" "D''";
+"C" lambda=0.5 dorm=0;
+"D" lambda=0.5 dorm=0;
+"C'" lambda=0.5 dorm=0;
+"D'" lambda=0.5 dorm=0;
+"C''" lambda=0.5 dorm=0;
+"D''" lambda=0.5 dorm=0;
+
diff --git a/examples/dft/symmetry4.dft b/examples/dft/symmetry4.dft
new file mode 100644
index 000000000..f401fc241
--- /dev/null
+++ b/examples/dft/symmetry4.dft
@@ -0,0 +1,12 @@
+toplevel "A";
+"A" and "B" "B'" "C" "C'";
+"B" and "D" "E";
+"B'" and "D'" "E'";
+"C" or "F";
+"C'" or "F'";
+"D" lambda=0.5 dorm=0;
+"E" lambda=0.5 dorm=0;
+"D'" lambda=0.5 dorm=0;
+"E'" lambda=0.5 dorm=0;
+"F" lambda=0.5 dorm=0;
+"F'" lambda=0.5 dorm=0;
diff --git a/examples/dft/symmetry5.dft b/examples/dft/symmetry5.dft
new file mode 100644
index 000000000..18f700bb2
--- /dev/null
+++ b/examples/dft/symmetry5.dft
@@ -0,0 +1,21 @@
+toplevel "A";
+"A" and "BA" "BB" "BC" "BD" "BE" "BF";
+"BA" and "CA" "DA";
+"BB" and "CB" "DB";
+"BC" and "CC" "DC";
+"BD" and "CD" "DD";
+"BE" and "CE" "DE";
+"BF" and "CF" "DF";
+"CA" lambda=0.5 dorm=0;
+"DA" lambda=0.5 dorm=0;
+"CB" lambda=0.5 dorm=0;
+"DB" lambda=0.5 dorm=0;
+"CC" lambda=0.5 dorm=0;
+"DC" lambda=0.5 dorm=0;
+"CD" lambda=0.5 dorm=0;
+"DD" lambda=0.5 dorm=0;
+"CE" lambda=0.5 dorm=0;
+"DE" lambda=0.5 dorm=0;
+"CF" lambda=0.5 dorm=0;
+"DF" lambda=0.5 dorm=0;
+
diff --git a/src/builder/ExplicitDFTModelBuilder.cpp b/src/builder/ExplicitDFTModelBuilder.cpp
index 4900845a5..f14408190 100644
--- a/src/builder/ExplicitDFTModelBuilder.cpp
+++ b/src/builder/ExplicitDFTModelBuilder.cpp
@@ -15,78 +15,57 @@ namespace storm {
         }
         
         template <typename ValueType>
-        ExplicitDFTModelBuilder<ValueType>::ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft) : mDft(dft), mStates(((mDft.stateVectorSize() / 64) + 1) * 64, INITIAL_BUCKETSIZE) {
+        ExplicitDFTModelBuilder<ValueType>::ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTIndependentSymmetries const& symmetries) : mDft(dft), mStates(((mDft.stateVectorSize() / 64) + 1) * 64, INITIAL_BUCKETSIZE) {
             // stateVectorSize is bound for size of bitvector
 
-            // Find symmetries
-            // TODO activate
-            // Currently using hack to test
-            std::vector<std::vector<size_t>> symmetriesTmp;
-            std::vector<size_t> vecB;
-            std::vector<size_t> vecC;
-            std::vector<size_t> vecD;
-            if (false) {
-                for (size_t i = 0; i < mDft.nrElements(); ++i) {
-                    std::string name = mDft.getElement(i)->name();
-                    size_t id = mDft.getElement(i)->id();
-                    if (boost::starts_with(name, "B")) {
-                        vecB.push_back(id);
-                    } else if (boost::starts_with(name, "C")) {
-                        vecC.push_back(id);
-                    } else if (boost::starts_with(name, "D")) {
-                        vecD.push_back(id);
-                    }
-                }
-                symmetriesTmp.push_back(vecB);
-                symmetriesTmp.push_back(vecC);
-                symmetriesTmp.push_back(vecD);
-                std::cout << "Found the following symmetries:" << std::endl;
-                for (auto const& symmetry : symmetriesTmp) {
-                    for (auto const& elem : symmetry) {
-                        std::cout << elem << " -> ";
-                    }
-                    std::cout << std::endl;
-                }
-            } else {
-                vecB.push_back(mDft.getTopLevelIndex());
-            }
-            mStateGenerationInfo = std::make_shared<storm::storage::DFTStateGenerationInfo>(mDft.buildStateGenerationInfo(vecB, symmetriesTmp));
+            mStateGenerationInfo = std::make_shared<storm::storage::DFTStateGenerationInfo>(mDft.buildStateGenerationInfo(symmetries));
         }
 
 
         template <typename ValueType>
         std::shared_ptr<storm::models::sparse::Model<ValueType>> ExplicitDFTModelBuilder<ValueType>::buildModel(LabelOptions const& labelOpts) {
             // Initialize
-            DFTStatePointer state = std::make_shared<storm::storage::DFTState<ValueType>>(mDft, *mStateGenerationInfo, newIndex++);
-            mStates.findOrAdd(state->status(), state->getId());
-
-            std::queue<DFTStatePointer> stateQueue;
-            stateQueue.push(state);
-
             bool deterministicModel = false;
             ModelComponents modelComponents;
             std::vector<uint_fast64_t> tmpMarkovianStates;
             storm::storage::SparseMatrixBuilder<ValueType> transitionMatrixBuilder(0, 0, 0, false, !deterministicModel, 0);
+            
 
+            if(mergeFailedStates) {
+                newIndex++;
+                
+                transitionMatrixBuilder.newRowGroup(failedIndex);
+                transitionMatrixBuilder.addNextValue(failedIndex, failedIndex, storm::utility::one<ValueType>());
+                STORM_LOG_TRACE("Added self loop for " << failedIndex);
+                modelComponents.exitRates.push_back(storm::utility::one<ValueType>());
+                tmpMarkovianStates.push_back(failedIndex);        
+            }
+            
+            DFTStatePointer state = std::make_shared<storm::storage::DFTState<ValueType>>(mDft, *mStateGenerationInfo, newIndex++);
+            mStates.findOrAdd(state->status(), state->getId());
+            initialStateIndex = state->getId();
+            std::queue<DFTStatePointer> stateQueue;
+            stateQueue.push(state);
+            
             // Begin model generation
             bool deterministic = exploreStates(stateQueue, transitionMatrixBuilder, tmpMarkovianStates, modelComponents.exitRates);
-            STORM_LOG_DEBUG("Generated " << mStates.size() << " states");
+            STORM_LOG_DEBUG("Generated " << mStates.size() + (mergeFailedStates ? 1 : 0) << " states");
             STORM_LOG_DEBUG("Model is " << (deterministic ? "deterministic" : "non-deterministic"));
 
+            size_t stateSize = mStates.size() + (mergeFailedStates ? 1 : 0);
             // Build Markov Automaton
-            modelComponents.markovianStates = storm::storage::BitVector(mStates.size(), tmpMarkovianStates);
+            modelComponents.markovianStates = storm::storage::BitVector(stateSize, tmpMarkovianStates);
             // Build transition matrix
-            modelComponents.transitionMatrix = transitionMatrixBuilder.build();
+            modelComponents.transitionMatrix = transitionMatrixBuilder.build(stateSize, stateSize);
             STORM_LOG_DEBUG("Transition matrix: " << std::endl << modelComponents.transitionMatrix);
             STORM_LOG_DEBUG("Exit rates: " << modelComponents.exitRates);
             STORM_LOG_DEBUG("Markovian states: " << modelComponents.markovianStates);
-            assert(!deterministic || modelComponents.transitionMatrix.getRowCount() == modelComponents.transitionMatrix.getColumnCount());
-
+            
             // Build state labeling
-            modelComponents.stateLabeling = storm::models::sparse::StateLabeling(mStates.size());
+            modelComponents.stateLabeling = storm::models::sparse::StateLabeling(mStates.size() + (mergeFailedStates ? 1 : 0));
             // Initial state is always first state without any failure
             modelComponents.stateLabeling.addLabel("init");
-            modelComponents.stateLabeling.addLabelToState("init", 0);
+            modelComponents.stateLabeling.addLabelToState("init", initialStateIndex);
             // Label all states corresponding to their status (failed, failsafe, failed BE)
             if(labelOpts.buildFailLabel) {
                 modelComponents.stateLabeling.addLabel("failed");
@@ -103,10 +82,13 @@ namespace storm {
                 }
             }
 
+            if(mergeFailedStates) {
+                modelComponents.stateLabeling.addLabelToState("failed", failedIndex);
+            }
             for (auto const& stateIdPair : mStates) {
                 storm::storage::BitVector state = stateIdPair.first;
                 size_t stateId = stateIdPair.second;
-                if (labelOpts.buildFailLabel && mDft.hasFailed(state, *mStateGenerationInfo)) {
+                if (!mergeFailedStates && labelOpts.buildFailLabel && mDft.hasFailed(state, *mStateGenerationInfo)) {
                     modelComponents.stateLabeling.addLabelToState("failed", stateId);
                 }
                 if (labelOpts.buildFailSafeLabel && mDft.isFailsafe(state, *mStateGenerationInfo)) {
@@ -151,13 +133,12 @@ namespace storm {
         template <typename ValueType>
         bool ExplicitDFTModelBuilder<ValueType>::exploreStates(std::queue<DFTStatePointer>& stateQueue, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<uint_fast64_t>& markovianStates, std::vector<ValueType>& exitRates) {
 
-            assert(exitRates.empty());
-            assert(markovianStates.empty());
             // TODO Matthias: set Markovian states directly as bitvector?
-            std::map<size_t, ValueType> outgoingTransitions;
+            std::map<uint_fast64_t, ValueType> outgoingTransitions;
             size_t rowOffset = 0; // Captures number of non-deterministic choices
             bool deterministic = true;
-
+            
+            
             while (!stateQueue.empty()) {
                 // Initialization
                 outgoingTransitions.clear();
@@ -189,6 +170,7 @@ namespace storm {
 
                 // Let BE fail
                 while (smallest < failableCount) {
+                    assert(!mDft.hasFailed(state));
                     STORM_LOG_TRACE("exploring from: " << mDft.getStateString(state));
 
                     // Construct new state as copy from original one
@@ -207,11 +189,23 @@ namespace storm {
                             queues.propagateFailure(parent);
                         }
                     }
+                    for (DFTRestrictionPointer restr : nextBE->restrictions()) {
+                        queues.checkRestrictionLater(restr);
+                    }
 
                     while (!queues.failurePropagationDone()) {
                         DFTGatePointer next = queues.nextFailurePropagation();
                         next->checkFails(*newState, queues);
                     }
+                    
+                    while(!queues.restrictionChecksDone()) {
+                        DFTRestrictionPointer next = queues.nextRestrictionCheck();
+                        next->checkFails(*newState, queues);
+                    } 
+                    
+                    if(newState->isInvalid()) {
+                        continue;
+                    }
 
                     while (!queues.failsafePropagationDone()) {
                         DFTGatePointer next = queues.nextFailsafePropagation();
@@ -226,19 +220,52 @@ namespace storm {
                     // Update failable dependencies
                     newState->updateFailableDependencies(nextBE->id());
                     
-                    size_t newStateId;
-                    if (mStates.contains(newState->status())) {
-                        // State already exists
-                        newStateId = mStates.getValue(newState->status());
-                        STORM_LOG_TRACE("State " << mDft.getStateString(newState) << " already exists");
+                    bool dftFailed = newState->hasFailed(mDft.getTopLevelIndex());
+                    uint_fast64_t newStateId;
+                    if(dftFailed && mergeFailedStates) {
+                        newStateId = failedIndex;
                     } else {
-                        // New state
-                        newState->setId(newIndex++);
-                        newStateId = mStates.findOrAdd(newState->status(), newState->getId());
-                        STORM_LOG_TRACE("New state " << mDft.getStateString(newState));
-
-                        // Add state to search queue
-                        stateQueue.push(newState);
+                        // Order state by symmetry
+                        STORM_LOG_TRACE("Check for symmetry: " << mDft.getStateString(newState));
+                        bool changed = newState->orderBySymmetry();
+                        STORM_LOG_TRACE("State " << (changed ? "changed to " : "did not change") << (changed ? mDft.getStateString(newState) : ""));
+                        
+                        // Check if state already exists
+                        if (mStates.contains(newState->status())) {
+                            // State already exists
+                            newStateId = mStates.getValue(newState->status());
+                            STORM_LOG_TRACE("State " << mDft.getStateString(newState) << " with id " << newStateId << " already exists");
+                            
+                            // Check if possible pseudo state can be created now
+                            if (!changed && newStateId >= OFFSET_PSEUDO_STATE) {
+                                newStateId = newStateId - OFFSET_PSEUDO_STATE;
+                                assert(newStateId < mPseudoStatesMapping.size());
+                                if (mPseudoStatesMapping[newStateId] == 0) {
+                                    // Create pseudo state now
+                                    newState->setId(newIndex++);
+                                    mPseudoStatesMapping[newStateId] = newState->getId();
+                                    newStateId = newState->getId();
+                                    mStates.setOrAdd(newState->status(), newStateId);
+                                    stateQueue.push(newState);
+                                    STORM_LOG_TRACE("Now create state " << mDft.getStateString(newState) << " with id " << newStateId);
+                                }
+                            }
+                        } else {
+                            // New state
+                            if (changed) {
+                                // Remember to create state later on
+                                newState->setId(mPseudoStatesMapping.size() + OFFSET_PSEUDO_STATE);
+                                mPseudoStatesMapping.push_back(0);
+                                newStateId = mStates.findOrAdd(newState->status(), newState->getId());
+                                STORM_LOG_TRACE("New state for later creation: " << mDft.getStateString(newState));
+                            } else {
+                                // Create new state
+                                newState->setId(newIndex++);
+                                newStateId = mStates.findOrAdd(newState->status(), newState->getId());
+                                STORM_LOG_TRACE("New state: " << mDft.getStateString(newState));
+                                stateQueue.push(newState);
+                            }
+                        }
                     }
 
                     // Set transitions
@@ -249,10 +276,11 @@ namespace storm {
                         STORM_LOG_TRACE("Added transition from " << state->getId() << " to " << newStateId << " with probability " << dependency->probability());
 
                         if (!storm::utility::isOne(dependency->probability())) {
+                            // TODO Matthias: use symmetry as well
                             // Add transition to state where dependency was unsuccessful
                             DFTStatePointer unsuccessfulState = std::make_shared<storm::storage::DFTState<ValueType>>(*state);
                             unsuccessfulState->letDependencyBeUnsuccessful(smallest-1);
-                            size_t unsuccessfulStateId;
+                            uint_fast64_t unsuccessfulStateId;
                             if (mStates.contains(unsuccessfulState->status())) {
                                 // Unsuccessful state already exists
                                 unsuccessfulStateId = mStates.getValue(unsuccessfulState->status());
@@ -304,6 +332,38 @@ namespace storm {
                 if (hasDependencies) {
                     rowOffset--; // One increment to many
                 } else {
+                    // Try to merge pseudo states with their instantiation
+                    // TODO Matthias: improve?
+                    auto it = outgoingTransitions.begin();
+                    while (it != outgoingTransitions.end()) {
+                        if (it->first >= OFFSET_PSEUDO_STATE) {
+                            uint_fast64_t newId = it->first - OFFSET_PSEUDO_STATE;
+                            assert(newId < mPseudoStatesMapping.size());
+                            if (mPseudoStatesMapping[newId] > 0) {
+                                // State exists already
+                                newId = mPseudoStatesMapping[newId];
+                                auto itFind = outgoingTransitions.find(newId);
+                                if (itFind != outgoingTransitions.end()) {
+                                    // Add probability from pseudo state to instantiation
+                                    itFind->second += it->second;
+                                    STORM_LOG_TRACE("Merged pseudo state " << newId << " adding rate " << it->second << " to total rate of " << itFind->second);
+                                } else {
+                                    // Only change id
+                                    outgoingTransitions.emplace(newId, it->second);
+                                    STORM_LOG_TRACE("Instantiated pseudo state " << newId << " with rate " << it->second);
+                                }
+                                // Remove pseudo state
+                                auto itErase = it;
+                                ++it;
+                                outgoingTransitions.erase(itErase);
+                            } else {
+                                ++it;
+                            }
+                        } else {
+                            ++it;
+                        }
+                    }
+                    
                     // Add all probability transitions
                     for (auto it = outgoingTransitions.begin(); it != outgoingTransitions.end(); ++it)
                     {
@@ -317,6 +377,11 @@ namespace storm {
                 assert(exitRates.size()-1 == state->getId());
 
             } // end while queue
+            assert(std::find(mPseudoStatesMapping.begin(), mPseudoStatesMapping.end(), 0) == mPseudoStatesMapping.end());
+            
+            // Replace pseudo states in matrix
+            transitionMatrixBuilder.replaceColumns(mPseudoStatesMapping, OFFSET_PSEUDO_STATE);
+            
             assert(!deterministic || rowOffset == 0);
             return deterministic;
         }
diff --git a/src/builder/ExplicitDFTModelBuilder.h b/src/builder/ExplicitDFTModelBuilder.h
index a64efaddd..13b3eaf0d 100644
--- a/src/builder/ExplicitDFTModelBuilder.h
+++ b/src/builder/ExplicitDFTModelBuilder.h
@@ -1,13 +1,13 @@
 #ifndef EXPLICITDFTMODELBUILDER_H
 #define	EXPLICITDFTMODELBUILDER_H
 
-#include "../storage/dft/DFT.h"
-
 #include <src/models/sparse/StateLabeling.h>
 #include <src/models/sparse/StandardRewardModel.h>
 #include <src/models/sparse/Model.h>
 #include <src/storage/SparseMatrix.h>
 #include <src/storage/BitVectorHashMap.h>
+#include <src/storage/dft/DFT.h>
+#include <src/storage/dft/SymmetricUnits.h>
 #include <boost/container/flat_set.hpp>
 #include <boost/optional/optional.hpp>
 #include <stack>
@@ -23,6 +23,8 @@ namespace storm {
             using DFTElementCPointer = std::shared_ptr<storm::storage::DFTElement<ValueType> const>;
             using DFTGatePointer = std::shared_ptr<storm::storage::DFTGate<ValueType>>;
             using DFTStatePointer = std::shared_ptr<storm::storage::DFTState<ValueType>>;
+            using DFTRestrictionPointer = std::shared_ptr<storm::storage::DFTRestriction<ValueType>>;
+
 
             // A structure holding the individual components of a model.
             struct ModelComponents {
@@ -45,13 +47,16 @@ namespace storm {
             };
             
             const size_t INITIAL_BUCKETSIZE = 20000;
-
-          
+            const uint_fast64_t OFFSET_PSEUDO_STATE = UINT_FAST64_MAX / 2;
             
             storm::storage::DFT<ValueType> const& mDft;
             std::shared_ptr<storm::storage::DFTStateGenerationInfo> mStateGenerationInfo;
-            storm::storage::BitVectorHashMap<uint32_t> mStates;
+            storm::storage::BitVectorHashMap<uint_fast64_t> mStates;
+            std::vector<uint_fast64_t> mPseudoStatesMapping;
             size_t newIndex = 0;
+            bool mergeFailedStates = true;
+            size_t failedIndex = 0;
+            size_t initialStateIndex = 0;
 
         public:
               struct LabelOptions {
@@ -60,7 +65,7 @@ namespace storm {
                 std::set<std::string> beLabels = {};
             };
             
-            ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft);
+            ExplicitDFTModelBuilder(storm::storage::DFT<ValueType> const& dft, storm::storage::DFTIndependentSymmetries const& symmetries);
 
             std::shared_ptr<storm::models::sparse::Model<ValueType>> buildModel(LabelOptions const& labelOpts);
 
diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
index 2ea8fa459..5425044b2 100644
--- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
+++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
@@ -584,11 +584,13 @@ namespace storm {
             
             // Do some sanity checks to establish some required properties.
             RewardModelType const& rewardModel = this->getModel().getRewardModel(checkTask.isRewardModelSet() ? checkTask.getRewardModel() : "");
-            // TODO use maybeStates instead of all states
-            std::vector<ValueType> stateRewardValues = rewardModel.getTotalRewardVector(trueStates.getNumberOfSetBits(), this->getModel().getTransitionMatrix(), trueStates);
 
             STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::IllegalArgumentException, "Input model does not have a reward model.");
-            std::vector<ValueType> result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates, stateRewardValues, checkTask.isOnlyInitialStatesRelevantSet());
+            std::vector<ValueType> result = computeReachabilityRewards(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), targetStates,
+                                                                       [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) {
+                                                                           return rewardModel.getTotalRewardVector(numberOfRows, transitionMatrix, maybeStates);
+                                                                       },
+                                                                       checkTask.isOnlyInitialStatesRelevantSet());
 
             // Construct check result.
             std::unique_ptr<CheckResult> checkResult(new ExplicitQuantitativeCheckResult<ValueType>(result));
@@ -597,6 +599,17 @@ namespace storm {
 
         template<typename SparseDtmcModelType>
         std::vector<typename SparseDtmcModelType::ValueType> SparseDtmcEliminationModelChecker<SparseDtmcModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::vector<ValueType>& stateRewardValues, bool computeForInitialStatesOnly) {
+            return computeReachabilityRewards(probabilityMatrix, backwardTransitions, initialStates, targetStates,
+                                              [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) {
+                                                  std::vector<ValueType> result(numberOfRows);
+                                                  storm::utility::vector::selectVectorValues(result, maybeStates, stateRewardValues);
+                                                  return result;
+                                              },
+                                              computeForInitialStatesOnly);
+        }
+
+        template<typename SparseDtmcModelType>
+        std::vector<typename SparseDtmcModelType::ValueType> SparseDtmcEliminationModelChecker<SparseDtmcModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly) {
 
             uint_fast64_t numberOfStates = probabilityMatrix.getRowCount();
             
@@ -639,9 +652,9 @@ namespace storm {
                 storm::storage::SparseMatrix<ValueType> submatrixTransposed = submatrix.transpose();
                 
                 // Project the state reward vector to all maybe-states.
-                storm::storage::BitVector phiStates(numberOfStates, true);
+                std::vector<ValueType> stateRewardValues = totalStateRewardVectorGetter(submatrix.getRowCount(), probabilityMatrix, maybeStates);
 
-                std::vector<ValueType> subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, phiStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates));
+                std::vector<ValueType> subresult = computeReachabilityValues(submatrix, stateRewardValues, submatrixTransposed, newInitialStates, computeForInitialStatesOnly, trueStates, targetStates, probabilityMatrix.getConstrainedRowSumVector(maybeStates, targetStates));
                 storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, subresult);
             }
             
diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h
index 76e939d66..33b61a2c4 100644
--- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h
+++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.h
@@ -91,6 +91,8 @@ namespace storm {
             
             static std::vector<ValueType> computeLongRunValues(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& maybeStates, bool computeResultsForInitialStatesOnly, std::vector<ValueType>& stateValues);
             
+            static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, bool computeForInitialStatesOnly);
+
             static std::vector<ValueType> computeReachabilityValues(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType>& values, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& initialStates, bool computeResultsForInitialStatesOnly, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& oneStepProbabilitiesToTarget);
             
             static std::shared_ptr<StatePriorityQueue<ValueType>> createStatePriorityQueue(boost::optional<std::vector<uint_fast64_t>> const& stateDistances, storm::storage::FlexibleSparseMatrix<ValueType> const& transitionMatrix, storm::storage::FlexibleSparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType>& oneStepProbabilities, storm::storage::BitVector const& states);
diff --git a/src/parser/DFTGalileoParser.cpp b/src/parser/DFTGalileoParser.cpp
index d11f5ae3b..65ad82437 100644
--- a/src/parser/DFTGalileoParser.cpp
+++ b/src/parser/DFTGalileoParser.cpp
@@ -102,6 +102,8 @@ namespace storm {
                         success = builder.addPandElement(name, childNames);
                     } else if (tokens[1] == "wsp" || tokens[1] == "csp") {
                         success = builder.addSpareElement(name, childNames);
+                    } else if (tokens[1] == "seq") {
+                        success = builder.addSequenceEnforcer(name, childNames);
                     } else if (tokens[1] == "fdep") {
                         success = builder.addDepElement(name, childNames, storm::utility::one<ValueType>());
                     } else if (boost::starts_with(tokens[1], "pdep=")) {
diff --git a/src/storage/BitVectorHashMap.cpp b/src/storage/BitVectorHashMap.cpp
index 9237dfde3..07eefee1d 100644
--- a/src/storage/BitVectorHashMap.cpp
+++ b/src/storage/BitVectorHashMap.cpp
@@ -140,6 +140,11 @@ namespace storm {
             return findOrAddAndGetBucket(key, value).first;
         }
         
+        template<class ValueType, class Hash1, class Hash2>
+        void BitVectorHashMap<ValueType, Hash1, Hash2>::setOrAdd(storm::storage::BitVector const& key, ValueType const& value) {
+            setOrAddAndGetBucket(key, value);
+        }
+        
         template<class ValueType, class Hash1, class Hash2>
         std::pair<ValueType, std::size_t> BitVectorHashMap<ValueType, Hash1, Hash2>::findOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value) {
             // If the load of the map is too high, we increase the size.
@@ -161,6 +166,25 @@ namespace storm {
             }
         }
         
+        template<class ValueType, class Hash1, class Hash2>
+        std::size_t BitVectorHashMap<ValueType, Hash1, Hash2>::setOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value) {
+            // If the load of the map is too high, we increase the size.
+            if (numberOfElements >= loadFactor * *currentSizeIterator) {
+                this->increaseSize();
+            }
+            
+            std::tuple<bool, std::size_t, bool> flagBucketTuple = this->findBucketToInsert<true>(key);
+            STORM_LOG_ASSERT(!std::get<2>(flagBucketTuple), "Failed to find bucket for insertion.");
+            if (!std::get<0>(flagBucketTuple)) {
+                // Insert the new bits into the bucket.
+                buckets.set(std::get<1>(flagBucketTuple) * bucketSize, key);
+                occupied.set(std::get<1>(flagBucketTuple));
+                ++numberOfElements;
+            }
+            values[std::get<1>(flagBucketTuple)] = value;
+            return std::get<1>(flagBucketTuple);
+        }
+        
         template<class ValueType, class Hash1, class Hash2>
         ValueType BitVectorHashMap<ValueType, Hash1, Hash2>::getValue(storm::storage::BitVector const& key) const {
             std::pair<bool, std::size_t> flagBucketPair = this->findBucket(key);
diff --git a/src/storage/BitVectorHashMap.h b/src/storage/BitVectorHashMap.h
index d50dac59e..18ca52186 100644
--- a/src/storage/BitVectorHashMap.h
+++ b/src/storage/BitVectorHashMap.h
@@ -66,6 +66,15 @@ namespace storm {
              * @return The found value if the key is already contained in the map and the provided new value otherwise.
              */
             ValueType findOrAdd(storm::storage::BitVector const& key, ValueType const& value);
+            
+            /*!
+             * Sets the given key value pain in the map. If the key is found in the map, the corresponding value is
+             * overwritten with the given value. Otherwise, the key is inserted with the given value.
+             *
+             * @param key The key to search or insert.
+             * @param value The value to set.
+             */
+            void setOrAdd(storm::storage::BitVector const& key, ValueType const& value);
 
             /*!
              * Searches for the given key in the map. If it is found, the mapped-to value is returned. Otherwise, the
@@ -79,6 +88,16 @@ namespace storm {
              */
             std::pair<ValueType, std::size_t> findOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value);
             
+            /*!
+             * Sets the given key value pain in the map. If the key is found in the map, the corresponding value is
+             * overwritten with the given value. Otherwise, the key is inserted with the given value.
+             *
+             * @param key The key to search or insert.
+             * @param value The value to set.
+             * @return The index of the bucket into which the key was inserted.
+             */
+            std::size_t setOrAddAndGetBucket(storm::storage::BitVector const& key, ValueType const& value);
+            
             /*!
              * Retrieves the key stored in the given bucket (if any) and the value it is mapped to.
              *
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index d41e2ef2c..fa95f3930 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -220,6 +220,73 @@ namespace storm {
             return lastColumn;
         }
         
+        // Debug method for printing the current matrix
+        template<typename ValueType>
+        void print(std::vector<typename SparseMatrix<ValueType>::index_type> const& rowGroupIndices, std::vector<MatrixEntry<typename SparseMatrix<ValueType>::index_type, typename SparseMatrix<ValueType>::value_type>> const& columnsAndValues, std::vector<typename SparseMatrix<ValueType>::index_type> const& rowIndications) {
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::index_type endRows;
+            // Iterate over all row groups.
+            for (typename SparseMatrix<ValueType>::index_type group = 0; group < rowGroupIndices.size(); ++group) {
+                std::cout << "\t---- group " << group << "/" << (rowGroupIndices.size() - 1) << " ---- " << std::endl;
+                endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size();
+                // Iterate over all rows in a row group
+                for (typename SparseMatrix<ValueType>::index_type i = rowGroupIndices[group]; i < endGroups; ++i) {
+                    endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size();
+                    // Print the actual row.
+                    std::cout << "Row " << i << " (" << rowIndications[i] << " - " << endRows << ")" << ": ";
+                    for (typename SparseMatrix<ValueType>::index_type pos = rowIndications[i]; pos < endRows; ++pos) {
+                        std::cout << "(" << columnsAndValues[pos].getColumn() << ": " << columnsAndValues[pos].getValue() << ") ";
+                    }
+                    std::cout << std::endl;
+                }
+            }
+        }
+        
+        template<typename ValueType>
+        bool SparseMatrixBuilder<ValueType>::replaceColumns(std::vector<index_type> const& replacements, index_type offset) {
+            bool changed = false;
+            index_type maxColumn = 0;
+            for (auto& elem : columnsAndValues) {
+                if (elem.getColumn() >= offset) {
+                    elem.setColumn(replacements[elem.getColumn() - offset]);
+                    changed = true;
+                }
+                maxColumn = std::max(maxColumn, elem.getColumn());
+            }
+            assert(changed || highestColumn == maxColumn);
+            highestColumn = maxColumn;
+            assert(changed || lastColumn == columnsAndValues[columnsAndValues.size() - 1].getColumn());
+            lastColumn = columnsAndValues[columnsAndValues.size() - 1].getColumn();
+            
+            if (changed) {
+                fixColumns();
+            }
+            return changed;
+        }
+
+        template<typename ValueType>
+        void SparseMatrixBuilder<ValueType>::fixColumns() {
+            // Sort columns per row
+            typename SparseMatrix<ValueType>::index_type endGroups;
+            typename SparseMatrix<ValueType>::index_type endRows;
+            for (index_type group = 0; group < rowGroupIndices.size(); ++group) {
+                endGroups = group < rowGroupIndices.size()-1 ? rowGroupIndices[group+1] : rowIndications.size();
+                for (index_type i = rowGroupIndices[group]; i < endGroups; ++i) {
+                    endRows = i < rowIndications.size()-1 ? rowIndications[i+1] : columnsAndValues.size();
+                    // Sort the row
+                    std::sort(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                              [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                  return a.getColumn() < b.getColumn();
+                              });
+                    // Assert no equal elements
+                    assert(std::is_sorted(columnsAndValues.begin() + rowIndications[i], columnsAndValues.begin() + endRows,
+                                          [](MatrixEntry<index_type, value_type> const& a, MatrixEntry<index_type, value_type> const& b) {
+                                              return a.getColumn() <= b.getColumn();
+                                          }));
+                }
+            }
+        }
+        
         template<typename ValueType>
         SparseMatrix<ValueType>::rows::rows(iterator begin, index_type entryCount) : beginIterator(begin), entryCount(entryCount) {
             // Intentionally left empty.
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index f9d9916c4..4892dbe32 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -217,6 +217,17 @@ namespace storm {
              */
             index_type getLastColumn() const;
             
+            /*!
+             * Replaces all columns with id > offset according to replacements.
+             * Every state  with id offset+i is replaced by the id in replacements[i].
+             * Afterwards the columns are sorted.
+             *
+             * @param replacements Mapping indicating the replacements from offset+i -> value of i.
+             * @param offset Offset to add to each id in vector index.
+             * @return True if replacement took place, False if nothing changed.
+             */
+            bool replaceColumns(std::vector<index_type> const& replacements, index_type offset);
+            
         private:
             // A flag indicating whether a row count was set upon construction.
             bool initialRowCountSet;
@@ -277,6 +288,11 @@ namespace storm {
             // Stores the currently active row group. This is used for correctly constructing the row grouping of the
             // matrix.
             index_type currentRowGroup;
+            
+            /*!
+             * Fixes the matrix by sorting the columns to gain increasing order again.
+             */
+            void fixColumns();
         };
         
         /*!
diff --git a/src/storage/dft/DFT.cpp b/src/storage/dft/DFT.cpp
index 3fbcd5ee1..eef83928d 100644
--- a/src/storage/dft/DFT.cpp
+++ b/src/storage/dft/DFT.cpp
@@ -68,20 +68,77 @@ namespace storm {
         }
 
         template<typename ValueType>
-        DFTStateGenerationInfo DFT<ValueType>::buildStateGenerationInfo(std::vector<size_t> const& subTreeRoots, std::vector<std::vector<size_t>> const& symmetries) const {
+        DFTStateGenerationInfo DFT<ValueType>::buildStateGenerationInfo(storm::storage::DFTIndependentSymmetries const& symmetries) const {
             // Use symmetry
             // Collect all elements in the first subtree
             // TODO make recursive to use for nested subtrees
-            
+
             DFTStateGenerationInfo generationInfo(nrElements());
 
             // Perform DFS and insert all elements of subtree sequentially
             size_t stateIndex = 0;
             std::queue<size_t> visitQueue;
-            std::set<size_t> visited;
-            visitQueue.push(subTreeRoots[0]);
-            stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
-            
+            storm::storage::BitVector visited(nrElements(), false);
+
+            if (symmetries.groups.empty()) {
+                // Perform DFS for whole tree
+                visitQueue.push(mTopLevelIndex);
+                stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
+            } else {
+                for (auto const& symmetryGroup : symmetries.groups) {
+                    assert(!symmetryGroup.second.empty());
+
+                    // Perform DFS for first subtree of each symmetry
+                    visitQueue.push(symmetryGroup.first);
+                    size_t groupIndex = stateIndex;
+                    stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
+                    size_t offset = stateIndex - groupIndex;
+                    
+                    // Mirror symmetries
+                    size_t noSymmetricElements = symmetryGroup.second.front().size();
+                    assert(noSymmetricElements > 1);
+
+                    for (std::vector<size_t> symmetricElements : symmetryGroup.second) {
+                        assert(symmetricElements.size() == noSymmetricElements);
+
+                        // Initialize for original element
+                        size_t originalElement = symmetricElements[0];
+                        size_t index = generationInfo.getStateIndex(originalElement);
+                        size_t activationIndex = isRepresentative(originalElement) ? generationInfo.getSpareActivationIndex(originalElement) : 0;
+                        size_t usageIndex = mElements[originalElement]->isSpareGate() ? generationInfo.getSpareUsageIndex(originalElement) : 0;
+                        
+                        // Mirror symmetry for each element
+                        for (size_t i = 1; i < symmetricElements.size(); ++i) {
+                            size_t symmetricElement = symmetricElements[i];
+                            visited.set(symmetricElement);
+
+                            generationInfo.addStateIndex(symmetricElement, index + offset * i);
+                            stateIndex += 2;
+
+                            assert((activationIndex > 0) == isRepresentative(symmetricElement));
+                            if (activationIndex > 0) {
+                                generationInfo.addSpareActivationIndex(symmetricElement, activationIndex + offset * i);
+                                ++stateIndex;
+                            }
+
+                            assert((usageIndex > 0) == mElements[symmetricElement]->isSpareGate());
+                            if (usageIndex > 0) {
+                                generationInfo.addSpareUsageIndex(symmetricElement, usageIndex + offset * i);
+                                stateIndex += generationInfo.usageInfoBits();
+                            }
+                        }
+                    }
+
+                    // Store starting indices of symmetry groups
+                    std::vector<size_t> symmetryIndices;
+                    for (size_t i = 0; i < noSymmetricElements; ++i) {
+                        symmetryIndices.push_back(groupIndex + i * offset);
+                    }
+                    generationInfo.addSymmetry(offset, symmetryIndices);
+                }
+            }
+
+            // TODO symmetries in dependencies?
             // Consider dependencies
             for (size_t idDependency : getDependencies()) {
                 std::shared_ptr<DFTDependency<ValueType> const> dependency = getDependency(idDependency);
@@ -90,47 +147,61 @@ namespace storm {
                 visitQueue.push(dependency->dependentEvent()->id());
             }
             stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
-            
-            assert(stateIndex = mStateVectorSize);
+
+            // Visit all remaining states
+            for (size_t i = 0; i < visited.size(); ++i) {
+                if (!visited[i]) {
+                    visitQueue.push(i);
+                    stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
+                }
+            }
 
             STORM_LOG_TRACE(generationInfo);
+            assert(stateIndex == mStateVectorSize);
+            assert(visited.full());
             
             return generationInfo;
         }
         
         template<typename ValueType>
-        size_t DFT<ValueType>::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, std::set<size_t>& visited, size_t stateIndex) const {
+        size_t DFT<ValueType>::generateStateInfo(DFTStateGenerationInfo& generationInfo, size_t id, storm::storage::BitVector& visited, size_t stateIndex) const {
+            assert(!visited[id]);
+            visited.set(id);
+            
+            // Reserve bits for element
+            generationInfo.addStateIndex(id, stateIndex);
+            stateIndex += 2;
+            
+            if (isRepresentative(id)) {
+                generationInfo.addSpareActivationIndex(id, stateIndex);
+                ++stateIndex;
+            }
+            
+            if (mElements[id]->isSpareGate()) {
+                generationInfo.addSpareUsageIndex(id, stateIndex);
+                stateIndex += generationInfo.usageInfoBits();
+            }
+            
+            return stateIndex;
+        }
+        
+        template<typename ValueType>
+        size_t DFT<ValueType>::performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, storm::storage::BitVector& visited, size_t stateIndex) const {
             while (!visitQueue.empty()) {
                 size_t id = visitQueue.front();
                 visitQueue.pop();
-                if (visited.count(id) == 1) {
+                if (visited[id]) {
                     // Already visited
                     continue;
                 }
-                visited.insert(id);
-                DFTElementPointer element = mElements[id];
+                stateIndex = generateStateInfo(generationInfo, id, visited, stateIndex);
                 
                 // Insert children
-                if (element->isGate()) {
-                    for (auto const& child : std::static_pointer_cast<DFTGate<ValueType>>(element)->children()) {
+                if (mElements[id]->isGate()) {
+                    for (auto const& child : std::static_pointer_cast<DFTGate<ValueType>>(mElements[id])->children()) {
                         visitQueue.push(child->id());
                     }
                 }
-                
-                // Reserve bits for element
-                generationInfo.addStateIndex(id, stateIndex);
-                stateIndex += 2;
-                
-                if (isRepresentative(id)) {
-                    generationInfo.addSpareActivationIndex(id, stateIndex);
-                    ++stateIndex;
-                }
-                
-                if (element->isSpareGate()) {
-                    generationInfo.addSpareUsageIndex(id, stateIndex);
-                    stateIndex += generationInfo.usageInfoBits();
-                }
-                
             }
             return stateIndex;
         }
@@ -192,10 +263,10 @@ namespace storm {
                     stream << "\t** " << storm::storage::toChar(state->getElementState(elem->id()));
                     if(elem->isSpareGate()) {
                         size_t useId = state->uses(elem->id());
-                        if(state->isActive(useId)) {
-                            stream << " actively ";
+                        if(useId == elem->id() || state->isActive(useId)) {
+                            stream << "actively ";
                         }
-                        stream << " using " << useId;
+                        stream << "using " << useId;
                     }
                 }
                 stream << std::endl;
@@ -216,8 +287,8 @@ namespace storm {
                     if(elem->isSpareGate()) {
                         stream << "[";
                         size_t useId = state->uses(elem->id());
-                        if(state->isActive(useId)) {
-                            stream << " actively ";
+                        if(useId == elem->id() || state->isActive(useId)) {
+                            stream << "actively ";
                         }
                         stream << "using " << useId << "]";
                     }
diff --git a/src/storage/dft/DFT.h b/src/storage/dft/DFT.h
index a0059d81b..453db1443 100644
--- a/src/storage/dft/DFT.h
+++ b/src/storage/dft/DFT.h
@@ -10,6 +10,7 @@
 #include <boost/iterator/counting_iterator.hpp>
 
 #include "DFTElements.h"
+#include "elements/DFTRestriction.h"
 #include "../BitVector.h"
 #include "SymmetricUnits.h"
 #include "../../utility/math.h"
@@ -40,6 +41,7 @@ namespace storm {
             std::map<size_t, size_t> mSpareUsageIndex; // id spare -> index first bit in state
             std::map<size_t, size_t> mSpareActivationIndex; // id spare representative -> index in state
             std::vector<size_t> mIdToStateIndex; // id -> index first bit in state
+            std::vector<std::pair<size_t, std::vector<size_t>>> mSymmetries; // pair (lenght of symmetry group, vector indicating the starting points of the symmetry groups)
             
         public:
             
@@ -78,6 +80,24 @@ namespace storm {
                 return mSpareActivationIndex.at(id);
             }
             
+            size_t addSymmetry(size_t lenght, std::vector<size_t>& startingIndices) {
+                mSymmetries.push_back(std::make_pair(lenght, startingIndices));
+            }
+            
+            size_t getSymmetrySize() const {
+                return mSymmetries.size();
+            }
+            
+            size_t getSymmetryLength(size_t pos) const {
+                assert(pos < mSymmetries.size());
+                return mSymmetries[pos].first;
+            }
+            
+            std::vector<size_t> const& getSymmetryIndices(size_t pos) const {
+                assert(pos < mSymmetries.size());
+                return mSymmetries[pos].second;
+            }
+            
             friend std::ostream& operator<<(std::ostream& os, DFTStateGenerationInfo const& info) {
                 os << "Id to state index:" << std::endl;
                 for (size_t id = 0; id < info.mIdToStateIndex.size(); ++id) {
@@ -91,6 +111,14 @@ namespace storm {
                 for (auto pair : info.mSpareActivationIndex) {
                     os << pair.first << " -> " << pair.second << std::endl;
                 }
+                os << "Symmetries:" << std::endl;
+                for (auto pair : info.mSymmetries) {
+                    os << "Length: " << pair.first << ", starting indices: ";
+                    for (size_t index : pair.second) {
+                        os << index << ", ";
+                    }
+                    os << std::endl;
+                }
                 return os;
             }
             
@@ -125,9 +153,11 @@ namespace storm {
         public:
             DFT(DFTElementVector const& elements, DFTElementPointer const& tle);
             
-            DFTStateGenerationInfo buildStateGenerationInfo(std::vector<size_t> const& subTreeRoots, std::vector<std::vector<size_t>> const& symmetries) const;
+            DFTStateGenerationInfo buildStateGenerationInfo(storm::storage::DFTIndependentSymmetries const& symmetries) const;
+            
+            size_t generateStateInfo(DFTStateGenerationInfo& generationInfo, size_t id, storm::storage::BitVector& visited, size_t stateIndex) const;
 
-            size_t performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, std::set<size_t>& visited, size_t stateIndex) const;
+            size_t performStateGenerationInfoDFS(DFTStateGenerationInfo& generationInfo, std::queue<size_t>& visitQueue, storm::storage::BitVector& visited, size_t stateIndex) const;
             
             size_t stateVectorSize() const {
                 return mStateVectorSize;
diff --git a/src/storage/dft/DFTBuilder.cpp b/src/storage/dft/DFTBuilder.cpp
index 1979aea3b..213fce0df 100644
--- a/src/storage/dft/DFTBuilder.cpp
+++ b/src/storage/dft/DFTBuilder.cpp
@@ -27,6 +27,7 @@ namespace storm {
                         childElement->addParent(gate);
                     } else {
                         // Child not found -> find first dependent event to assure that child is dependency
+                        // TODO: Not sure whether this is the intended behaviour?
                         auto itFind = mElements.find(child + "_1");
                         assert(itFind != mElements.end());
                         assert(itFind->second->isDependency());
@@ -34,6 +35,17 @@ namespace storm {
                     }
                 }
             }
+            
+            for(auto& elem : mRestrictionChildNames) {
+                for(auto const& childName : elem.second) {
+                    auto itFind = mElements.find(childName);
+                    assert(itFind != mElements.end());
+                    DFTElementPointer childElement = itFind->second;
+                    assert(!childElement->isDependency() && !childElement->isRestriction());
+                    elem.first->pushBackChild(childElement);
+                    childElement->addRestriction(elem.first);
+                }
+            }
 
             // Initialize dependencies
             for (auto& dependency : mDependencies) {
@@ -45,6 +57,17 @@ namespace storm {
                 dependentEvent->addIngoingDependency(dependency);
             }
 
+//            for (auto& restriction : mRestrictions) {
+//                std::set<DFTGatePointer> parentsOfRestrictedElements;
+//                for (auto& child : restriction->children()) {
+//                    for(DFTGatePointer& parent : child->parents()) {
+//                        parentsOfRestrictedElements.insert(parent);
+//                    }
+//                }
+//
+//
+//            }
+
             // Sort elements topologically
             // compute rank
             for (auto& elem : mElements) {
@@ -82,6 +105,34 @@ namespace storm {
             return elem->rank();
         }
 
+        template<typename ValueType>
+        bool DFTBuilder<ValueType>::addRestriction(std::string const& name, std::vector<std::string> const& children, DFTElementType tp) {
+            if (children.size() <= 1) {
+                STORM_LOG_ERROR("Sequence enforcers require at least two children");
+            }
+            if (mElements.count(name) != 0) {
+                return false;
+            }
+            DFTRestrictionPointer restr;
+            switch (tp) {
+                case DFTElementType::SEQ:
+                    restr = std::make_shared<DFTSeq < ValueType>>
+                    (mNextId++, name);
+                    break;
+                case DFTElementType::MUTEX:
+                    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Gate type not supported.");
+                    break;
+                default:
+                    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Gate type not known.");
+                    break;
+            }
+
+            mElements[name] = restr;
+            mRestrictionChildNames[restr] = children;
+            mRestrictions.push_back(restr);
+            return true;
+        }
+
         template<typename ValueType>
         bool DFTBuilder<ValueType>::addStandardGate(std::string const& name, std::vector<std::string> const& children, DFTElementType tp) {
             assert(children.size() > 0);
@@ -108,11 +159,11 @@ namespace storm {
                    break;
                 case DFTElementType::BE:
                 case DFTElementType::VOT:
+                case DFTElementType::PDEP:
                     // Handled separately
                     STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Gate type handled separately.");
                 case DFTElementType::CONSTF:
                 case DFTElementType::CONSTS:
-                case DFTElementType::PDEP:
                     STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Gate type not supported.");
                 default:
                     STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Gate type not known.");
diff --git a/src/storage/dft/DFTBuilder.h b/src/storage/dft/DFTBuilder.h
index 17a16cb49..8ef0d6cf1 100644
--- a/src/storage/dft/DFTBuilder.h
+++ b/src/storage/dft/DFTBuilder.h
@@ -3,6 +3,7 @@
 #define	DFTBUILDER_H
 
 #include "DFTElements.h"
+#include "elements/DFTRestriction.h"
 #include <iostream>
 #include <unordered_map>
 #include <map>
@@ -21,13 +22,16 @@ namespace storm {
             using DFTGatePointer = std::shared_ptr<DFTGate<ValueType>>;
             using DFTGateVector = std::vector<DFTGatePointer>;
             using DFTDependencyPointer = std::shared_ptr<DFTDependency<ValueType>>;
+            using DFTRestrictionPointer = std::shared_ptr<DFTRestriction<ValueType>>;
 
         private:
             std::size_t mNextId = 0;
             std::string topLevelIdentifier;
             std::unordered_map<std::string, DFTElementPointer> mElements;
             std::unordered_map<DFTElementPointer, std::vector<std::string>> mChildNames;
+            std::unordered_map<DFTRestrictionPointer, std::vector<std::string>> mRestrictionChildNames;
             std::vector<DFTDependencyPointer> mDependencies;
+            std::vector<DFTRestrictionPointer> mRestrictions;
             
         public:
             DFTBuilder() {
@@ -53,9 +57,19 @@ namespace storm {
             bool addSpareElement(std::string const& name, std::vector<std::string> const& children) {
                 return addStandardGate(name, children, DFTElementType::SPARE);
             }
+
+            bool addSequenceEnforcer(std::string const& name, std::vector<std::string> const& children) {
+                return addRestriction(name, children, DFTElementType::SEQ);
+            }
+
+            bool addMutex(std::string const& name, std::vector<std::string> const& children) {
+                return addRestriction(name, children, DFTElementType::MUTEX);
+            }
             
             bool addDepElement(std::string const& name, std::vector<std::string> const& children, ValueType probability) {
-                assert(children.size() > 1);
+                if(children.size() <= 1) {
+                    STORM_LOG_ERROR("Dependencies require at least two children");
+                }
                 if(mElements.count(name) != 0) {
                     // Element with that name already exists.
                     return false;
@@ -138,7 +152,8 @@ namespace storm {
             unsigned computeRank(DFTElementPointer const& elem);
             
             bool addStandardGate(std::string const& name, std::vector<std::string> const& children, DFTElementType tp);
-            
+            bool addRestriction(std::string const& name, std::vector<std::string> const& children, DFTElementType tp);
+
             enum class topoSortColour {WHITE, BLACK, GREY}; 
             
             void topoVisit(DFTElementPointer const& n, std::map<DFTElementPointer, topoSortColour, OrderElementsById<ValueType>>& visited, DFTElementVector& L);
diff --git a/src/storage/dft/DFTElements.cpp b/src/storage/dft/DFTElements.cpp
index e06b3fe7d..7bb8c2cf5 100644
--- a/src/storage/dft/DFTElements.cpp
+++ b/src/storage/dft/DFTElements.cpp
@@ -10,6 +10,7 @@ namespace storm {
             if (state.dontCare(mId)) {
                 return false;
             }
+            
             // Check that no outgoing dependencies can be triggered anymore
             for (DFTDependencyPointer dependency : mOutgoingDependencies) {
                 if (state.isOperational(dependency->dependentEvent()->id()) && state.isOperational(dependency->triggerEvent()->id())) {
@@ -22,6 +23,9 @@ namespace storm {
                     return false;
                 }
             }
+            
+            
+        
 
             state.setDontCare(mId);
             return true;
@@ -52,7 +56,7 @@ namespace storm {
 
         template<typename ValueType>
         std::vector<size_t> DFTElement<ValueType>::independentSubDft() const {
-            std::cout << "INDEPENDENT SUBTREE CALL " << this->id() << std::endl;
+            //std::cout << "INDEPENDENT SUBTREE CALL " << this->id() << std::endl;
             std::vector<size_t> res;
             res.push_back(this->id());
             return res;
diff --git a/src/storage/dft/DFTElements.h b/src/storage/dft/DFTElements.h
index 22c819f84..3236a5551 100644
--- a/src/storage/dft/DFTElements.h
+++ b/src/storage/dft/DFTElements.h
@@ -26,6 +26,8 @@ namespace storm {
         
         template<typename ValueType>
         class DFTDependency;
+        template<typename ValueType>
+        class DFTRestriction;
 
         template<typename ValueType>
         class DFTElement {
@@ -34,6 +36,9 @@ namespace storm {
             using DFTGateVector = std::vector<DFTGatePointer>;
             using DFTDependencyPointer = std::shared_ptr<DFTDependency<ValueType>>;
             using DFTDependencyVector = std::vector<DFTDependencyPointer>;
+            using DFTRestrictionPointer = std::shared_ptr<DFTRestriction<ValueType>>;
+            using DFTRestrictionVector = std::vector<DFTRestrictionPointer>;
+
 
         protected:
             size_t mId;
@@ -41,6 +46,8 @@ namespace storm {
             size_t mRank = -1;
             DFTGateVector mParents;
             DFTDependencyVector mOutgoingDependencies;
+            DFTRestrictionVector mRestrictions;
+
 
         public:
             DFTElement(size_t id, std::string const& name) :
@@ -95,6 +102,10 @@ namespace storm {
             virtual bool isDependency() const {
                 return false;
             }
+            
+            virtual bool isRestriction() const {
+                return false;
+            }
 
             virtual void setId(size_t newId) {
                 mId = newId;
@@ -118,6 +129,15 @@ namespace storm {
                 }
             }
 
+            bool addRestriction(DFTRestrictionPointer const& e) {
+                if (std::find(mRestrictions.begin(), mRestrictions.end(), e) != mRestrictions.end()) {
+                    return false;
+                } else {
+                    mRestrictions.push_back(e);
+                    return true;
+                }
+            }
+
              bool hasOnlyStaticParents() const {
                 for(auto const& parent : mParents) {
                     if(!isStaticGateType(parent->type())) {
@@ -139,6 +159,18 @@ namespace storm {
             DFTGateVector const& parents() const {
                 return mParents;
             }
+            
+            bool hasRestrictions() const {
+                return !mRestrictions.empty();
+            }
+            
+            size_t nrRestrictions() const {
+                return mRestrictions.size();
+            }
+            
+            DFTRestrictionVector const& restrictions() const {
+                return mRestrictions;
+            }
 
             std::vector<size_t> parentIds() const {
                 std::vector<size_t> res;
@@ -351,7 +383,11 @@ namespace storm {
                         queues.propagateFailure(parent);
                     }
                 }
+                for(std::shared_ptr<DFTRestriction<ValueType>> restr : this->mRestrictions) {
+                    queues.checkRestrictionLater(restr);
+                }
                 state.setFailed(this->mId);
+                // TODO can this be moved towards DFTSpare?
                 if (this->isSpareGate()) {
                     this->finalizeSpare(state);
                 }
@@ -376,12 +412,12 @@ namespace storm {
             }
 
             /**
-             * Finish failed/failsafe spare gate by activating the children and setting the useIndex to zero.
+             * Finish failed/failsafe spare gate by activating the children and setting the useIndex to the spare id.
              * This prevents multiple fail states with different usages or activations.
              * @param state The current state.
              */
             void finalizeSpare(DFTState<ValueType>& state) const {
-                state.setUses(this->mId, 0);
+                state.setUses(this->mId, this->mId);
                 for (auto child : this->children()) {
                     if (!state.isActive(child->id())) {
                         state.activate(child->id());
diff --git a/src/storage/dft/DFTState.cpp b/src/storage/dft/DFTState.cpp
index 7f34e03ef..be5c33060 100644
--- a/src/storage/dft/DFTState.cpp
+++ b/src/storage/dft/DFTState.cpp
@@ -242,6 +242,37 @@ namespace storm {
             }
             return false;
         }
+        
+        template<typename ValueType>
+        bool DFTState<ValueType>::orderBySymmetry() {
+            bool changed = false;
+            for (size_t pos = 0; pos < mStateGenerationInfo.getSymmetrySize(); ++pos) {
+                // Check each symmetry
+                size_t length = mStateGenerationInfo.getSymmetryLength(pos);
+                std::vector<size_t> symmetryIndices = mStateGenerationInfo.getSymmetryIndices(pos);
+                // Sort symmetry group in decreasing order by bubble sort
+                // TODO use better algorithm?
+                size_t tmp, elem1, elem2;
+                size_t n = symmetryIndices.size();
+                do {
+                    tmp = 0;
+                    for (size_t i = 1; i < n; ++i) {
+                        elem1 = mStatus.getAsInt(symmetryIndices[i-1], length);
+                        elem2 = mStatus.getAsInt(symmetryIndices[i], length);
+                        if (elem1 < elem2) {
+                            // Swap elements
+                            mStatus.setFromInt(symmetryIndices[i-1], length, elem2);
+                            mStatus.setFromInt(symmetryIndices[i], length, elem1);
+                            tmp = i;
+                            changed = true;
+                        }
+                    }
+                    n = tmp;
+                } while (n > 0);
+            }
+            return changed;
+        }
+
 
         // Explicitly instantiate the class.
         template class DFTState<double>;
diff --git a/src/storage/dft/DFTState.h b/src/storage/dft/DFTState.h
index 7a26e10ab..eccf7f01d 100644
--- a/src/storage/dft/DFTState.h
+++ b/src/storage/dft/DFTState.h
@@ -164,6 +164,12 @@ namespace storm {
              */
             void letDependencyBeUnsuccessful(size_t index = 0);
             
+            /**
+             * Order the state in decreasing order using the symmetries.
+             * @return True, if elements were swapped, false if nothing changed.
+             */
+            bool orderBySymmetry();
+            
             std::string getCurrentlyFailableString() const {
                 std::stringstream stream;
                 if (nrFailableDependencies() > 0) {
diff --git a/src/storage/dft/DFTStateSpaceGenerationQueues.h b/src/storage/dft/DFTStateSpaceGenerationQueues.h
index b881a220b..6ded7cb38 100644
--- a/src/storage/dft/DFTStateSpaceGenerationQueues.h
+++ b/src/storage/dft/DFTStateSpaceGenerationQueues.h
@@ -15,6 +15,8 @@ namespace storm {
         class DFTGate;
         template<typename ValueType>
         class DFTElement;
+        template<typename ValueType>
+        class DFTRestriction;
 
 
         template<typename ValueType>
@@ -24,11 +26,14 @@ namespace storm {
             using DFTElementVector = std::vector<DFTElementPointer>;
             using DFTGatePointer = std::shared_ptr<DFTGate<ValueType>>;
             using DFTGateVector = std::vector<DFTGatePointer>;
+            using DFTRestrictionPointer = std::shared_ptr<DFTRestriction<ValueType>>;
+            using DFTRestrictionVector = std::vector<DFTRestrictionPointer>;
 
             std::priority_queue<DFTGatePointer, DFTGateVector, OrderElementsByRank<ValueType>> failurePropagation;
             DFTGateVector failsafePropagation;
             DFTElementVector dontcarePropagation;
             DFTElementVector activatePropagation;
+            DFTRestrictionVector restrictionChecks;
 
         public:
             void propagateFailure(DFTGatePointer const& elem) {
@@ -45,6 +50,21 @@ namespace storm {
                 return next;
             }
             
+            bool restrictionChecksDone() const {
+                return restrictionChecks.empty();
+            }
+            
+            DFTRestrictionPointer nextRestrictionCheck() {
+                assert(!restrictionChecksDone());
+                DFTRestrictionPointer next = restrictionChecks.back();
+                restrictionChecks.pop_back();
+                return next;
+            }
+            
+            void checkRestrictionLater(DFTRestrictionPointer const& restr) {
+                restrictionChecks.push_back(restr);
+            }
+            
             bool failsafePropagationDone() const {
                 return failsafePropagation.empty();
             }
diff --git a/src/storage/dft/elements/DFTRestriction.h b/src/storage/dft/elements/DFTRestriction.h
index 8bd586deb..46b6e23cc 100644
--- a/src/storage/dft/elements/DFTRestriction.h
+++ b/src/storage/dft/elements/DFTRestriction.h
@@ -6,8 +6,195 @@ namespace storm {
     namespace storage {
         template<typename ValueType>
         class DFTRestriction : public DFTElement<ValueType> {
+            using DFTElementPointer = std::shared_ptr<DFTElement<ValueType>>;
+            using DFTElementVector = std::vector<DFTElementPointer>;
+        protected:
+            DFTElementVector mChildren;
+
+        public:
+            DFTRestriction(size_t id, std::string const& name, DFTElementVector const& children) :
+            DFTElement<ValueType>(id, name), mChildren(children)
+                    {}
+
+            virtual ~DFTRestriction() {}
+
+            void pushBackChild(DFTElementPointer elem) {
+                return mChildren.push_back(elem);
+            }
+
+            size_t nrChildren() const override {
+                return mChildren.size();
+            }
+
+            DFTElementVector const& children() const {
+                return mChildren;
+            }
+
+            virtual bool isRestriction() const override {
+                return true;
+            }
+
+
+            virtual std::string typestring() const = 0;
+
+            virtual void checkFails(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const = 0;
+
+            virtual void checkFailsafe(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const = 0;
+
+            virtual void extendSpareModule(std::set<size_t>& elementsInSpareModule) const override {
+                DFTElement<ValueType>::extendSpareModule(elementsInSpareModule);
+                for(auto const& child : mChildren) {
+                    if(elementsInSpareModule.count(child->id()) == 0) {
+                        elementsInSpareModule.insert(child->id());
+                        child->extendSpareModule(elementsInSpareModule);
+                    }
+                }
+            }
+
+            virtual std::vector<size_t> independentUnit() const override {
+                std::set<size_t> unit = {this->mId};
+                for(auto const& child : mChildren) {
+                    child->extendUnit(unit);
+                }
+                return std::vector<size_t>(unit.begin(), unit.end());
+            }
+
+            virtual void extendUnit(std::set<size_t>& unit) const override {
+                DFTElement<ValueType>::extendUnit(unit);
+                for(auto const& child : mChildren) {
+                    child->extendUnit(unit);
+                }
+            }
+
+            virtual std::vector<size_t> independentSubDft() const override {
+                auto prelRes = DFTElement<ValueType>::independentSubDft();
+                if(prelRes.empty()) {
+                    // No elements (especially not this->id) in the prelimanry result, so we know already that it is not a subdft.
+                    return prelRes;
+                }
+                std::set<size_t> unit(prelRes.begin(), prelRes.end());
+                std::vector<size_t> pids = this->parentIds();
+                for(auto const& child : mChildren) {
+                    child->extendSubDft(unit, pids);
+                    if(unit.empty()) {
+                        // Parent in the subdft, ie it is *not* a subdft
+                        break;
+                    }
+                }
+                return std::vector<size_t>(unit.begin(), unit.end());
+            }
+
+            virtual void extendSubDft(std::set<size_t>& elemsInSubtree, std::vector<size_t> const& parentsOfSubRoot) const override {
+                if(elemsInSubtree.count(this->id()) > 0) return;
+                DFTElement<ValueType>::extendSubDft(elemsInSubtree, parentsOfSubRoot);
+                if(elemsInSubtree.empty()) {
+                    // Parent in the subdft, ie it is *not* a subdft
+                    return;
+                }
+                for(auto const& child : mChildren) {
+                    child->extendSubDft(elemsInSubtree, parentsOfSubRoot);
+                    if(elemsInSubtree.empty()) {
+                        // Parent in the subdft, ie it is *not* a subdft
+                        return;
+                    }
+                }
+            }
+
+
+            virtual std::string toString() const override {
+                std::stringstream stream;
+                stream << "{" << this->name() << "} " << typestring() << "( ";
+                typename DFTElementVector::const_iterator it = mChildren.begin();
+                stream << (*it)->name();
+                ++it;
+                while(it != mChildren.end()) {
+                    stream <<  ", " << (*it)->name();
+                    ++it;
+                }
+                stream << ")";
+                return stream.str();
+            }
+
+            virtual bool checkDontCareAnymore(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const override {
+                return false;
+            }
+
+
+        protected:
+
+            void fail(DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const {
+                state.markAsInvalid();
+            }
+
+            void failsafe(DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const {
+
+            }
+
+            bool hasFailsafeChild(DFTState<ValueType>& state) const {
+                for(auto const& child : mChildren) {
+                    if(state.isFailsafe(child->id()))
+                    {
+                        return true;
+                    }
+                }
+                return false;
+            }
+
+            bool hasFailedChild(DFTState<ValueType>& state) const {
+                for(auto const& child : mChildren) {
+                    if(state.hasFailed(child->id())) {
+                        return true;
+                    }
+                }
+                return false;
+            }
+
+
+
+
+        };
+
+        template<typename ValueType>
+        class DFTSeq : public DFTRestriction<ValueType> {
+
+
+        public:
+            DFTSeq(size_t id, std::string const& name, std::vector<std::shared_ptr<DFTElement<ValueType>>> const& children = {}) :
+                    DFTRestriction<ValueType>(id, name, children)
+            {}
+
+            void checkFails(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const override {
+                assert(queues.failurePropagationDone());
+                bool childOperationalBefore = false;
+                for(auto const& child : this->mChildren)
+                {
+                    if(!state.hasFailed(child->id())) {
+                        childOperationalBefore = true;
+                    } else if(childOperationalBefore && state.hasFailed(child->id())){
+                        this->fail(state, queues);
+                        return;
+                    }
+                }
+
+            }
+
+            void checkFailsafe(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) const override {
+                
+
+            }
             
+            bool checkDontCareAnymore(storm::storage::DFTState<ValueType>& state, DFTStateSpaceGenerationQueues<ValueType>& queues) {
+                
+            }
+
+            virtual DFTElementType type() const override {
+                return DFTElementType::SEQ;
+            }
+
+            std::string typestring() const override {
+                return "SEQ";
+            }
         };
-        
+
     }
 }
\ No newline at end of file
diff --git a/src/storm-dyftee.cpp b/src/storm-dyftee.cpp
index 2160ed210..ba50d90f3 100644
--- a/src/storm-dyftee.cpp
+++ b/src/storm-dyftee.cpp
@@ -28,17 +28,19 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false)
     std::vector<std::shared_ptr<storm::logic::Formula>> parsedFormulas = storm::parseFormulasForExplicit(property);
     std::vector<std::shared_ptr<const storm::logic::Formula>> formulas(parsedFormulas.begin(), parsedFormulas.end());
     assert(formulas.size() == 1);
-    std::cout << "parsed formula." << std::endl;
+    std::cout << "Parsed formula." << std::endl;
+    std::map<size_t, std::vector<std::vector<size_t>>> emptySymmetry;
+    storm::storage::DFTIndependentSymmetries symmetries(emptySymmetry);
     if(symred) {
         std::cout << dft.getElementsString() << std::endl;
         auto colouring = dft.colourDFT();
-        storm::storage::DFTIndependentSymmetries symmetries = dft.findSymmetries(colouring);
-        std::cout << symmetries;
+        symmetries = dft.findSymmetries(colouring);
+        std::cout << "Symmetries: " << symmetries << std::endl;
     }
         
     // Building Markov Automaton
     std::cout << "Building Model..." << std::endl;
-    storm::builder::ExplicitDFTModelBuilder<ValueType> builder(dft);
+    storm::builder::ExplicitDFTModelBuilder<ValueType> builder(dft, symmetries);
     typename storm::builder::ExplicitDFTModelBuilder<ValueType>::LabelOptions labeloptions; // TODO initialize this with the formula
     std::shared_ptr<storm::models::sparse::Model<ValueType>> model = builder.buildModel(labeloptions);
     std::cout << "Built Model" << std::endl;
@@ -48,11 +50,10 @@ void analyzeDFT(std::string filename, std::string property, bool symred = false)
     
     std::cout << "Bisimulation..." << std::endl;
     
-    if (model->isOfType(storm::models::ModelType::Ctmc)) {
+    if (model->getNumberOfStates() > 500 && model->isOfType(storm::models::ModelType::Ctmc)) {
             model =  storm::performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas, storm::storage::BisimulationType::Weak)->template as<storm::models::sparse::Ctmc<ValueType>>();
     }
     
-    
     model->printModelInformationToStream(std::cout);
     
     // Model checking