diff --git a/PEPit/constraint.py b/PEPit/constraint.py
index fd28e74e..7aa64ce3 100644
--- a/PEPit/constraint.py
+++ b/PEPit/constraint.py
@@ -8,6 +8,7 @@ class Constraint(object):
 
     Attributes:
         name (str): A name set through the set_name method. None is no name is given.
+        activated (bool): A boolean flag used to activate/deactivate the Constraint in the PEP.
         expression (Expression): The :class:`Expression` that is compared to 0.
         equality_or_inequality (str): "equality" or "inequality". Encodes the type of constraint.
         _value (float): numerical value of `self.expression` obtained after solving the PEP via SDP solver.
@@ -56,13 +57,16 @@ def __init__(self,
             AssertionError: if provided `equality_or_inequality` argument is neither "equality" nor "inequality".
 
         """
-        # Initialize name of the constraint
-        self.name = None
+        # Initialize the activated attribute to True
+        self.activated = True
 
         # Update the counter
         self.counter = Constraint.counter
         Constraint.counter += 1
 
+        # Initialize name of the constraint
+        self.name = "Constraint {}".format(self.counter)
+
         # Store the underlying expression
         self.expression = expression
 
@@ -92,6 +96,20 @@ def get_name(self):
         """
         return self.name
 
+    def activate(self):
+        """
+        Activate the use of the Constraint.
+
+        """
+        self.activated = True
+
+    def deactivate(self):
+        """
+        Deactivate the use of the Constraint.
+
+        """
+        self.activated = False
+
     def eval(self):
         """
         Compute, store and return the value of the underlying :class:`Expression` of this :class:`Constraint`.
diff --git a/PEPit/pep.py b/PEPit/pep.py
index 82653f09..8fd32905 100644
--- a/PEPit/pep.py
+++ b/PEPit/pep.py
@@ -37,8 +37,11 @@ class PEP(object):
         wrapper_name (str): name of the used wrapper.
         wrapper (Wrapper): :class:`Wrapper` object that interfaces between the :class:`PEP` and the solver.
 
-        _list_of_constraints_sent_to_wrapper (list): list of :class:`Constraint` objects sent to the wrapper.
-        _list_of_psd_sent_to_wrapper (list): list of :class:`PSDMatrix` objects sent to the wrapper.
+        _list_of_prepared_constraints (list): list of :class:`Constraint` objects ready to be sent to the wrapper.
+        _list_of_prepared_psd (list): list of :class:`PSDMatrix` objects ready to be sent to the wrapper.
+
+        _list_of_constraints_sent_to_wrapper (list): list of :class:`Constraint` objects actually sent to the wrapper.
+        _list_of_psd_sent_to_wrapper (list): list of :class:`PSDMatrix` objects actually sent to the wrapper.
 
         objective (Expression): the expression to be maximized by the solver.
                                 It is set by the method `solve`. And should not be updated otherwise.
@@ -87,6 +90,8 @@ def __init__(self):
 
         # Initialize lists of constraints that will be sent to the wrapper to solve the SDP.
         # Those lists should not be updated by hand, only the solve method does update them.
+        self._list_of_prepared_constraints = list()
+        self._list_of_prepared_psd = list()
         self._list_of_constraints_sent_to_wrapper = list()
         self._list_of_psd_sent_to_wrapper = list()
 
@@ -330,6 +335,13 @@ def solve(self, wrapper="cvxpy", return_primal_or_dual="dual", safe_mode=True, v
             float: Worst case guarantee of the PEP.
 
         """
+        # Prepare the list of constraints to be given to the wrapper
+        # This runs only once!
+        # Then the method solve directly interacts with the solvers without browsing all PEPit's data.
+        if self.objective is None:
+            self._prepare_constraints(verbose=verbose)
+
+        # Prepare wrapper and solver
         wrapper_name = wrapper.lower()
 
         # Check that the solver is installed, if it is not, switch to CVXPY.
@@ -356,53 +368,24 @@ def solve(self, wrapper="cvxpy", return_primal_or_dual="dual", safe_mode=True, v
         self.wrapper_name = wrapper_name
         self.wrapper = wrapper
 
-        # Call the internal solve methods, which formulates and solves the PEP via the SDP solver.
+        # Call the internal solving methods, which formulates and solves the PEP via the SDP solver.
         out = self._solve_with_wrapper(wrapper, verbose, return_primal_or_dual,
                                        dimension_reduction_heuristic,
                                        eig_regularization, tol_dimension_reduction, **kwargs)
         return out
 
-    def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
-                            dimension_reduction_heuristic=None, eig_regularization=1e-3, tol_dimension_reduction=1e-4,
-                            **kwargs):
+    def _prepare_constraints(self, verbose=1):
         """
-        Internal solve method. Translate the :class:`PEP` to an SDP, and solve it via the wrapper.
+        Prepare the lists of scalar and lmi constraints that can be used.
+        Those are stored in the attributes `_list_of_prepared_constraints` and  `_list_of_prepared_psd`.
 
         Args:
-            wrapper (Wrapper): Interface to the solver.
             verbose (int, optional): Level of information details to print
-                                     (Override the CVXPY solver verbose parameter).
-
-                            - 0: No verbose at all
-                            - 1: PEPit information is printed but not CVXPY's
-                            - 2: Both PEPit and solver details are printed
-            return_primal_or_dual (str, optional): If "dual", it returns a worst-case upper bound of the PEP
-                                                   (dual value of the objective).
-                                                   If "primal", it returns a worst-case lower bound of the PEP
-                                                   (primal value of the objective).
-                                                   Default is "dual".
-                                                   Note both value should be almost the same by strong duality.
-            dimension_reduction_heuristic (str, optional): An heuristic to reduce the dimension of the solution
-                                                           (rank of the Gram matrix). Set to None to deactivate
-                                                           it (default value). Available heuristics are:
-
-                                                            - "trace": minimize :math:`Tr(G)`
-                                                            - "logdet{an integer n}": minimize
-                                                              :math:`\\log\\left(\\mathrm{Det}(G)\\right)`
-                                                              using n iterations of local approximation problems.
-
-            eig_regularization (float, optional): The regularization we use to make
-                                                  :math:`G + \\mathrm{eig_regularization}I_d \succ 0`.
-                                                  (only used when "dimension_reduction_heuristic" is not None)
-                                                  The default value is 1e-5.
-            tol_dimension_reduction (float, optional): The error tolerance in the heuristic minimization problem.
-                                                       Precisely, the second problem minimizes "optimal_value - tol"
-                                                       (only used when "dimension_reduction_heuristic" is not None)
-                                                       The default value is 1e-5.
-            kwargs (keywords, optional): Additional solver-specific arguments.
+                                     (Override the solver verbose parameter).
 
-        Returns:
-            float: Worst-case guarantee of the PEP.
+                                     - 0: No verbose at all
+                                     - 1: PEPit information is printed but not solver's
+                                     - 2: Both PEPit and solver details are printed
 
         """
 
@@ -423,15 +406,9 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         # Create an expression that serve for the objective (min of the performance measures)
         self.objective = Expression(is_leaf=True)
 
-        # Report the creation of variables (G, F)
-        if verbose:
-            print('(PEPit) Setting up the problem:'
-                  ' size of the Gram matrix: {}x{}'.format(Point.counter, Point.counter))
-        wrapper.set_main_variables()
-
         # Initialize the lists of constraints sent to wrapper
-        self._list_of_constraints_sent_to_wrapper = list()
-        self._list_of_psd_sent_to_wrapper = list()
+        self._list_of_prepared_constraints = list()
+        self._list_of_prepared_psd = list()
 
         # Defining performance metrics
         # Note maximizing the minimum of all the performance metrics
@@ -440,8 +417,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         for performance_metric in self.list_of_performance_metrics:
             assert isinstance(performance_metric, Expression)
             performance_metric_constraint = (self.objective <= performance_metric)
-            wrapper.send_constraint_to_solver(performance_metric_constraint)
-            self._list_of_constraints_sent_to_wrapper.append(performance_metric_constraint)
+            self._list_of_prepared_constraints.append(performance_metric_constraint)
 
         if verbose:
             print('(PEPit) Setting up the problem:'
@@ -451,8 +427,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
         if verbose:
             print('(PEPit) Setting up the problem: Adding initial conditions and general constraints ...')
         for condition in self.list_of_constraints:
-            wrapper.send_constraint_to_solver(condition)
-            self._list_of_constraints_sent_to_wrapper.append(condition)
+            self._list_of_prepared_constraints.append(condition)
         if verbose:
             print('(PEPit) Setting up the problem:'
                   ' initial conditions and general constraints ({} constraint(s) added)'.format(
@@ -463,8 +438,11 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
             if verbose:
                 print('(PEPit) Setting up the problem: {} lmi constraint(s) added'.format(len(self.list_of_psd)))
             for psd_counter, psd_matrix in enumerate(self.list_of_psd):
-                wrapper.send_lmi_constraint_to_solver(psd_counter, psd_matrix)
-                self._list_of_psd_sent_to_wrapper.append(psd_matrix)
+                # Print a message if verbose mode activated
+                if verbose > 0:
+                    print('\t\t Size of PSD matrix {}: {}x{}'.format(psd_counter + 1, *psd_matrix.shape))
+                # Add the PSD matrix to the list of prepared PSDs
+                self._list_of_prepared_psd.append(psd_matrix)
 
         # Defining class constraints
         if verbose:
@@ -479,8 +457,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
                       'scalar constraint(s) ...')
 
             for constraint in function.list_of_class_constraints:
-                wrapper.send_constraint_to_solver(constraint)
-                self._list_of_constraints_sent_to_wrapper.append(constraint)
+                self._list_of_prepared_constraints.append(constraint)
 
             if verbose:
                 print('\t\t\tFunction', function_counter, ':', len(function.list_of_class_constraints),
@@ -492,8 +469,11 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
                           'lmi constraint(s) ...')
 
                 for psd_counter, psd_matrix in enumerate(function.list_of_class_psd):
-                    wrapper.send_lmi_constraint_to_solver(psd_counter, psd_matrix)
-                    self._list_of_psd_sent_to_wrapper.append(psd_matrix)
+                    # Print a message if verbose mode activated
+                    if verbose > 0:
+                        print('\t\t Size of PSD matrix {}: {}x{}'.format(psd_counter + 1, *psd_matrix.shape))
+                    # Add the PSD matrix to the list of prepared PSDs
+                    self._list_of_prepared_psd.append(psd_matrix)
 
                 if verbose:
                     print('\t\t\tFunction', function_counter, ':', len(function.list_of_class_psd),
@@ -513,8 +493,7 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
                           'scalar constraint(s) ...')
 
                 for constraint in function.list_of_constraints:
-                    wrapper.send_constraint_to_solver(constraint)
-                    self._list_of_constraints_sent_to_wrapper.append(constraint)
+                    self._list_of_prepared_constraints.append(constraint)
 
                 if verbose:
                     print('\t\t\tFunction', function_counter, ':', len(function.list_of_constraints),
@@ -526,8 +505,11 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
                           'lmi constraint(s) ...')
 
                 for psd_counter, psd_matrix in enumerate(function.list_of_psd):
-                    wrapper.send_lmi_constraint_to_solver(psd_counter, psd_matrix)
-                    self._list_of_psd_sent_to_wrapper.append(psd_matrix)
+                    # Print a message if verbose mode activated
+                    if verbose > 0:
+                        print('\t\t Size of PSD matrix {}: {}x{}'.format(psd_counter + 1, *psd_matrix.shape))
+                    # Add the PSD matrix to the list of prepared PSDs
+                    self._list_of_prepared_psd.append(psd_matrix)
 
                 if verbose:
                     print('\t\t\tFunction', function_counter, ':', len(function.list_of_psd),
@@ -545,12 +527,77 @@ def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
                 print('\t\t\tPartition', partition_counter, 'with', partition.get_nb_blocks(),
                       'blocks: Adding', len(partition.list_of_constraints), 'scalar constraint(s)...')
             for constraint in partition.list_of_constraints:
-                wrapper.send_constraint_to_solver(constraint)
-                self._list_of_constraints_sent_to_wrapper.append(constraint)
+                self._list_of_prepared_constraints.append(constraint)
             if verbose:
                 print('\t\t\tPartition', partition_counter, 'with', partition.get_nb_blocks(),
                       'blocks:', len(partition.list_of_constraints), 'scalar constraint(s) added')
 
+    def _solve_with_wrapper(self, wrapper, verbose=1, return_primal_or_dual="dual",
+                            dimension_reduction_heuristic=None, eig_regularization=1e-3, tol_dimension_reduction=1e-4,
+                            **kwargs):
+        """
+        Internal solve method. Translate the :class:`PEP` to an SDP, and solve it via the wrapper.
+
+        Args:
+            wrapper (Wrapper): Interface to the solver.
+            verbose (int, optional): Level of information details to print
+                                     (Override the CVXPY solver verbose parameter).
+
+                            - 0: No verbose at all
+                            - 1: PEPit information is printed but not CVXPY's
+                            - 2: Both PEPit and solver details are printed
+            return_primal_or_dual (str, optional): If "dual", it returns a worst-case upper bound of the PEP
+                                                   (dual value of the objective).
+                                                   If "primal", it returns a worst-case lower bound of the PEP
+                                                   (primal value of the objective).
+                                                   Default is "dual".
+                                                   Note both value should be almost the same by strong duality.
+            dimension_reduction_heuristic (str, optional): An heuristic to reduce the dimension of the solution
+                                                           (rank of the Gram matrix). Set to None to deactivate
+                                                           it (default value). Available heuristics are:
+
+                                                            - "trace": minimize :math:`Tr(G)`
+                                                            - "logdet{an integer n}": minimize
+                                                              :math:`\\log\\left(\\mathrm{Det}(G)\\right)`
+                                                              using n iterations of local approximation problems.
+
+            eig_regularization (float, optional): The regularization we use to make
+                                                  :math:`G + \\mathrm{eig_regularization}I_d \succ 0`.
+                                                  (only used when "dimension_reduction_heuristic" is not None)
+                                                  The default value is 1e-5.
+            tol_dimension_reduction (float, optional): The error tolerance in the heuristic minimization problem.
+                                                       Precisely, the second problem minimizes "optimal_value - tol"
+                                                       (only used when "dimension_reduction_heuristic" is not None)
+                                                       The default value is 1e-5.
+            kwargs (keywords, optional): Additional solver-specific arguments.
+
+        Returns:
+            float: Worst-case guarantee of the PEP.
+
+        """
+
+        # Report the creation of variables (G, F)
+        if verbose:
+            print('(PEPit) Setting up the problem:'
+                  ' size of the Gram matrix: {}x{}'.format(Point.counter, Point.counter))
+        wrapper.set_main_variables()
+
+        # Determine the lists of constraints and LMIs sent to wrapper
+        self._list_of_constraints_sent_to_wrapper = [constraint for constraint in self._list_of_prepared_constraints if constraint.activated]
+        self._list_of_psd_sent_to_wrapper = [psd for psd in self._list_of_prepared_psd if psd.activated]
+        
+        # Clean dual variables before solving
+        for constraint in self._list_of_prepared_constraints:
+            constraint._dual_variable_value = 0.
+        for psd in self._list_of_prepared_psd:
+            psd._dual_variable_value = np.zeros_like(psd.matrix_of_expressions)
+
+        # Send constraints to the wrapper
+        for constraint in self._list_of_constraints_sent_to_wrapper:
+            wrapper.send_constraint_to_solver(constraint)
+        for psd in self._list_of_psd_sent_to_wrapper:
+            wrapper.send_lmi_constraint_to_solver(psd)
+
         # Instantiate the problem
         if verbose:
             print('(PEPit) Compiling SDP')
@@ -704,13 +751,13 @@ def check_feasibility(self, wc_value, verbose=1):
                 print(message)
 
         # Get the max value of all transgression of the constraints
-        if self._list_of_constraints_sent_to_wrapper:
+        if self._list_of_prepared_constraints:
             max_constraint_error = np.max(
                 [constraint.eval()
-                 for constraint in self._list_of_constraints_sent_to_wrapper
+                 for constraint in self._list_of_prepared_constraints
                  if constraint.equality_or_inequality == "inequality"]
                 + [np.abs(constraint.eval())
-                   for constraint in self._list_of_constraints_sent_to_wrapper
+                   for constraint in self._list_of_prepared_constraints
                    if constraint.equality_or_inequality == "equality"]
             )
             if verbose:
@@ -755,7 +802,7 @@ def check_feasibility(self, wc_value, verbose=1):
         # Scalar constraints
         # Dual of inequality constraints >= 0
         inequality_constraint_dual_values = [constraint.eval_dual()
-                                             for constraint in self._list_of_constraints_sent_to_wrapper
+                                             for constraint in self._list_of_prepared_constraints
                                              if constraint.equality_or_inequality == "inequality"]
         if inequality_constraint_dual_values:
             inequality_constraint_dual_min_value = np.min(inequality_constraint_dual_values)
@@ -765,7 +812,7 @@ def check_feasibility(self, wc_value, verbose=1):
                     message += " up to an error of {}".format(-inequality_constraint_dual_min_value)
                 print(message)
         # + <expression, dual> <= 0
-        for constraint in self._list_of_constraints_sent_to_wrapper:
+        for constraint in self._list_of_prepared_constraints:
             constraints_combination += constraint.eval_dual() * constraint.expression
 
         # Proof reconstruction
diff --git a/PEPit/psd_matrix.py b/PEPit/psd_matrix.py
index e3314a84..c319b372 100644
--- a/PEPit/psd_matrix.py
+++ b/PEPit/psd_matrix.py
@@ -9,6 +9,7 @@ class PSDMatrix(object):
 
     Attributes:
         name (str): A name set through the set_name method. None is no name is given.
+        activated (bool): A boolean flag used to activate/deactivate the LMI constraint in the PEP.
         matrix_of_expressions (Iterable of Iterable of Expression): a square matrix of :class:`Expression` objects.
         shape (tuple of ints): the shape of the underlying matrix of :class:`Expression` objects.
         _value (2D ndarray of floats): numerical values of :class:`Expression` objects
@@ -61,6 +62,9 @@ def __init__(self,
             TypeError: if provided matrix does not contain only Expressions and / or scalar values.
 
         """
+        # Initialize the activated attribute to True
+        self.activated = True
+
         # Initialize name of the psd matrix
         self.name = name
 
@@ -95,6 +99,20 @@ def get_name(self):
         """
         return self.name
 
+    def activate(self):
+        """
+        Activate the use of the LMI constraint.
+
+        """
+        self.activated = True
+
+    def deactivate(self):
+        """
+        Deactivate the use of the LMI constraint.
+
+        """
+        self.activated = False
+
     @staticmethod
     def _store(matrix_of_expressions):
         """
diff --git a/PEPit/wrapper.py b/PEPit/wrapper.py
index 7e224e4d..27697aa2 100644
--- a/PEPit/wrapper.py
+++ b/PEPit/wrapper.py
@@ -89,13 +89,12 @@ def send_constraint_to_solver(self, constraint):
         """
         raise NotImplementedError("This method must be overwritten in children classes")
 
-    def send_lmi_constraint_to_solver(self, psd_counter, psd_matrix):
+    def send_lmi_constraint_to_solver(self, psd_matrix):
         """
         Transfer a PEPit :class:`PSDMatrix` (LMI constraint) to the solver
         and add it the tracking lists.
 
         Args:
-            psd_counter (int): a counter useful for the verbose mode.
             psd_matrix (PSDMatrix): a matrix of expressions that is constrained to be PSD.
 
         """
diff --git a/PEPit/wrappers/cvxpy_wrapper.py b/PEPit/wrappers/cvxpy_wrapper.py
index 2605e1ff..5e248265 100644
--- a/PEPit/wrappers/cvxpy_wrapper.py
+++ b/PEPit/wrappers/cvxpy_wrapper.py
@@ -135,18 +135,17 @@ def send_constraint_to_solver(self, constraint):
             # Raise an exception otherwise
             raise ValueError('The attribute \'equality_or_inequality\' of a constraint object'
                              ' must either be \'equality\' or \'inequality\'.'
-                             'Got {}'.format(constraint.equality_or_inequality))
+                             '{} got {}'.format(constraint.get_name(), constraint.equality_or_inequality))
 
         # Add the corresponding CVXPY constraint to the list of constraints to be sent to CVXPY
         self._list_of_solver_constraints.append(cvxpy_constraint)
 
-    def send_lmi_constraint_to_solver(self, psd_counter, psd_matrix):
+    def send_lmi_constraint_to_solver(self, psd_matrix):
         """
         Transform a PEPit :class:`PSDMatrix` into a CVXPY symmetric PSD matrix
         and add the 2 formats of the constraints into the tracking lists.
 
         Args:
-            psd_counter (int): a counter useful for the verbose mode.
             psd_matrix (PSDMatrix): a matrix of expressions that is constrained to be PSD.
 
         """
@@ -170,10 +169,6 @@ def send_lmi_constraint_to_solver(self, psd_counter, psd_matrix):
             for j in range(psd_matrix.shape[1]):
                 cvxpy_constraints_list.append(M[i, j] == self._expression_to_solver(psd_matrix[i, j]))
 
-        # Print a message if verbose mode activated
-        if self.verbose > 0:
-            print('\t\t Size of PSD matrix {}: {}x{}'.format(psd_counter + 1, *psd_matrix.shape))
-
         # Add the corresponding CVXPY constraints to the list of constraints to be sent to CVXPY
         self._list_of_solver_constraints += cvxpy_constraints_list
 
diff --git a/PEPit/wrappers/mosek_wrapper.py b/PEPit/wrappers/mosek_wrapper.py
index 1cc2427d..36c495a6 100644
--- a/PEPit/wrappers/mosek_wrapper.py
+++ b/PEPit/wrappers/mosek_wrapper.py
@@ -159,14 +159,13 @@ def send_constraint_to_solver(self, constraint, track=True):
             # Raise an exception otherwise
             raise ValueError('The attribute \'equality_or_inequality\' of a constraint object'
                              ' must either be \'equality\' or \'inequality\'.'
-                             'Got {}'.format(constraint.equality_or_inequality))
+                             '{} got {}'.format(constraint.get_name(), constraint.equality_or_inequality))
 
-    def send_lmi_constraint_to_solver(self, psd_counter, psd_matrix):
+    def send_lmi_constraint_to_solver(self, psd_matrix):
         """
         Transfer a PEPit :class:`PSDMatrix` (LMI constraint) to MOSEK and add it the tracking lists.
 
         Args:
-            psd_counter (int): a counter useful for the verbose mode.
             psd_matrix (PSDMatrix): a matrix of expressions that is constrained to be PSD.
 
         """
@@ -204,10 +203,6 @@ def send_lmi_constraint_to_solver(self, psd_counter, psd_matrix):
                 self.task.putaijlist(np.full(a_i.shape, nb_cons), a_i, a_val)
                 self.task.putconbound(nb_cons, mosek.boundkey.fx, -alpha_val, -alpha_val)
 
-        # Print a message if verbose mode activated
-        if self.verbose > 0:
-            print('\t\t Size of PSD matrix {}: {}x{}'.format(psd_counter + 1, *psd_matrix.shape))
-
     def _recover_dual_values(self):
         """
         Recover all dual variables from solver.
diff --git a/tests/test_constraints.py b/tests/test_constraints.py
index e286ebb8..9eaccca8 100644
--- a/tests/test_constraints.py
+++ b/tests/test_constraints.py
@@ -74,7 +74,7 @@ def test_counter(self):
 
     def test_name(self):
 
-        self.assertIsNone(self.initial_condition.get_name())
+        self.assertEqual(self.initial_condition.get_name(), "Constraint 0")
         self.assertIsNone(self.performance_metric.get_name())
 
         self.initial_condition.set_name("init")
@@ -93,6 +93,13 @@ def test_equality_inequality(self):
             self.assertIsInstance(self.problem.list_of_constraints[i].equality_or_inequality, str)
             self.assertIn(self.problem.list_of_constraints[i].equality_or_inequality, {'equality', 'inequality'})
 
+    def test_activated(self):
+        self.assertIs(self.initial_condition.activated, True)
+        self.initial_condition.deactivate()
+        self.assertIs(self.initial_condition.activated, False)
+        self.initial_condition.activate()
+        self.assertIs(self.initial_condition.activated, True)
+
     def test_eval(self):
 
         for i in range(len(self.func.list_of_class_constraints)):
diff --git a/tests/test_pep.py b/tests/test_pep.py
index ab24ba3b..51b14ed6 100644
--- a/tests/test_pep.py
+++ b/tests/test_pep.py
@@ -160,12 +160,40 @@ def test_lmi_constraints_in_real_problem(self):
     def test_consistency(self):
 
         # Solve twice the same problem in a row and verify the two lists of constraints have same length.
+        self.assertEqual(len(self.problem._list_of_prepared_constraints), 0)
+        self.assertEqual(len(self.problem._list_of_constraints_sent_to_wrapper), 0)
+
+        # Run solve
         _ = self.problem.solve(verbose=self.verbose)
-        l1 = self.problem._list_of_constraints_sent_to_wrapper
+        lp1 = self.problem._list_of_prepared_constraints
+        ls1 = self.problem._list_of_constraints_sent_to_wrapper
+
+        # Rerun solve
         _ = self.problem.solve(verbose=self.verbose)
-        l2 = self.problem._list_of_constraints_sent_to_wrapper
+        lp2 = self.problem._list_of_prepared_constraints
+        ls2 = self.problem._list_of_constraints_sent_to_wrapper
 
-        self.assertEqual(len(l1), len(l2))
+        # Deactivate one constraint, then rerun solve
+        self.problem.list_of_constraints[0].deactivate()
+        _ = self.problem.solve(verbose=self.verbose)
+        lp3 = self.problem._list_of_prepared_constraints
+        ls3 = self.problem._list_of_constraints_sent_to_wrapper
+
+        # Reactivate the constraint, then rerun solve
+        self.problem.list_of_constraints[0].activate()
+        _ = self.problem.solve(verbose=self.verbose)
+        lp4 = self.problem._list_of_prepared_constraints
+        ls4 = self.problem._list_of_constraints_sent_to_wrapper
+
+        # Compare lengths
+        self.assertEqual(len(lp1), len(lp2))
+        self.assertEqual(len(lp1), len(lp3))
+        self.assertEqual(len(lp1), len(lp4))
+
+        self.assertEqual(len(lp1), len(ls1))
+        self.assertEqual(len(lp2), len(ls2))
+        self.assertEqual(len(lp3), len(ls3)+1)
+        self.assertEqual(len(lp4), len(ls4))
 
     def test_dimension_reduction(self):
 
diff --git a/tests/test_psd_matrix.py b/tests/test_psd_matrix.py
index 26e7b5e8..0a97cd69 100644
--- a/tests/test_psd_matrix.py
+++ b/tests/test_psd_matrix.py
@@ -103,6 +103,13 @@ def test_getitem(self):
         self.assertIs(self.psd1[0, 0], self.expr1)
         self.assertIs(self.psd2[0, 1], self.expr3)
 
+    def test_activated(self):
+        self.assertIs(self.psd1.activated, True)
+        self.psd1.deactivate()
+        self.assertIs(self.psd1.activated, False)
+        self.psd1.activate()
+        self.assertIs(self.psd1.activated, True)
+
     def test_eval(self):
         
         # The PEP has not been solved yet, so no value is accessible.
diff --git a/tests/test_wrappers.py b/tests/test_wrappers.py
index f8a1ee80..5473312a 100644
--- a/tests/test_wrappers.py
+++ b/tests/test_wrappers.py
@@ -112,6 +112,7 @@ def test_dual_sign_in_equality_constraints(self):
                            1 == (self.x0 - self.xs) ** 2,
                            ]:
             self.problem.list_of_constraints = [constraint]
+            self.problem._prepare_constraints(verbose=self.verbose)
             self.problem.solve(verbose=self.verbose, wrapper=self.wrapper)
             elements_of_proof.append(constraint.eval_dual() * constraint.expression)