diff --git a/src/storage/dft/DFT.cpp b/src/storage/dft/DFT.cpp
index e9b4c659e..dd771f230 100644
--- a/src/storage/dft/DFT.cpp
+++ b/src/storage/dft/DFT.cpp
@@ -110,22 +110,35 @@ namespace storm {
                 visitQueue.push(mTopLevelIndex);
                 stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
             } else {
-                for (auto const& symmetryGroup : symmetries.groups) {
-                    assert(!symmetryGroup.second.empty());
+                for (size_t symmetryIndex : symmetries.sortedSymmetries) {
+                    assert(!visited[symmetryIndex]);
+                    auto const& symmetryGroup = symmetries.groups.at(symmetryIndex);
+                    assert(!symmetryGroup.empty());
 
                     // Insert all elements of first subtree of each symmetry
                     size_t groupIndex = stateIndex;
-                    for (std::vector<size_t> const& symmetryElement : symmetryGroup.second) {
-                        stateIndex = generateStateInfo(generationInfo, symmetryElement[0], visited, stateIndex);
+                    for (std::vector<size_t> const& symmetryElement : symmetryGroup) {
+                        if (visited[symmetryElement[0]]) {
+                            groupIndex = std::min(groupIndex, generationInfo.getStateIndex(symmetryElement[0]));
+                        } else {
+                            stateIndex = generateStateInfo(generationInfo, symmetryElement[0], visited, stateIndex);
+                        }
                     }
                     size_t offset = stateIndex - groupIndex;
                     
                     // Mirror symmetries
-                    size_t noSymmetricElements = symmetryGroup.second.front().size();
+                    size_t noSymmetricElements = symmetryGroup.front().size();
                     assert(noSymmetricElements > 1);
 
-                    for (std::vector<size_t> symmetricElements : symmetryGroup.second) {
+                    for (std::vector<size_t> symmetricElements : symmetryGroup) {
                         assert(symmetricElements.size() == noSymmetricElements);
+                        if (visited[symmetricElements[1]]) {
+                            // Elements already mirrored
+                            for (size_t index : symmetricElements) {
+                                assert(visited[index]);
+                            }
+                            continue;
+                        }
 
                         // Initialize for original element
                         size_t originalElement = symmetricElements[0];
@@ -181,6 +194,8 @@ namespace storm {
                     stateIndex = performStateGenerationInfoDFS(generationInfo, visitQueue, visited, stateIndex);
                 }
             }
+            
+            generationInfo.generateSymmetries(symmetries);
 
             STORM_LOG_TRACE(generationInfo);
             assert(stateIndex == mStateVectorSize);
diff --git a/src/storage/dft/DFTStateGenerationInfo.h b/src/storage/dft/DFTStateGenerationInfo.h
index 4b3ef92c9..d4bbdee8b 100644
--- a/src/storage/dft/DFTStateGenerationInfo.h
+++ b/src/storage/dft/DFTStateGenerationInfo.h
@@ -10,7 +10,7 @@ namespace storm {
             std::vector<size_t> mIdToStateIndex; // id -> index first bit in state
             std::map<size_t, std::vector<size_t>> mSeqRestrictionPreElements; // id -> list of restriction pre elements;
             std::map<size_t, std::vector<size_t>> mSeqRestrictionPostElements; // id -> list of restriction post elements;
-            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)
+            std::vector<std::pair<size_t, std::vector<size_t>>> mSymmetries; // pair (length of symmetry group, vector indicating the starting points of the symmetry groups)
 
         public:
 
@@ -70,6 +70,32 @@ namespace storm {
             void addSymmetry(size_t length, std::vector<size_t>& startingIndices) {
                 mSymmetries.push_back(std::make_pair(length, startingIndices));
             }
+            
+            /**
+             * Generate more symmetries by combining two symmetries
+             */
+            void generateSymmetries(storm::storage::DFTIndependentSymmetries const& symmetries) {
+                for (size_t i = 0; i < mSymmetries.size(); ++i) {
+                    size_t childStart = mSymmetries[i].second[0];
+                    size_t childLength = mSymmetries[i].first;
+                    for (size_t j = i+1; j < mSymmetries.size(); ++j) {
+                        size_t parentStart = mSymmetries[j].second[0];
+                        size_t parentLength = mSymmetries[j].first;
+                        // Check if child lies in parent
+                        if (parentStart <= childStart && childStart + childLength < parentStart + parentLength) {
+                            // Get symmetric start by applying the bijection
+                            std::vector<size_t> newStarts;
+                            for (size_t symmetryStarts : mSymmetries[i].second) {
+                                newStarts.push_back(symmetryStarts + parentLength);
+                            }
+                            // Insert after child
+                            mSymmetries.insert(mSymmetries.begin() + i + 1, std::make_pair(childLength, newStarts));
+                            ++i;
+                            break;
+                        }
+                    }
+                }
+            }
 
             size_t getSymmetrySize() const {
                 return mSymmetries.size();
diff --git a/src/storage/dft/SymmetricUnits.h b/src/storage/dft/SymmetricUnits.h
index 58a70e492..33c8efeae 100644
--- a/src/storage/dft/SymmetricUnits.h
+++ b/src/storage/dft/SymmetricUnits.h
@@ -1,19 +1,123 @@
 #pragma once
 
+#include "src/utility/macros.h"
+
 namespace storm {
     namespace storage {
         struct DFTIndependentSymmetries {
             std::map<size_t, std::vector<std::vector<size_t>>> groups;
             
+            std::vector<size_t> sortedSymmetries;
+            
+            bool existsInFirstSymmetry(size_t index, size_t value) {
+                for (std::vector<size_t> symmetry : groups[index]) {
+                    if (symmetry.front() == value) {
+                        return true;
+                    }
+                }
+                return false;
+            }
+            
+            bool existsInSymmetry(size_t index, size_t value) {
+                for (std::vector<size_t> symmetry : groups[index]) {
+                    for (size_t index : symmetry) {
+                        if (index == value) {
+                            return true;
+                        }
+                    }
+                }
+                return false;
+            }
+            
+            /**
+             * Apply symmetry and get bijection. Result is symmetry(value)[index]
+             *
+             * @param symmetry The symmetry.
+             * @param value The value to apply the bijection to.
+             * @param index The index of the symmetry group to apply.
+             * @return Symmetric value, -1 if the bijection could not be applied.
+             */
+            int applySymmetry(std::vector<std::vector<size_t>> const& symmetry, size_t value, size_t index) const {
+                for (std::vector<size_t> element : symmetry) {
+                    if (element[0] == value) {
+                        return element[index];
+                    }
+                }
+                return -1;
+            }
+            
+            std::vector<std::vector<size_t>> createSymmetry(std::vector<std::vector<size_t>> parentSymmetry, std::vector<std::vector<size_t>> childSymmetry, size_t index) {
+                std::vector<std::vector<size_t>> result;
+                for (std::vector<size_t> childSymmetry : childSymmetry) {
+                    std::vector<size_t> symmetry;
+                    for (size_t child : childSymmetry) {
+                        int bijectionValue = applySymmetry(parentSymmetry, child, index);
+                        if (bijectionValue >= 0) {
+                            symmetry.push_back(bijectionValue);
+                        } else {
+                            return result;
+                        }
+                    }
+                    result.push_back(symmetry);
+                }
+                return result;
+            }
+            
+            void sortHierarchical(size_t parent, std::vector<size_t>& candidates) {
+                std::vector<size_t> children;
+                for (int i = candidates.size() - 1; i >= 0; --i) {
+                    size_t currentRoot = candidates[i];
+                    if (existsInSymmetry(parent, currentRoot)) {
+                        // Is child
+                        STORM_LOG_TRACE(currentRoot << " is child of " << parent);
+                        children.insert(children.begin(), currentRoot);
+                        candidates.erase(candidates.begin()+i);
+                    }
+                }
+                for (size_t i = 0; i < children.size(); ++i) {
+                    std::vector<std::vector<size_t>> possibleSymmetry = createSymmetry(groups.at(parent), groups.at(children[i]), 1);
+                    for (size_t j = i + 1; j < children.size(); ++j) {
+                        if (possibleSymmetry == groups.at(children[j])) {
+                            STORM_LOG_TRACE("Child " << children[j] << " ignored as created by symmetries " << parent << " and " << children[j]);
+                            groups.erase(children[j]);
+                            children.erase(children.begin() + j);
+                        }
+                    }
+                }
+                while (!children.empty()) {
+                    // Insert largest element
+                    size_t largestChild = children.back();
+                    children.pop_back();
+                    sortHierarchical(largestChild, children);
+                    sortedSymmetries.push_back(largestChild);
+                }
+            }
+            
             DFTIndependentSymmetries(std::map<size_t, std::vector<std::vector<size_t>>> groups) : groups(groups) {
+                std::vector<size_t> sortedGroups;
+                for (auto const& cl : groups) {
+                    sortedGroups.push_back(cl.first);
+                }
+                // Sort by length of symmetry
+                std::sort(sortedGroups.begin(), sortedGroups.end(), [&](const size_t left, const size_t right) {
+                    return groups.at(left).size() < groups.at(right).size();
+                });
                 
+                // Sort hierarchical
+                while (!sortedGroups.empty()) {
+                    // Insert largest element
+                    size_t currentRoot = sortedGroups.back();
+                    sortedGroups.pop_back();
+                    sortHierarchical(currentRoot, sortedGroups);
+                    sortedSymmetries.push_back(currentRoot);
+                }
             }
         };
         
         inline std::ostream& operator<<(std::ostream& os, DFTIndependentSymmetries const& s)  {
-            for(auto const& cl : s.groups) {
-                os << "Symmetry group for " << cl.first << std::endl;
-                for(auto const& eqClass : cl.second) {
+            for(size_t index : s.sortedSymmetries) {
+                os << "Symmetry group for " << index << std::endl;
+                for(auto const& eqClass : s.groups.at(index)) {
                     for(auto const& i : eqClass) {
                         os << i << " ";
                     }