diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3997cde5..6cdf2ffc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,17 @@
 # Release history
 
+## Version 1.1.3 [unpublished]
+
+### Changes
+
+- Implemented the DRSC algorithm for SumEquals constraints. A side-effect of this implementation is that you can now
+  ask for multiple points with opt.ask() while using this type of constraint (previously you could only ask for one).
+- Dropped support for Python 3.9
+
+### Bugfixes
+
+- Fixed TypeError when combining categorical and numerical dimensions in get_Brownie_Bee_Pareto
+
 ## Version 1.1.2 [unpublished]
 
 ### Changes
diff --git a/ProcessOptimizer/optimizer/optimizer.py b/ProcessOptimizer/optimizer/optimizer.py
index 9dd8e5ce..229e1bf9 100644
--- a/ProcessOptimizer/optimizer/optimizer.py
+++ b/ProcessOptimizer/optimizer/optimizer.py
@@ -269,9 +269,11 @@ def __init__(
             )
 
         # check if regressor
-        if not is_regressor(base_estimator) and base_estimator is not None:
-            raise ValueError("%s has to be a regressor." % base_estimator)
-
+        try:
+            if base_estimator is not None and not is_regressor(base_estimator):
+                raise ValueError("%s has to be a regressor." % base_estimator)
+        except AttributeError as me:
+            raise ValueError("%s has to be a regressor." % base_estimator) from me
         # treat per second acqusition function specially
         is_multi_regressor = isinstance(base_estimator, MultiOutputRegressor)
         if "ps" in self.acq_func and not is_multi_regressor:
@@ -1309,12 +1311,12 @@ def __MinimalDistance(X, Y):
     # This function calls NSGAII to estimate the Pareto Front
     def NSGAII(self, MU=40):
         from ._NSGA2 import NSGAII
-        
+
         if MU % 4 != 0:
             raise ValueError(
                 "Number of simulated points must be divisible by 4 for the NSGAII algorithm"
             )
-            
+
         pop, logbook, front = NSGAII(
             self.n_objectives,
             self.__ObjectiveGP,
diff --git a/ProcessOptimizer/space/DRSC.py b/ProcessOptimizer/space/DRSC.py
new file mode 100644
index 00000000..adc2b893
--- /dev/null
+++ b/ProcessOptimizer/space/DRSC.py
@@ -0,0 +1,412 @@
+from typing import List, Tuple, Callable, Optional
+
+import numpy as np
+
+from scipy.optimize import linprog
+
+
+class DRSCGenerator:
+    """
+    Dirichlet-Rescale-Constraints (DRSC) Algorithm
+    Generates n-dimensional vectors with fixed sum satisfying linear and nonlinear constraints.
+    
+    This implementation is based on the paper "Generating Random Vectors 
+    satisfying Linear and Nonlinear Constraints" by Rick S. H. Willemsen, 
+    Wilco van den Heuvel and Michel van de Velden. It can be found here:
+    https://arxiv.org/abs/2501.16936
+    """
+    
+    def __init__(
+        self, 
+        n: int,
+        bounds: np.ndarray,
+        linear_constraints: Optional[List[Tuple[np.ndarray, float]]] = None,
+        nonlinear_constraints: Optional[List[Callable]] = None,
+        epsilon_ineq: float = 1e-3,
+        epsilon_eq: float = 1e-2,
+        max_restarts: int = 10000
+    ):
+        """
+        Initialize the DRSC generator.
+        
+        Parameters
+        ----------
+        * `n` [int]: 
+            Dimension of vectors to be generated
+            
+        * `bounds` [np.ndarray]: 
+            Array of shape (n, 2) with [lower, upper] bounds for each dimension
+            (normalized to sum = 1 simplex)
+        
+        * `linear_constraints`: 
+            List of (a, b) tuples representing a^T x <= b
+            
+        * `nonlinear_constraints`: 
+            List of functions f(x) that should be <= 0
+            
+        * `epsilon_ineq`: 
+            Feasibility tolerance for inequality constraints
+            
+        * `epsilon_eq`:
+            Feasibility tolerance for equality constraints
+        
+        * `max_restarts`:
+            Maximum number of restarts before giving up
+        """
+        self.n = n
+        self.original_bounds = np.array(bounds)
+        self.linear_constraints = linear_constraints or []
+        self.nonlinear_constraints = nonlinear_constraints or []
+        self.epsilon_ineq = epsilon_ineq
+        self.epsilon_eq = epsilon_eq
+        self.max_restarts = max_restarts
+        
+        # Compute optimal dimension ordering based on constrained dimension width
+        self._sort_order, self._inverse_order = self._compute_optimal_order()
+        
+        # Reorder bounds and constraints according to sort order
+        self._sorted_bounds = self.original_bounds[self._sort_order]
+        self._sorted_linear_constraints = self._reorder_linear_constraints()
+        self._sorted_nonlinear_constraints = self._reorder_nonlinear_constraints()
+        
+        # Find induced simplices using sorted dimensions (Algorithm 3, line 1)
+        self._thetas = self._find_induced_simplices()
+
+    
+    def _compute_optimal_order(self) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Compute optimal dimension ordering based on width.
+        
+        Dimensions with smaller widths (less flexible) are processed first.
+        This maximizes sampling efficiency by leaving flexible dimensions
+        to absorb variation from random sampling.
+        
+        Returns
+        -------
+            `sort_order` [np.ndarray]: 
+                Indices to sort dimensions by width (ascending)
+            `inverse_order`[np.ndarray]: 
+                Indices to restore original dimension order
+        """
+        widths = self.original_bounds[:, 1] - self.original_bounds[:, 0]
+        sort_order = np.argsort(widths)
+        inverse_order = np.argsort(sort_order)
+        
+        return sort_order, inverse_order
+    
+    
+    def _reorder_linear_constraints(self) -> List[Tuple[np.ndarray, float]]:
+        """Reorder linear constraint coefficients according to sort order."""
+        reordered = []
+        for a, b in self.linear_constraints:
+            a_reordered = a[self._sort_order]
+            reordered.append((a_reordered, b))
+        return reordered
+    
+    
+    def _reorder_nonlinear_constraints(self) -> List[Callable]:
+        """Wrap nonlinear constraints to work with sorted dimensions."""
+        wrapped = []
+        for f in self.nonlinear_constraints:
+            def wrapped_f(x_sorted, func=f):
+                # Convert from sorted order to original order before evaluating
+                x_original = x_sorted[self._inverse_order]
+                return func(x_original)
+            wrapped.append(wrapped_f)
+        return wrapped
+    
+    
+    def _find_induced_simplices(self) -> np.ndarray:
+        """
+        Find induced simplices Si for each dimension using LP formulation (1)-(7).
+        Uses sorted dimension order for optimal efficiency.
+        
+        For each dimension i, find theta_i such that the induced simplex
+        S_i = {x in S : x_i >= theta_i} covers part of the infeasible region
+        while not overlapping with previous induced simplices.
+        
+        The key insight from the paper is that induced simplices should cover
+        the INFEASIBLE region, so that when a point lands there, we can 
+        transform it back to the standard simplex.
+        
+        Returns
+        -------
+            `thetas` [np.ndarray]:
+                Array of theta values for each dimension
+        """
+        thetas = np.zeros(self.n)
+        
+        for i in range(self.n):
+            theta_i = self._solve_lp_for_dimension(i, thetas)
+            thetas[i] = theta_i
+        
+        return thetas
+    
+    
+    def _solve_lp_for_dimension(
+            self,
+            dim: int,
+            previous_thetas: np.ndarray
+        ) -> float:
+        """
+        Solve LP formulation (1)-(7) for dimension dim (in sorted order).
+        
+        The goal is to find the largest theta_i such that the induced simplex
+        S_i = {x : x_i >= theta_i, sum(x) = 1, x >= 0} lies within the 
+        infeasible region AND doesn't overlap with previous induced simplices.
+        
+        According to the paper, we maximize x_i subject to:
+        - sum(x) = 1 (on the simplex)
+        - 0 <= x <= 1 (standard simplex bounds)
+        - All linear constraints in J_lin are satisfied
+        - x_k < theta_k for all k < i (non-overlapping with previous simplices)
+        
+        Parameters
+        ----------
+            `dim` [int]:
+                Current dimension index (in sorted order)
+            `previous_thetas` [np.ndarray]:
+                Array of theta values from previous dimensions
+            
+        Returns
+        -------
+            `theta_i` [float]:
+                theta_i value for the induced simplex
+        """
+        # Objective: maximize x_i (we'll negate for minimization)
+        c = np.zeros(self.n)
+        c[dim] = -1  # Negative because linprog minimizes
+        
+        # Equality constraint: sum(x) = 1
+        A_eq = np.ones((1, self.n))
+        b_eq = np.array([1.0])
+        
+        # Inequality constraints
+        A_ub = []
+        b_ub = []
+        
+        # Box constraints: l <= x <= u (from the sorted bounds)
+        for i in range(self.n):
+            # Lower bound: x[i] >= lower  →  -x[i] <= -lower
+            a_lower = np.zeros(self.n)
+            a_lower[i] = -1.0
+            A_ub.append(a_lower)
+            b_ub.append(-self._sorted_bounds[i, 0] + self.epsilon_ineq)
+            
+            # Upper bound: x[i] <= upper
+            a_upper = np.zeros(self.n)
+            a_upper[i] = 1.0
+            A_ub.append(a_upper)
+            b_ub.append(self._sorted_bounds[i, 1] + self.epsilon_ineq)
+        
+        # Additional linear constraints (already reordered)
+        for a, b in self._sorted_linear_constraints:
+            A_ub.append(a)
+            b_ub.append(b + self.epsilon_ineq)
+        
+        # Non-overlapping constraints: x_k < theta_k for k < dim
+        # This ensures the new induced simplex doesn't overlap with previous ones
+        # We enforce x[k] <= theta_k - epsilon (strictly less than)
+        for k in range(dim):
+            if previous_thetas[k] > 0:
+                constraint = np.zeros(self.n)
+                constraint[k] = 1.0
+                A_ub.append(constraint)
+                # x[k] <= theta_k - small_epsilon to ensure non-overlap
+                b_ub.append(previous_thetas[k] - 1e-9)
+        
+        # Convert to arrays
+        A_ub = np.array(A_ub) if A_ub else None
+        b_ub = np.array(b_ub) if b_ub else None
+        
+        # Bounds: 0 <= x_i <= 1
+        bounds = [(0, 1) for _ in range(self.n)]
+        
+        # Solve LP
+        result = linprog(
+            c, 
+            A_ub=A_ub, 
+            b_ub=b_ub, 
+            A_eq=A_eq, 
+            b_eq=b_eq, 
+            bounds=bounds, 
+            method='highs',
+        )
+        
+        if result.success:
+            return result.x[dim]
+        else:
+            return 0.0
+    
+    
+    def _sample_flat_dirichlet(self) -> np.ndarray:
+        """
+        Sample from flat Dirichlet distribution (uniform on standard simplex).
+        """
+        y = np.random.exponential(1.0, self.n)
+        return y / np.sum(y)
+    
+    
+    def _check_all_constraints(self, x: np.ndarray) -> bool:
+        """
+        Check if x satisfies ALL constraints.
+        
+        For the DRSC algorithm to work correctly, we must check:
+        1. Box bounds (essential - feasible region is subset of standard simplex)
+        2. Additional linear constraints (if any)
+        3. Nonlinear constraints (if any)
+        
+        Parameters
+        ----------
+        * `x` [np.ndarray]:
+            Point to check (in sorted dimension order)
+            
+        Returns
+        -------
+        * bool
+            True if all constraints are satisfied
+        """
+        # Check box bounds - ESSENTIAL
+        # The standard simplex is [0,1]^n with sum=1, but our feasible region
+        # is defined by tighter bounds
+        for i in range(self.n):
+            if x[i] < self._sorted_bounds[i, 0] - self.epsilon_ineq:
+                return False
+            if x[i] > self._sorted_bounds[i, 1] + self.epsilon_ineq:
+                return False
+        
+        # Check additional linear constraints (if any)
+        for a, b in self._sorted_linear_constraints:
+            if np.dot(a, x) > b + self.epsilon_ineq:
+                return False
+        
+        # Check nonlinear constraints (if any)
+        for f in self._sorted_nonlinear_constraints:
+            if f(x) > self.epsilon_ineq:
+                return False
+        
+        return True
+    
+    
+    def _is_in_induced_simplex(self, x: np.ndarray) -> Optional[int]:
+        """
+        Check if x is in any induced simplex.
+        
+        Returns
+        -------
+            * Index of the induced simplex, or None if not in any
+        """
+        for idx, (theta, vertices, dim) in enumerate(self.induced_simplices):
+            if x[dim] >= theta:
+                return idx
+        return None
+    
+    
+    def _find_containing_induced_simplex(self, x: np.ndarray) -> Optional[int]:
+        """
+        Find which induced simplex contains point x.
+        
+        An induced simplex S_i is defined as {x in S : x_i >= theta_i}.
+        Due to the non-overlapping construction, a point can be in at most
+        one induced simplex.
+        
+        Returns
+        -------
+            * Index of the induced simplex containing x, or None if not in any
+        """
+        for i in range(self.n):
+            if self._thetas[i] > 0 and x[i] >= self._thetas[i]:
+                return i
+        return None
+    
+    
+    def _affine_transform_to_standard_simplex(self, x: np.ndarray, dim: int) -> np.ndarray:
+        """
+        Apply affine transformation from induced simplex S_dim to standard simplex.
+        
+        Parameters
+        ----------
+            * `x` [np.ndarray]: 
+                Point in induced simplex (sorted order)
+            * `dim` [int]: 
+                Index of the induced simplex dimension
+            
+        Returns
+        -------
+            * `x_new` [np.ndarray]:
+                Transformed point in standard simplex (sorted order)
+        """
+        theta = self._thetas[dim]
+        
+        # Create the transformed point
+        x_new = x.copy()
+        x_new[dim] = x[dim] - theta
+    
+        # Rescale to sum to 1
+        # The sum is now 1 - theta, so we divide by (1 - theta)
+        scale = 1.0 - theta
+        if scale > 1e-10:
+            x_new = x_new / scale
+        
+        return x_new
+    
+    
+    def generate(self) -> np.ndarray:
+        """
+        Generate a single vector using the DRSC algorithm (Algorithm 3).
+        
+        Returns
+        -------
+            * `x` [np.ndarray]:
+                n-dimensional vector with sum=1 satisfying all constraints. 
+                Vector is returned in ORIGINAL dimension order.
+        """
+        restarts = 0
+        max_transforms_per_sample = 1000
+        
+        while restarts < self.max_restarts:
+            # Line 2: Sample from flat Dirichlet (in sorted dimension space)
+            x = self._sample_flat_dirichlet()
+            transforms = 0
+
+            while transforms < max_transforms_per_sample:
+                # Line 3: Check if all constraints satisfied
+                if self._check_all_constraints(x):
+                    return x[self._inverse_order]
+                
+                # Line 5: Check if x is inside an induced simplex
+                simplex_idx = self._find_containing_induced_simplex(x)
+                
+                if simplex_idx is not None:
+                    # Line 6: Apply affine transformation
+                    x = self._affine_transform_to_standard_simplex(x, simplex_idx)
+                    transforms += 1
+                    
+                    # Safety check: if values are getting too large, restart
+                    if np.any(np.abs(x) > 1e6):
+                        restarts += 1
+                        break  # Will trigger restart
+                else:
+                    # Line 8: Restart - violated constraint but not in any induced simplex
+                    restarts += 1
+                    break
+        
+        raise RuntimeError(
+            f"Failed to generate valid vector after {self.max_restarts} restarts. "
+            f"The constraint region may be very small or empty."
+        )
+    
+    def generate_sample(self, size: int) -> np.ndarray:
+        """
+        Generate multiple vectors.
+        
+        Parameters
+        ----------
+            * `size`[int]: 
+                Number of vectors to generate
+            
+        Returns
+        -------
+            * Array of shape (size, n) with vectors in ORIGINAL dimension order
+        """
+        return np.array([self.generate() for _ in range(size)])
\ No newline at end of file
diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index 760bea09..c6e33bd6 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -1,12 +1,16 @@
+from typing import Union, List, Optional, Callable, Tuple
+
+import numpy as np
+
 from sklearn.utils import check_random_state
+
 from .space import Real, Integer, Categorical, Space, Task
-import numpy as np
-from scipy import linalg
-from typing import Union, List
+from .DRSC import DRSCGenerator
+
 
 class Constraints:
     def __init__(self, constraints_list, space):
-        """Constraints used when sampling for the aqcuisiiton function
+        """Constraints used when sampling for the acquisition function
 
         Parameters
         ----------
@@ -118,7 +122,7 @@ def rvs(
 
             n_samples_candidates += n_samples
             if n_samples_candidates > 100000 and len(rows) < 100:
-                # If we have less than a 1/10.000 succes rate on the sampling
+                # If we have less than a 1/10.000 success rate on the sampling
                 # we throw an error
                 raise RuntimeError(
                     '''Could not find valid samples in constrained space.
@@ -134,7 +138,8 @@ def sumequal_sampling(
             n_samples: int = 1,
             random_state: Union[int, np.random.RandomState, None] = None,
         ) -> List:
-        """Draw samples that respect SumEquals constraints.
+        """Draw samples that respect SumEquals constraints using the DRSC
+        algorithm.
 
         The samples are in the original space. They need to be transformed
         before being passed to a model or minimizer by `space.transform()`.
@@ -153,156 +158,102 @@ def sumequal_sampling(
         * `points`: [list of lists, shape=(n_points, n_dims)]
            Points sampled from the space.
         """
-        
-        def null_space(A, rcond=None) -> np.ndarray:
-            """Helper function to calculate the null space of a matrix
-
-            Parameters
-            ----------
-            A : numpy array
-                The matrix to calculate the null space of
-            rcond : float, optional
-                The tolerance for determining the effective rank of A.
-
-            Returns
-            -------
-            Q : numpy array
-                The null space of A, as a matrix with orthonormal columns.
-            """
-            
-            u, s, vh = np.linalg.svd(A, full_matrices=True) # A = u*s*vh
-            M, N = u.shape[0], vh.shape[1] # M = number of rows, N = number of columns
-            if rcond is None: # default value of rcond
-                rcond = np.finfo(s.dtype).eps * max(M, N) # machine precision times max dimension
-            tol = np.amax(s) * rcond # tolerance for singular values
-            num = np.sum(s > tol, dtype=int) # number of singular values greater than tol
-            Q = vh[num:,:].T.conj() # columns of vh corresponding to singular values greater than tol
-            return Q # return the null space of A
+        return self._drsc_sampling(n_samples, random_state)
 
-        rng = check_random_state(random_state)
-                
-        # Find a point on the plane defined by the SumEquals constraint where
-        # A + B + ... = value. We do this by asking where the diagonal between
-        # the origin and A_max, B_max, ... intersects the plane.
-        d = len(self.sum_equals[0].dimensions)
-        origin = np.array(
-            [self.space.bounds[dim][0] 
-             for dim in self.sum_equals[0].dimensions]
-        )
-        delta = np.array(
-            [self.space.bounds[dim][1] - self.space.bounds[dim][0] 
-             for dim in self.sum_equals[0].dimensions]
-        )
         
-        A = np.zeros((d,d))
-        B = np.zeros(d)
-        # Row representing the sum constraint
-        A[0,:] = 1
-        B[0] = self.sum_equals[0].value
-        # Rows that define the linear equation for the diagonal along the 
-        # constrained dimensions
-        for i in range(1,d):
-            A[i, 0] = -delta[i]/delta[0]
-            A[i, i] = 1
-            B[i] = origin[i]
-        # Identify the point that lies on the constraint plane and on the diagonal
-        point = np.linalg.solve(A, B)
-        # Use the fact that the vector [1, 1, ...] (a 1 for each constrained 
-        # dimension) is normal to the plane defined by A + B + ... to build
-        # basis-vectors inside the plane, using the null_space function
-        N = np.array(np.ones(d))
-        ns = null_space(N[np.newaxis, :])
-        # We only need to simulate points up to a distance of half the diagonal
-        # from the origin to A_max, B_max, etc.
-        sim_distance = np.sqrt(np.sum(delta**2)) / 2
+    def _drsc_sampling(self, n_samples, random_state):
+        """Sample using DRSC algorithm."""
+        rng = check_random_state(random_state)
+        constraint = self.sum_equals[0]
+        constrained_dims = constraint.dimensions
+        d = len(constrained_dims)
+        # Get bounds for constrained dimensions
+        bounds = np.array([self.space.bounds[dim] for dim in constrained_dims])
         
-        # To avoid "clustering" of points in the constrained plane, we will 
-        # create samples using low discrepancy quasirandom sequences, see:
-        # http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
-        # for background on this method
+        # Get or create DRSC generator (handles sorting internally)
+        drsc_gen = constraint.get_drsc_generator(self.space)
         
-        # Helper function for calculating the generalized golden ratio
-        def phi(d):
-            x = 2.0
-            for i in range(10): 
-                x = pow(1+x,1/(d+1)) 
-            return x
-        # Golden ratio for our present dimensionality (the constrained space)
-        g = phi(d-1)
-        alpha = np.zeros(d-1)
-        for j in range(d-1):
-            alpha[j] = pow(1/g, j+1) %1
-        vec_comp = np.zeros((1, d-1))
-        # Choose seed (starting location) in the normalized space
-        seed = 0.5
+        # Set numpy random seed for DRSC
+        np.random.seed(rng.randint(0, 2**31))
         
-        # Build a list of samples
         samples = []
-        j = 0
-        while len(samples) < n_samples:
-            # Generate next step in the sequence
-            vec_comp = (seed + alpha*(j+1)) %1
-            j += 1
-            # Center these components on zero
-            vec_comp = vec_comp - 0.5
-            # Simulate lengths of each null_space vector to add to our point
-            vec_comp = vec_comp * sim_distance * 2
-            sample_candidate = point[:, None].T + ns @ vec_comp.T
-            # Generate the correct shape 
-            sample_candidate = sample_candidate[0]
-            # Check that the candidate is inside the original parameter space
-            inspace = [
-                (sample_candidate[i] >= self.space.bounds[dim][0]) and
-                (sample_candidate[i] <= self.space.bounds[dim][1])
-                for i, dim in enumerate(self.sum_equals[0].dimensions)
-            ]
-            # Only accept the candidate if it is in our space
-            if all(inspace):
-                samples.append(sample_candidate)      
-        # Convert the list of arrays to a list of lists
-        samples = [arr.tolist() for arr in samples]
         
-        # Create settings for the dimensions that are not part of the constraint
-        if d < self.space.n_dims:
-            remaining_dimensions = [
-                i for i in range(self.space.n_dims) 
-                if i not in self.sum_equals[0].dimensions
-            ]
-            # Convert our list of samples to an array
-            samples = np.array(samples)
-            # Expand the sample array to make space for settings of the 
-            # unconstrained dimensions
-            for i in remaining_dimensions:
-                samples = np.insert(
-                    samples,
-                    i, 
-                    np.zeros(len(samples)), 
-                    axis=1
-                )
-            # Convert back to list of lists
-            samples = samples.tolist()
-            
-            # Generate random settings across all factors
-            full_sample = self.space.rvs(n_samples=n_samples, random_state=rng)
-            # Sort the settings in each column, which will ensure that the
-            # unconstrained settings are distributed in a space-filling way too
-            transposed_sample = list(zip(*full_sample))
-            sort_trans_sample = [sorted(inner_tuple) for inner_tuple in transposed_sample]
-            # Place the settings back in a list of lists
-            full_sample = list(map(list, zip(*sort_trans_sample)))
+        for _ in range(n_samples):
+            # Generate sample on simplex, then convert to original space
+            x_simplex = drsc_gen.generate()
+            x_constrained = x_simplex * constraint.value              
             
-            # Overwrite the random setting values for the constrained factors
-            for j in remaining_dimensions:
-                for i in range(len(samples)):
-                    samples[i][j] = full_sample[i][j]
+            # If the dimensions in the space have very different widths we risk 
+            # violating the space bounds due to the numerical tolerance of the
+            # DRSC algorithm. Fix this, while preserving the sum
+            x_constrained = self._fix_bounds_preserving_sum(
+                x_constrained, 
+                bounds, 
+                constraint.value
+            )
+                        
+            # Build full sample
+            if d < self.space.n_dims:
+                # Generate random settings across all factors
+                full_sample = self.space.rvs(n_samples=1, random_state=rng)[0]
+                # Overwrite the random setting values for the constrained factors
+                for i, dim_idx in enumerate(constrained_dims):
+                    full_sample[dim_idx] = x_constrained[i]   
+                samples.append(full_sample)
+            else:
+                samples.append(x_constrained.tolist())
+        
+        return samples
+    
+   
+    def _fix_bounds_preserving_sum(self, x, bounds, target_sum, max_iterations=100):
+        """
+        Adjust x to satisfy bounds while preserving the sum constraint.
+        
+        When a value exceeds its bounds, the excess is redistributed to other
+        dimensions that have room to absorb it.
+        """
+        x = x.copy()
+        # Rescale proportionally to restore sum of x to target value
+        x = x * (target_sum / np.sum(x))
+        
+        for _ in range(max_iterations):
+            all_satisfied = True
             
-            # Shuffle the order of the samples, otherwise the unconstrained
-            # settings will be returned in a sorted order. Use seeding to 
-            # provide consistent initial samples
-            rng2 = np.random.default_rng(seed=42)
-            rng2.shuffle(samples)
+            for i in range(len(x)):
+                if x[i] < bounds[i, 0]:
+                    deficit = bounds[i, 0] - x[i]
+                    x[i] = bounds[i, 0]
                     
-        return samples
+                    # Subtract from other dimensions that have room
+                    for j in range(len(x)):
+                        if j != i and deficit > 0:
+                            room = x[j] - bounds[j, 0]
+                            transfer = min(room, deficit)
+                            x[j] -= transfer
+                            deficit -= transfer
+                    
+                    all_satisfied = False
+                    
+                elif x[i] > bounds[i, 1]:
+                    excess = x[i] - bounds[i, 1]
+                    x[i] = bounds[i, 1]
+                    
+                    # Add to other dimensions that have room
+                    for j in range(len(x)):
+                        if j != i and excess > 0:
+                            room = bounds[j, 1] - x[j]
+                            transfer = min(room, excess)
+                            x[j] += transfer
+                            excess -= transfer
+                    
+                    all_satisfied = False
+            
+            if all_satisfied:
+                break
+        
+        return x
+
 
     def validate_sample(self, sample: List) -> bool:
         """ Validates a sample of parameter values in regards to the
@@ -320,7 +271,7 @@ def validate_sample(self, sample: List) -> bool:
         """
 
         # We iterate through all the dimensions and check the the type of
-        # constriants that are applied to a single dimensions, i.e Single,
+        # constraints that are applied to a single dimensions, i.e Single,
         # Exclusive and Inclusive
         #
         # We iterate through all samples which corresponds to number of
@@ -334,7 +285,7 @@ def validate_sample(self, sample: List) -> bool:
             # Inclusive constraints.
             # Check if there is a least one inclusive constraint:
             if self.inclusive[dim]:
-                # We go through all inlcusive constraints for this dimension
+                # We go through all inclusive constraints for this dimension
                 # and if the value is not found to be included in any of the
                 # bounds of the inclusive constraints we return false.
                 value_is_valid = False
@@ -352,11 +303,11 @@ def validate_sample(self, sample: List) -> bool:
                 if not constraint.validate_constraint(sample[dim]):
                     return False
 
-        # We iterate through sum constriants
+        # We iterate through sum constraints
         for constraint in self.sum:
             if not constraint.validate_sample(sample):
                 return False
-        # We iterate through sum_equals constriants
+        # We iterate through sum_equals constraints
         for constraint in self.sum_equals:
             if not constraint.validate_sample(sample):
                 return False
@@ -364,8 +315,7 @@ def validate_sample(self, sample: List) -> bool:
         for constraint in self.conditional:
             if not constraint.validate_sample(sample):
                 return False
-        # If we we did not find any violaiton of the constraints we return
-        # True.
+        # If we did not find any violation of the constraints we return True
         return True
 
     def __repr__(self):
@@ -496,7 +446,7 @@ def __init__(self, dimension, bounds, dimension_type):
             consist of floats.
 
             For 'categorical' dimensions the tuple must be of length < number
-            of dimensions and lenght > 1.
+            of dimensions and length > 1.
 
             The tuple can contain any combination of str, int or float
 
@@ -563,7 +513,7 @@ def __init__(self, dimension, bounds, dimension_type):
             consist of floats.
 
             For 'categorical' dimensions the tuple must be of length < number
-            of dimensions and lenght > 1.
+            of dimensions and length > 1.
 
             The tuple can contain any combination of str, int or float
 
@@ -616,7 +566,7 @@ def __init__(self, dimensions, value, less_than=True):
         Parameters
         ----------
         * `dimensions` [list of ints]:
-            A list of integers coresponding to the index of the dimensions that
+            A list of integers corresponding to the index of the dimensions that
             should be summed
 
         * `value` [float or int]:
@@ -663,7 +613,13 @@ def __eq__(self, other):
             return False
 
 class SumEquals():
-    def __init__(self, dimensions: List[int], value: Union[float, int]):
+    def __init__(
+        self, 
+        dimensions: List[int], 
+        value: Union[float, int],
+        linear_constraints: Optional[List[Tuple[np.ndarray, float]]] = None,
+        nonlinear_constraints: Optional[List[Callable]] = None,
+    ):
         """Constraint class of type SumEquals.
 
         This constraint enforces that the sum of all values drawn for the
@@ -678,6 +634,12 @@ def __init__(self, dimensions: List[int], value: Union[float, int]):
 
         * `value` [float or int]:
             The value for which the sum should be equal to.
+            
+        * `linear_constraints` [list of (a, b) tuples, optional]:
+            Additional linear constraints a^T x <= b (only used with DRSC)
+            
+        * `nonlinear_constraints` [list of callables, optional]:
+            Additional nonlinear constraints f(x) <= 0 (only used with DRSC)
 
         """
 
@@ -694,9 +656,104 @@ def __init__(self, dimensions: List[int], value: Union[float, int]):
 
         self.dimensions = tuple(dimensions)
         self.value = value
+        self.linear_constraints = linear_constraints or []
+        self.nonlinear_constraints = nonlinear_constraints or []
+        # Lazy initialization of DRSC generator
+        self._drsc_generator = None
 
         self.validate_sample = self._validate_sample
-
+        
+    def get_drsc_generator(self, space) -> DRSCGenerator:
+        """
+        Get or create DRSC generator with space information.
+        
+        Parameters
+        ----------
+            space: Space object containing dimension bounds
+            
+        Returns
+        -------
+            DRSCGenerator instance configured for this constraint
+        """
+        if self._drsc_generator is None:
+            # Extract bounds for constrained dimensions
+            bounds = [space.bounds[dim] for dim in self.dimensions]
+            
+            # Normalize bounds to simplex (sum=1)
+            normalized_bounds = [(low / self.value, high / self.value)  for low, high in bounds]
+            
+            # Create DRSC generator
+            self._drsc_generator = DRSCGenerator(
+                n=len(self.dimensions),
+                bounds=normalized_bounds,
+                linear_constraints=self._convert_linear_constraints(bounds),
+                nonlinear_constraints=self._wrap_nonlinear_constraints(bounds),
+            )
+        return self._drsc_generator
+    
+    def _convert_linear_constraints(self, bounds):
+        """
+        Convert user-provided linear constraints to normalized form for use 
+        in the DRSC format.
+        """
+        # DRSC works on unit simplex, need to normalize constraints
+        converted = []
+        
+        # Add box constraints from dimension bounds
+        for i, (low, high) in enumerate(bounds):
+            # x_i >= low becomes constraint on simplex
+            # After scaling: x_simplex[i] * value >= low
+            # In simplex coords: -x[i] <= -low/value
+            a = np.zeros(len(self.dimensions))
+            a[i] = -1.0
+            b = -low / self.value
+            converted.append((a, b))
+            
+            # x_i <= high becomes: x[i] <= high/value
+            a = np.zeros(len(self.dimensions))
+            a[i] = 1.0
+            b = high / self.value
+            converted.append((a, b))
+        
+        # Add user-provided linear constraints (if any)
+        for a, b in self.linear_constraints:
+            # Extract only the constrained dimensions
+            a_constrained = np.array([a[dim] for dim in self.dimensions])
+            # Normalize by the sum value
+            a_normalized = a_constrained / self.value
+            b_normalized = b / self.value
+            converted.append((a_normalized, b_normalized))
+        
+        return converted
+    
+    def _wrap_nonlinear_constraints(self, bounds):
+        """Wrap nonlinear constraints to work on simplex."""
+        if not self.nonlinear_constraints:
+            return []
+        
+        wrapped = []
+        for f in self.nonlinear_constraints:
+            def wrapped_f(x_simplex):
+                # Convert from simplex to original space
+                x_original = self._simplex_to_original(x_simplex, bounds)
+                return f(x_original)
+            wrapped.append(wrapped_f)
+        return wrapped
+    
+    def _simplex_to_original(self, x_simplex, bounds):
+        """Convert from unit simplex coordinates to original space."""
+        x_original = np.zeros(len(self.dimensions))
+        for i, (low, high) in enumerate(bounds):
+            # Scale from simplex (sums to 1) to original space (sums to value)
+            scaled_value = x_simplex[i] * self.value
+            # Clip to dimension bounds to handle numerical tolerance and avoid errors
+            x_original[i] = np.clip(scaled_value, low, high)
+        return x_original
+    
+    def _original_to_simplex(self, x_original):
+        """Convert from original space to unit simplex."""
+        return np.array(x_original) / self.value
+    
     def _validate_sample(self, sample: List[int]) -> bool:
         # Returns True if sample does not violate the constraints aside from 
         # floating point errors
@@ -839,7 +896,7 @@ def check_constraints(space, constraints):
                 raise IndexError('Dimension index exceeds number of dimensions')
             for ind_dim in constraint.dimensions:
                 if isinstance(space.dimensions[ind_dim], Categorical):
-                    raise ValueError('Sum constraint can not be applid to categorical dimension: {}'.format(space.dimensions[ind_dim]))
+                    raise ValueError('Sum constraint can not be applied to categorical dimension: {}'.format(space.dimensions[ind_dim]))
         elif isinstance(constraint, SumEquals):
             # Check that there is only one SumEquals constraint being applied
             if len(constraints) > 1:
@@ -874,7 +931,9 @@ def check_constraints(space, constraints):
                 for constraint_if_false in constraint.if_false:
                     check_constraints(space, [constraint_if_false])
         else:
-            raise TypeError('Constraints must be of type "Single", "Exlusive", "Inclusive", "Sum", "SumEquals" or "Conditional". Got {}'.format(type(constraint)))
+            raise TypeError(
+                'Constraints must be of type "Single", "Exclusive", "Inclusive", "Sum", "SumEquals" or "Conditional". Got {}'.format(type(constraint))
+            )
 
 
 def check_dim_and_space(space, constraint):
@@ -941,7 +1000,7 @@ def check_value(dim, value):
             raise ValueError('Value {} exceeds bounds of space {}'.format(value, [dim.low, dim.high]))
     else:  # Categorical dimension.
         if value not in dim.categories:
-            raise ValueError('Categorical value {} is not in space with categoreis {}'.format(value, dim.categories))
+            raise ValueError('Categorical value {} is not in space with categories {}'.format(value, dim.categories))
 
 
 def check_is_constraint(constraint):
@@ -951,4 +1010,4 @@ def check_is_constraint(constraint):
     if not (isinstance(constraint, Single) or isinstance(constraint, Inclusive)
             or isinstance(constraint, Exclusive) or isinstance(constraint, Sum)
             or isinstance(constraint, Conditional)):
-        raise TypeError('Constraint must be of type Inclusive, Exlusive, Single, Sum or Conditional. Got {}'.format(type(constraint)))
+        raise TypeError('Constraint must be of type Inclusive, Exclusive, Single, Sum or Conditional. Got {}'.format(type(constraint)))
diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index 5124d96f..10e62dde 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -216,7 +216,7 @@ def test_SumEquals():
     assert not cons.validate_sample([3.00006, 2.0, "A"])    
     
     # Check again that only valid samples are drawn
-    samples = cons.sumequal_sampling(n_samples=1000)
+    samples = cons.sumequal_sampling(n_samples=10000)
     for sample in samples:
         factor_sum = np.sum(sample[0] + sample[1])
         assert np.isclose(factor_sum, cons.sum_equals[0].value)
@@ -225,10 +225,11 @@ def test_SumEquals():
     space = Space([[10.0, 20.0], [1000.0, 2000.0], [1, 5], ["A", "B"]])
     cons = Constraints([SumEquals((0, 1), 1234)], space)
     # Check again that only valid samples are drawn
-    samples = cons.sumequal_sampling(n_samples=1000)
+    samples = cons.sumequal_sampling(n_samples=10000)
     for sample in samples:
         factor_sum = np.sum(sample[0] + sample[1])
-        assert np.isclose(factor_sum, cons.sum_equals[0].value)
+        # The default tightness of the SumEquals sampling should be within 1e-4
+        assert np.isclose(factor_sum, cons.sum_equals[0].value, rtol=1e-4)
     # Check that the samples also have settings for the other dimensions
     assert len(samples[0]) == 4 and len(samples[999]) == 4
     # Check that the other dimensions have correct type
@@ -276,6 +277,49 @@ def test_SumEquals():
     # The values suggested for the constrained dimensions should be the same 
     # irrespective of their order in the dimension list
     assert x1[:3] == x2[1:]
+    
+    # Test that we can ask for all initial points as long as we use strategy=cl_min
+    assert opt.ask(5, strategy="cl_min")
+    # We can also ask for points beyond the initial ones without data
+    assert opt.ask(6, strategy="cl_min")
+    # Give the optimizer some data and test that we can ask for multiple points
+    # after we have data
+    y = 0
+    for _ in range(10):
+        x = opt.ask()
+        opt.tell(x, y)
+        y -= 0.1
+    assert opt.ask(3, strategy="cl_min")
+    
+    # Test that different seeds of the Optimizer provides different initial points
+    space = [
+        (50., 100.),
+        (1., 5.),
+        (3., 6.),
+        (1., 35.),
+        (1., 5.),
+        ("A", "B", "C"),
+        (1, 10),
+    ]
+    cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
+    # Build optimziers with different seeds
+    opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
+    opt1.set_constraints(cons)
+    opt2 = Optimizer(space, lhs=False, n_initial_points=10, random_state=2)
+    opt2.set_constraints(cons)
+    # Ask for the initial points
+    x1 = opt1.ask(10, strategy="cl_min")
+    x2 = opt2.ask(10, strategy="cl_min")
+    assert not x1 == x2
+    
+    # Check that seeding the generator directly inside each optimizer produces
+    # consistent points, despite the optimizers having different random states
+    constraint1 = opt1.get_constraints()
+    constraint2 = opt2.get_constraints()
+    # Ask for lots of points
+    x1 = constraint1.sumequal_sampling(n_samples=10000, random_state=31031988)
+    x2 = constraint2.sumequal_sampling(n_samples=10000, random_state=31031988)
+    assert x1 == x2
 
 @pytest.mark.fast_test
 def test_Conditional():
diff --git a/examples/design_qualifications/DRSC SumEqual implementation.py b/examples/design_qualifications/DRSC SumEqual implementation.py
new file mode 100644
index 00000000..a22fcba3
--- /dev/null
+++ b/examples/design_qualifications/DRSC SumEqual implementation.py	
@@ -0,0 +1,230 @@
+# -*- coding: utf-8 -*-
+"""
+A script demonstrating the features of the DRSC algorithm for use with the
+SumEquals constraint in ProcessOptimizer.
+
+Author: Morten Bormann Nielsen, Danish Technological Institute
+December 2025
+"""
+
+
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+from ProcessOptimizer import Optimizer
+from ProcessOptimizer.space import Real
+from ProcessOptimizer.space.constraints import SumEquals
+
+from ProcessOptimizer.model_systems import hart3
+from ProcessOptimizer.plots import plot_objective_1d
+
+#%% Setup of SumEquals sampling for very oblong space
+
+# The DRSC generator for SumEquals sampling is able to generate valid points
+# extremely fast
+space = [
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),        
+]
+# Define the constraint as all factor settings adding up to 100
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+# Extract the generator itself
+constraint = opt.get_constraints().sum_equals[0]
+drsc_gen = constraint.get_drsc_generator(opt.space)
+# Demonstrate speed of algorithm by generating 10,000 points. This call takes
+# about 12 seconds on an Intel i7-12850HX.
+start_time = time.time()
+x_simplex = drsc_gen.generate_sample(10000)
+print("Generated 10,000 points in the space in %s seconds" % (time.time() - start_time))
+
+#%% Demonstrate typical use for a real user with fewer points
+
+# Define the space
+space = [
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),        
+]
+# Set up the constraint
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
+opt = Optimizer(space, lhs=False, n_initial_points=10, random_state=31031988)
+opt.set_constraints(cons)
+# Ask for settings for the initial experiments
+print(time.strftime("Starting opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
+x = opt.ask(10, strategy="cl_min")
+print(time.strftime("Finished opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
+
+#%% Demonstrate that different seeds lead to different points
+
+opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
+opt1.set_constraints(cons)
+opt2 = Optimizer(space, lhs=False, n_initial_points=10, random_state=2)
+opt2.set_constraints(cons)
+
+x1 = opt1.ask(10, strategy="cl_min")
+x2 = opt2.ask(10, strategy="cl_min")
+
+x!=x2
+
+#%% Demonstrate that identical seeds lead to identical points
+
+constraint1 = opt1.get_constraints()
+constraint2 = opt2.get_constraints()
+
+x1 = constraint1.sumequal_sampling(n_samples=1000, random_state=31031988)
+x2 = constraint2.sumequal_sampling(n_samples=1000, random_state=31031988)
+
+x1 == x2
+
+#%% Demonstrate how to use SumEqual constraints when categoricals are present
+space = [
+    ("A", "B"),
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),
+    ("C", "D", "E"),
+    (1, 5),
+]
+cons = [SumEquals(dimensions=[1, 2, 3, 4, 5], value=100.0)]
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+print(time.strftime("Starting opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
+x = opt.ask(10, strategy="cl_min")
+print(time.strftime("Finished opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
+
+#%% Use the new sampling method with simulated data
+
+# Build Hartmann 3D ModelSystem object
+hart3_model = hart3.create_hart3(noise=True)
+# Define space
+space = [
+    Real(0., 1., name='x0'),
+    Real(0., 1., name='x1'),
+    Real(0., 1., name='x2'),
+]
+
+cons = [SumEquals(dimensions=[0, 1, 2], value=1.0)]
+opt = Optimizer(space, lhs=False, n_initial_points=20)
+opt.set_constraints(cons)
+
+# Run optimization, first through the initial points, then using EI  
+for _ in range(30):
+    x = opt.ask(5, strategy="cl_min")
+    y = [hart3_model.get_score(x) for x in x]
+    res = opt.tell(x, y)
+
+# Show the system
+plot_objective_1d(res, pars="expected_minimum")
+
+# Show the location of the sampled points in this experiment
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+cmap = mpl.colormaps["viridis"].resampled(len(opt.Xi))
+
+idx = np.arange(len(opt.Xi))
+x = np.array(opt.Xi)
+y = np.array(opt.yi)
+
+ax.scatter(x[:,0], x[:, 1], x[:, 2], c=idx, cmap=cmap, alpha=1)
+
+ax.set_xlim(0, 1)
+ax.set_ylim(0, 1)
+ax.set_zlim(0, 1)
+
+
+#%% Demonstrate the use of a non-linear SumEquals constraint 
+
+# Define space
+space = [
+    Real(0., 1., name='x0'),
+    Real(0., 1., name='x1'),
+    Real(0., 1., name='x2'),
+]
+
+sum_value = 1.0
+
+# Nonlinear constraint: x0 * x1 < 0.1
+# Convention: function returns value that should be <= 0 when SATISFIED
+def product_constraint(x):
+    """x[0] * x[1] <= 0.1"""
+    return x[0] * x[1] - 0.1
+
+# Create SumEquals constraint with nonlinear constraint
+cons = [
+    SumEquals(
+        dimensions=[0, 1, 2],
+        value=sum_value,
+        nonlinear_constraints=[product_constraint]
+    )
+]
+
+# Create optimizer
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+# Test sampling
+constraint = opt.get_constraints().sum_equals[0]
+drsc_gen = constraint.get_drsc_generator(opt.space)
+
+# Generate samples
+samples = drsc_gen.generate_sample(1000)
+
+# Convert to original space and verify constraints
+samples_original = samples * sum_value
+
+x = samples_original
+# Create a plot showing where the samples lie in the space
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+ax.scatter(x[:, 0], x[:, 1], x[:, 2],)
+ax.set_xlim(0, 1)
+ax.set_ylim(0, 1)
+ax.set_zlim(0, 1)
+
+
+#%% Demonstrate the use of three constrained dimensions and a categorical
+
+# Sample 1000 points and plot them to show it works
+space = [
+    (0., 10.),
+    (0., 10.),
+    (0., 10.),
+    ("A", "B", "C"),
+]
+
+cons = [SumEquals(dimensions=[0, 1, 2], value=15.0)]
+opt = Optimizer(space, lhs=False, n_initial_points=100)
+opt.set_constraints(cons)
+
+constraint = opt.get_constraints().sum_equals[0]
+
+x = opt.ask(100, strategy="cl_min")
+colors = []
+for xx in x:
+    if xx[3] == "A":
+        colors.append("r")
+    elif xx[3] == "B":
+        colors.append("g")
+    else:
+        colors.append("b")
+    
+# Create a plot showing the distribution of points in this space
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+[ax.scatter(x[0], x[1], x[2], c=colors[i]) for i, x in enumerate(x)]
+ax.set_xlim(0, 10)
+ax.set_ylim(0, 10)
+ax.set_zlim(0, 10)
diff --git a/pyproject.toml b/pyproject.toml
index 2e94aa42..0730b9d1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,7 +19,6 @@ classifiers = [
   "License :: OSI Approved :: BSD License",
   "Operating System :: OS Independent",
   "Programming Language :: Python :: 3",
-  "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",
   "Programming Language :: Python :: 3.12",