diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel
index 4ecad430ff..a8477d00e0 100644
--- a/ortools/sat/BUILD.bazel
+++ b/ortools/sat/BUILD.bazel
@@ -149,7 +149,9 @@ cc_library(
         ":sat_base",
         ":sat_parameters_cc_proto",
         ":sat_solver",
+        ":symmetry_util",
         ":util",
+        "//ortools/algorithms:sparse_permutation",
         "//ortools/base",
         "//ortools/base:status_macros",
         "//ortools/base:stl_util",
@@ -2133,6 +2135,7 @@ cc_library(
         ":integer_expr",
         ":linear_constraint",
         ":linear_constraint_manager",
+        ":linear_propagation",
         ":model",
         ":sat_base",
         ":sat_parameters_cc_proto",
@@ -2153,7 +2156,6 @@ cc_library(
         "//ortools/lp_data:lp_data_utils",
         "//ortools/lp_data:scattered_vector",
         "//ortools/lp_data:sparse",
-        "//ortools/lp_data:sparse_column",
         "//ortools/util:bitset",
         "//ortools/util:rev",
         "//ortools/util:saturated_arithmetic",
@@ -2161,6 +2163,7 @@ cc_library(
         "//ortools/util:time_limit",
         "@com_google_absl//absl/algorithm:container",
         "@com_google_absl//absl/container:flat_hash_map",
+        "@com_google_absl//absl/log",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/numeric:int128",
         "@com_google_absl//absl/strings",
diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto
index b57d2056ba..92e61fac35 100644
--- a/ortools/sat/cp_model.proto
+++ b/ortools/sat/cp_model.proto
@@ -105,9 +105,13 @@ message LinearConstraintProto {
 // The constraint target = vars[index].
 // This enforces that index takes one of the value in [0, vars_size()).
 message ElementConstraintProto {
-  int32 index = 1;
-  int32 target = 2;
-  repeated int32 vars = 3;
+  int32 index = 1;          // Legacy field.
+  int32 target = 2;         // Legacy field.
+  repeated int32 vars = 3;  // Legacy field.
+  // All expressions below must be affine function with at most one variable.
+  LinearExpressionProto linear_index = 4;
+  LinearExpressionProto linear_target = 5;
+  repeated LinearExpressionProto exprs = 6;
 }
 
 // This is not really a constraint. It is there so it can be referred by other
diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc
index 8645013cd0..576537f0fb 100644
--- a/ortools/sat/cp_model_checker.cc
+++ b/ortools/sat/cp_model_checker.cc
@@ -418,23 +418,60 @@ std::string ValidateElementConstraint(const CpModelProto& model,
                                       const ConstraintProto& ct) {
   const ElementConstraintProto& element = ct.element();
 
+  if ((element.has_linear_index() || element.has_linear_target() ||
+       !element.exprs().empty()) &&
+      (!element.vars().empty() || element.index() != 0 ||
+       element.target() != 0)) {
+    return absl::StrCat(
+        "Inconsistent element with both legacy and new format defined",
+        ProtobufShortDebugString(ct));
+  }
+
   // We need to be able to manipulate expression like "target - var" without
   // integer overflow.
-  LinearExpressionProto overflow_detection;
-  overflow_detection.add_vars(element.target());
-  overflow_detection.add_coeffs(1);
-  overflow_detection.add_vars(/*dummy*/ 0);
-  overflow_detection.add_coeffs(-1);
-  for (const int ref : element.vars()) {
-    overflow_detection.set_vars(1, ref);
-    if (PossibleIntegerOverflow(model, overflow_detection.vars(),
-                                overflow_detection.coeffs())) {
-      return absl::StrCat(
-          "Domain of the variables involved in element constraint may cause "
-          "overflow",
-          ProtobufShortDebugString(ct));
+  if (!element.vars().empty()) {
+    LinearExpressionProto overflow_detection;
+    overflow_detection.add_vars(element.target());
+    overflow_detection.add_coeffs(1);
+    overflow_detection.add_vars(/*dummy*/ 0);
+    overflow_detection.add_coeffs(-1);
+    for (const int ref : element.vars()) {
+      overflow_detection.set_vars(1, ref);
+      if (PossibleIntegerOverflow(model, overflow_detection.vars(),
+                                  overflow_detection.coeffs())) {
+        return absl::StrCat(
+            "Domain of the variables involved in element constraint may cause "
+            "overflow",
+            ProtobufShortDebugString(ct));
+      }
     }
   }
+
+  if (!element.exprs().empty() || element.has_linear_index() ||
+      element.has_linear_target()) {
+    RETURN_IF_NOT_EMPTY(
+        ValidateAffineExpression(model, element.linear_index()));
+    RETURN_IF_NOT_EMPTY(
+        ValidateAffineExpression(model, element.linear_target()));
+    for (const LinearExpressionProto& expr : element.exprs()) {
+      RETURN_IF_NOT_EMPTY(ValidateAffineExpression(model, expr));
+      LinearExpressionProto overflow_detection = ct.element().linear_target();
+      for (int i = 0; i < expr.vars_size(); ++i) {
+        overflow_detection.add_vars(expr.vars(i));
+        overflow_detection.add_coeffs(-expr.coeffs(i));
+      }
+      overflow_detection.set_offset(overflow_detection.offset() -
+                                    expr.offset());
+      if (PossibleIntegerOverflow(model, overflow_detection.vars(),
+                                  overflow_detection.coeffs(),
+                                  overflow_detection.offset())) {
+        return absl::StrCat(
+            "Domain of the variables involved in element constraint may cause "
+            "overflow");
+      }
+    }
+  }
+
   return "";
 }
 
@@ -1397,10 +1434,20 @@ class ConstraintChecker {
   }
 
   bool ElementConstraintIsFeasible(const ConstraintProto& ct) {
-    if (ct.element().vars().empty()) return false;
-    const int index = Value(ct.element().index());
-    if (index < 0 || index >= ct.element().vars_size()) return false;
-    return Value(ct.element().vars(index)) == Value(ct.element().target());
+    if (!ct.element().vars().empty()) {
+      const int index = Value(ct.element().index());
+      if (index < 0 || index >= ct.element().vars_size()) return false;
+      return Value(ct.element().vars(index)) == Value(ct.element().target());
+    }
+
+    if (!ct.element().exprs().empty()) {
+      const int index = LinearExpressionValue(ct.element().linear_index());
+      if (index < 0 || index >= ct.element().exprs_size()) return false;
+      return LinearExpressionValue(ct.element().exprs(index)) ==
+             LinearExpressionValue(ct.element().linear_target());
+    }
+
+    return false;
   }
 
   bool TableConstraintIsFeasible(const ConstraintProto& ct) {
diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc
index eaa97cdc80..cd7efad851 100644
--- a/ortools/sat/cp_model_expand.cc
+++ b/ortools/sat/cp_model_expand.cc
@@ -784,94 +784,81 @@ void ExpandLinMax(ConstraintProto* ct, PresolveContext* context) {
 }
 
 // A[V] == V means for all i, V == i => A_i == i
-void ExpandElementWithTargetEqualIndex(ConstraintProto* ct,
-                                       PresolveContext* context) {
+void ExpandElementWhenTargetShareVarWithIndex(ConstraintProto* ct,
+                                              PresolveContext* context) {
   const ElementConstraintProto& element = ct->element();
-  DCHECK_EQ(element.index(), element.target());
-
-  const int index_ref = element.index();
-  std::vector<int64_t> valid_indices;
-  for (const int64_t v : context->DomainOf(index_ref).Values()) {
-    if (!context->DomainContains(element.vars(v), v)) continue;
-    valid_indices.push_back(v);
-  }
-  if (valid_indices.size() < context->DomainOf(index_ref).Size()) {
-    if (!context->IntersectDomainWith(index_ref,
-                                      Domain::FromValues(valid_indices))) {
-      VLOG(1) << "No compatible variable domains in "
-                 "ExpandElementWithTargetEqualIndex()";
-      return;
-    }
-    context->UpdateRuleStats("element: reduced index domain");
-  }
 
-  for (const int64_t v : context->DomainOf(index_ref).Values()) {
-    const int var = element.vars(v);
-    if (context->MinOf(var) == v && context->MaxOf(var) == v) continue;
-    context->AddImplyInDomain(
-        context->GetOrCreateVarValueEncoding(index_ref, v), var, Domain(v));
+  const LinearExpressionProto& index = element.linear_index();
+  DCHECK_EQ(index.vars_size(), 1);
+  const int index_var = index.vars(0);
+  const LinearExpressionProto& target = element.linear_target();
+  DCHECK_EQ(target.vars_size(), 1);
+  DCHECK_EQ(target.vars(0), index_var);
+
+  for (const int64_t v : context->DomainOf(index_var).Values()) {
+    const int64_t index_value = AffineExpressionValueAt(index, v);
+    const int64_t target_value = AffineExpressionValueAt(target, v);
+    const LinearExpressionProto& expr = element.exprs(index_value);
+    ConstraintProto* imply = context->working_model->add_constraints();
+    imply->add_enforcement_literal(
+        context->GetOrCreateVarValueEncoding(index_var, v));
+    imply->mutable_linear()->add_domain(target_value);
+    imply->mutable_linear()->add_domain(target_value);
+    AddLinearExpressionToLinearConstraint(expr, 1, imply->mutable_linear());
+    context->CanonicalizeLinearConstraint(imply);
   }
+
   context->UpdateRuleStats(
-      "element: expanded with special case target = index");
+      "element: expanded when the index and the target share the same var");
   ct->Clear();
 }
 
 // Special case if the array of the element is filled with constant values.
 void ExpandConstantArrayElement(ConstraintProto* ct, PresolveContext* context) {
   const ElementConstraintProto& element = ct->element();
-  const int index_ref = element.index();
-  const int target_ref = element.target();
+  const LinearExpressionProto& index = element.linear_index();
+  DCHECK_EQ(index.vars_size(), 1);
+  const int index_var = index.vars(0);
+  const LinearExpressionProto& target = element.linear_target();
 
   // Index and target domain have been reduced before calling this function.
-  const Domain index_domain = context->DomainOf(index_ref);
-  const Domain target_domain = context->DomainOf(target_ref);
+  const Domain index_var_domain = context->DomainOf(index_var);
+  const Domain target_domain = context->DomainSuperSetOf(target);
 
   // This BoolOrs implements the deduction that if all index literals pointing
   // to the same value in the constant array are false, then this value is no
   // no longer valid for the target variable. They are created only for values
   // that have multiples literals supporting them.
-  // Order is not important.
-  absl::flat_hash_map<int64_t, BoolArgumentProto*> supports;
-  {
-    absl::flat_hash_map<int64_t, int> constant_var_values_usage;
-    for (const int64_t v : index_domain.Values()) {
-      DCHECK(context->IsFixed(element.vars(v)));
-      const int64_t value = context->MinOf(element.vars(v));
-      if (++constant_var_values_usage[value] == 2) {
-        // First time we cross > 1.
-        BoolArgumentProto* const support =
-            context->working_model->add_constraints()->mutable_bool_or();
-        const int target_literal =
-            context->GetOrCreateVarValueEncoding(target_ref, value);
-        support->add_literals(NegatedRef(target_literal));
-        supports[value] = support;
-      }
-    }
+  absl::btree_map<int64_t, std::vector<int>> supports;
+  for (const int64_t v : index_var_domain.Values()) {
+    const int64_t index_value = AffineExpressionValueAt(index, v);
+    const int64_t expr_value = context->FixedValue(element.exprs(index_value));
+    supports[expr_value].push_back(v);
   }
 
-  {
-    // While this is not strictly needed since all value in the index will be
-    // covered, it allows to easily detect this fact in the presolve.
-    auto* exactly_one =
-        context->working_model->add_constraints()->mutable_exactly_one();
-    for (const int64_t v : index_domain.Values()) {
+  // While this is not strictly needed since all value in the index will be
+  // covered, it allows to easily detect this fact in the presolve.
+  //
+  // TODO(user): Do we need the support part ? Is this discovered by probing?
+  auto* exactly_one =
+      context->working_model->add_constraints()->mutable_exactly_one();
+  for (const auto& [expr_value, support] : supports) {
+    const int target_literal =
+        context->GetOrCreateAffineValueEncoding(target, expr_value);
+    // not(indices supporting value) -> target != value
+    BoolArgumentProto* bool_or =
+        context->working_model->add_constraints()->mutable_bool_or();
+    bool_or->add_literals(NegatedRef(target_literal));
+    for (const int64_t v : support) {
       const int index_literal =
-          context->GetOrCreateVarValueEncoding(index_ref, v);
-      exactly_one->add_literals(index_literal);
+          context->GetOrCreateVarValueEncoding(index_var, v);
+      bool_or->add_literals(index_literal);
 
-      const int64_t value = context->MinOf(element.vars(v));
-      const auto& it = supports.find(value);
-      if (it != supports.end()) {
-        // The encoding literal for 'value' of the target_ref has been
-        // created before.
-        const int target_literal =
-            context->GetOrCreateVarValueEncoding(target_ref, value);
-        context->AddImplication(index_literal, target_literal);
-        it->second->add_literals(index_literal);
-      } else {
-        // Try to reuse the literal of the index.
-        context->InsertVarValueEncoding(index_literal, target_ref, value);
-      }
+      // index == v => target == value
+      context->AddImplication(index_literal, target_literal);
+
+      // Helps presolve.
+      exactly_one->add_literals(index_literal);
     }
   }
 
@@ -882,51 +869,36 @@ void ExpandConstantArrayElement(ConstraintProto* ct, PresolveContext* context) {
 // General element when the array contains non fixed variables.
 void ExpandVariableElement(ConstraintProto* ct, PresolveContext* context) {
   const ElementConstraintProto& element = ct->element();
-  const int index_ref = element.index();
-  const int target_ref = element.target();
-  const Domain index_domain = context->DomainOf(index_ref);
+  const LinearExpressionProto& index = element.linear_index();
+  DCHECK_EQ(index.vars_size(), 1);
+  const int index_var = index.vars(0);
+  const Domain index_var_domain = context->DomainOf(index_var);
+  const LinearExpressionProto& target = element.linear_target();
 
   BoolArgumentProto* exactly_one =
       context->working_model->add_constraints()->mutable_exactly_one();
 
-  for (const int64_t v : index_domain.Values()) {
-    const int var = element.vars(v);
-    const Domain var_domain = context->DomainOf(var);
-    const int index_lit = context->GetOrCreateVarValueEncoding(index_ref, v);
+  for (const int64_t v : index_var_domain.Values()) {
+    const int64_t index_value = AffineExpressionValueAt(index, v);
+    DCHECK_GE(index_value, 0);
+    DCHECK_LT(index_value, element.exprs_size());
+    const int index_lit = context->GetOrCreateVarValueEncoding(index_var, v);
     exactly_one->add_literals(index_lit);
 
-    if (var_domain.IsFixed()) {
-      context->AddImplyInDomain(index_lit, target_ref, var_domain);
-    } else {
-      // We make sure we only use positive ref.
-      //
-      // TODO(user): Get rid of this code once we accept affine in element
-      // constraint.
-      ConstraintProto* const ct = context->working_model->add_constraints();
-      ct->add_enforcement_literal(index_lit);
-      if (RefIsPositive(var)) {
-        ct->mutable_linear()->add_vars(var);
-        ct->mutable_linear()->add_coeffs(1);
-      } else {
-        ct->mutable_linear()->add_vars(NegatedRef(var));
-        ct->mutable_linear()->add_coeffs(-1);
-      }
-      if (RefIsPositive(target_ref)) {
-        ct->mutable_linear()->add_vars(target_ref);
-        ct->mutable_linear()->add_coeffs(-1);
-      } else {
-        ct->mutable_linear()->add_vars(NegatedRef(target_ref));
-        ct->mutable_linear()->add_coeffs(1);
-      }
-      ct->mutable_linear()->add_domain(0);
-      ct->mutable_linear()->add_domain(0);
+    ConstraintProto* const imply = context->working_model->add_constraints();
+    imply->add_enforcement_literal(index_lit);
+    imply->mutable_linear()->add_domain(0);
+    imply->mutable_linear()->add_domain(0);
+    AddLinearExpressionToLinearConstraint(target, -1, imply->mutable_linear());
+    AddLinearExpressionToLinearConstraint(ct->element().exprs(index_value), 1,
+                                          imply->mutable_linear());
+    context->CanonicalizeLinearConstraint(imply);
 
-      // Note that this should have been checked at model validation.
-      DCHECK(!PossibleIntegerOverflow(*context->working_model,
-                                      ct->mutable_linear()->vars(),
-                                      ct->mutable_linear()->coeffs()))
-          << google::protobuf::ShortFormat(*ct);
-    }
+    // Note that this should have been checked at model validation.
+    DCHECK(!PossibleIntegerOverflow(*context->working_model,
+                                    imply->mutable_linear()->vars(),
+                                    imply->mutable_linear()->coeffs()))
+        << google::protobuf::ShortFormat(*imply);
   }
 
   context->UpdateRuleStats("element: expanded");
@@ -936,68 +908,53 @@ void ExpandVariableElement(ConstraintProto* ct, PresolveContext* context) {
 void ExpandElement(ConstraintProto* ct, PresolveContext* context) {
   const ElementConstraintProto& element = ct->element();
 
-  const int index_ref = element.index();
-  const int target_ref = element.target();
-  const int size = element.vars_size();
+  const LinearExpressionProto& index = element.linear_index();
+  const LinearExpressionProto& target = element.linear_target();
+  const int size = element.exprs_size();
 
   // Reduce the domain of the index to be compatible with the array of
   // variables. Note that the element constraint is 0 based.
-  if (!context->IntersectDomainWith(index_ref, Domain(0, size - 1))) {
+  if (!context->IntersectDomainWith(index, Domain(0, size - 1))) {
     VLOG(1) << "Empty domain for the index variable in ExpandElement()";
     return;
   }
 
-  // Special case when index = target.
-  if (index_ref == target_ref) {
-    ExpandElementWithTargetEqualIndex(ct, context);
+  if (context->IsFixed(index)) {
+    ConstraintProto* const eq = context->working_model->add_constraints();
+    eq->mutable_linear()->add_domain(0);
+    eq->mutable_linear()->add_domain(0);
+    AddLinearExpressionToLinearConstraint(target, 1, eq->mutable_linear());
+    AddLinearExpressionToLinearConstraint(
+        ct->element().exprs(context->FixedValue(index)), -1,
+        eq->mutable_linear());
+    context->CanonicalizeLinearConstraint(eq);
+    context->UpdateRuleStats("element: expanded with fixed index");
+    ct->Clear();
     return;
   }
 
-  // Reduces the domain of the index and the target.
-  bool all_constants = true;
-  std::vector<int64_t> valid_indices;
-  const Domain index_domain = context->DomainOf(index_ref);
-  const Domain target_domain = context->DomainOf(target_ref);
-  Domain reached_domain;
-  for (const int64_t v : index_domain.Values()) {
-    const Domain var_domain = context->DomainOf(element.vars(v));
-    if (var_domain.IntersectionWith(target_domain).IsEmpty()) continue;
-
-    valid_indices.push_back(v);
-    reached_domain = reached_domain.UnionWith(var_domain);
-    if (var_domain.Min() != var_domain.Max()) {
-      all_constants = false;
-    }
-  }
-
-  if (valid_indices.size() < index_domain.Size()) {
-    if (!context->IntersectDomainWith(index_ref,
-                                      Domain::FromValues(valid_indices))) {
-      VLOG(1) << "No compatible variable domains in ExpandElement()";
-      return;
-    }
-
-    context->UpdateRuleStats("element: reduced index domain");
-  }
-
-  // We know the target_domain is not empty as this would have triggered the
-  // above check.
-  bool target_domain_changed = false;
-  if (!context->IntersectDomainWith(target_ref, reached_domain,
-                                    &target_domain_changed)) {
+  // Special case when index.var = target.var.
+  if (index.vars_size() == 1 && target.vars_size() == 1 &&
+      index.vars(0) == target.vars(0)) {
+    ExpandElementWhenTargetShareVarWithIndex(ct, context);
     return;
   }
 
-  if (target_domain_changed) {
-    context->UpdateRuleStats("element: reduced target domain");
+  // Checks if all elements are constant.
+  bool all_constants = true;
+  for (const int64_t v : context->DomainOf(index.vars(0)).Values()) {
+    const int64_t index_value = AffineExpressionValueAt(index, v);
+    if (!context->IsFixed(element.exprs(index_value))) {
+      all_constants = false;
+      break;
+    }
   }
 
   if (all_constants) {
     ExpandConstantArrayElement(ct, context);
-    return;
+  } else {
+    ExpandVariableElement(ct, context);
   }
-
-  ExpandVariableElement(ct, context);
 }
 
 // Adds clauses so that literals[i] true <=> encoding[values[i]] true.
diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc
index eb135cfc0e..b6b13994c3 100644
--- a/ortools/sat/cp_model_lns.cc
+++ b/ortools/sat/cp_model_lns.cc
@@ -933,8 +933,7 @@ NeighborhoodGeneratorHelper::GetSchedulingPrecedences(
   return result;
 }
 
-std::vector<std::vector<int>>
-NeighborhoodGeneratorHelper::GetRoutingPathVariables(
+std::vector<std::vector<int>> NeighborhoodGeneratorHelper::GetRoutingPaths(
     const CpSolverResponse& initial_solution) const {
   struct HeadAndArcLiteral {
     int head;
@@ -1919,9 +1918,11 @@ Neighborhood LocalBranchingLpBasedNeighborhoodGenerator::Generate(
   // TODO(user): Does the current solution can provide a warm-start for the
   // LP?
   auto* response_manager = model.GetOrCreate<SharedResponseManager>();
-  response_manager->InitializeObjective(local_cp_model);
-  LoadCpModel(local_cp_model, &model);
-  SolveLoadedCpModel(local_cp_model, &model);
+  {
+    response_manager->InitializeObjective(local_cp_model);
+    LoadCpModel(local_cp_model, &model);
+    SolveLoadedCpModel(local_cp_model, &model);
+  }
 
   // Update dtime.
   data.deterministic_time +=
@@ -2455,24 +2456,26 @@ Neighborhood RoutingRandomNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
   const std::vector<std::vector<int>> all_paths =
-      helper_.GetRoutingPathVariables(initial_solution);
+      helper_.GetRoutingPaths(initial_solution);
 
   // Collect all unique variables.
-  std::vector<int> variables_to_fix;
-  for (const auto& path : all_paths) {
-    variables_to_fix.insert(variables_to_fix.end(), path.begin(), path.end());
+  absl::flat_hash_set<int> all_path_variables;
+  for (auto& path : all_paths) {
+    all_path_variables.insert(path.begin(), path.end());
   }
-  gtl::STLSortAndRemoveDuplicates(&variables_to_fix);
-  GetRandomSubset(1.0 - data.difficulty, &variables_to_fix, random);
+  std::vector<int> fixed_variables(all_path_variables.begin(),
+                                   all_path_variables.end());
+  std::sort(fixed_variables.begin(), fixed_variables.end());
+  GetRandomSubset(1.0 - data.difficulty, &fixed_variables, random);
   return helper_.FixGivenVariables(
-      initial_solution, {variables_to_fix.begin(), variables_to_fix.end()});
+      initial_solution, {fixed_variables.begin(), fixed_variables.end()});
 }
 
 Neighborhood RoutingPathNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
   std::vector<std::vector<int>> all_paths =
-      helper_.GetRoutingPathVariables(initial_solution);
+      helper_.GetRoutingPaths(initial_solution);
 
   // Collect all unique variables.
   absl::flat_hash_set<int> all_path_variables;
@@ -2518,7 +2521,7 @@ Neighborhood RoutingFullPathNeighborhoodGenerator::Generate(
     const CpSolverResponse& initial_solution, SolveData& data,
     absl::BitGenRef random) {
   std::vector<std::vector<int>> all_paths =
-      helper_.GetRoutingPathVariables(initial_solution);
+      helper_.GetRoutingPaths(initial_solution);
   // Remove a corner case where all paths are empty.
   if (all_paths.empty()) {
     return helper_.NoNeighborhood();
@@ -2580,32 +2583,6 @@ Neighborhood RoutingFullPathNeighborhoodGenerator::Generate(
   return helper_.FixGivenVariables(initial_solution, fixed_variables);
 }
 
-Neighborhood RoutingStartsNeighborhoodGenerator::Generate(
-    const CpSolverResponse& initial_solution, SolveData& data,
-    absl::BitGenRef random) {
-  std::vector<std::vector<int>> all_paths =
-      helper_.GetRoutingPathVariables(initial_solution);
-  // TODO(user): Maybe enable some routes to be non empty?
-  if (all_paths.empty()) return helper_.NoNeighborhood();
-
-  // Collect all unique path variables, except the start and end of the paths.
-  // For circuit constraints, we approximate this with the arc going in and out
-  // of the smallest node. This works well for model where the zero node is an
-  // artificial node.
-  std::vector<int> variables_to_fix;
-  for (const auto& path : all_paths) {
-    for (const int ref : path) {
-      if (ref == path.front() || ref == path.back()) continue;
-      variables_to_fix.push_back(PositiveRef(ref));
-    }
-  }
-  gtl::STLSortAndRemoveDuplicates(&variables_to_fix);
-  GetRandomSubset(1.0 - data.difficulty, &variables_to_fix, random);
-
-  return helper_.FixGivenVariables(
-      initial_solution, {variables_to_fix.begin(), variables_to_fix.end()});
-}
-
 bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const {
   return (incomplete_solutions_->HasSolution() ||
           lp_solutions_->NumSolutions() > 0);
diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h
index ac69d2562b..7c2a096c4f 100644
--- a/ortools/sat/cp_model_lns.h
+++ b/ortools/sat/cp_model_lns.h
@@ -31,6 +31,7 @@
 #include "absl/types/span.h"
 #include "google/protobuf/arena.h"
 #include "ortools/sat/cp_model.pb.h"
+#include "ortools/sat/cp_model_solver_helpers.h"
 #include "ortools/sat/integer.h"
 #include "ortools/sat/model.h"
 #include "ortools/sat/sat_parameters.pb.h"
@@ -238,7 +239,7 @@ class NeighborhoodGeneratorHelper : public SubSolver {
   // self-looping arcs. Path are sorted, starting from the arc with the lowest
   // tail index, and going in sequence up to the last arc before the circuit is
   // closed. Each entry correspond to the arc literal on the circuit.
-  std::vector<std::vector<int>> GetRoutingPathVariables(
+  std::vector<std::vector<int>> GetRoutingPaths(
       const CpSolverResponse& initial_solution) const;
 
   // Returns all precedences extracted from the scheduling constraint and the
@@ -609,9 +610,10 @@ class LocalBranchingLpBasedNeighborhoodGenerator
  public:
   LocalBranchingLpBasedNeighborhoodGenerator(
       NeighborhoodGeneratorHelper const* helper, absl::string_view name,
-      ModelSharedTimeLimit* const global_time_limit)
+      ModelSharedTimeLimit* const global_time_limit, SharedClasses* shared)
       : NeighborhoodGenerator(name, helper),
-        global_time_limit_(global_time_limit) {
+        global_time_limit_(global_time_limit),
+        shared_(shared) {
     // Given that we spend time generating a good neighborhood it sounds
     // reasonable to spend a bit more time solving it too.
     deterministic_limit_ = 0.5;
@@ -622,6 +624,7 @@ class LocalBranchingLpBasedNeighborhoodGenerator
 
  private:
   ModelSharedTimeLimit* const global_time_limit_;
+  SharedClasses* const shared_;
 };
 
 // Helper method for the scheduling neighborhood generators. Returns a
@@ -785,7 +788,7 @@ class RoutingPathNeighborhoodGenerator : public NeighborhoodGenerator {
                         SolveData& data, absl::BitGenRef random) final;
 };
 
-// This routing based LNS generator aims at relaxing one full path, and make
+// This routing based LNS generator aims are relaxing one full path, and make
 // some room on the other paths to absorb the nodes of the relaxed path.
 //
 // In order to do so, it will relax the first and the last arc of each path in
@@ -802,19 +805,6 @@ class RoutingFullPathNeighborhoodGenerator : public NeighborhoodGenerator {
                         SolveData& data, absl::BitGenRef random) final;
 };
 
-// This routing based LNS generator performs like the
-// RoutingRandomNeighborhoodGenerator, but always relax the arcs going in and
-// out of the depot for routes constraints, and of the node with the minimal
-// index for circuit constraints.
-class RoutingStartsNeighborhoodGenerator : public NeighborhoodGenerator {
- public:
-  RoutingStartsNeighborhoodGenerator(NeighborhoodGeneratorHelper const* helper,
-                                     absl::string_view name)
-      : NeighborhoodGenerator(name, helper) {}
-  Neighborhood Generate(const CpSolverResponse& initial_solution,
-                        SolveData& data, absl::BitGenRef random) final;
-};
-
 // Generates a neighborhood by fixing the variables to solutions reported in
 // various repositories. This is inspired from RINS published in "Exploring
 // relaxation induced neighborhoods to improve MIP solutions" 2004 by E. Danna
diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc
index bca33c8ab2..bd046fd112 100644
--- a/ortools/sat/cp_model_postsolve.cc
+++ b/ortools/sat/cp_model_postsolve.cc
@@ -208,6 +208,14 @@ int64_t EvaluateLinearExpression(const LinearExpressionProto& expr,
   return value;
 }
 
+bool LinearExpressionIsFixed(const LinearExpressionProto& expr,
+                             const std::vector<Domain>& domains) {
+  for (const int var : expr.vars()) {
+    if (!domains[var].IsFixed()) return false;
+  }
+  return true;
+}
+
 }  // namespace
 
 // Compute the max of each expression, and assign it to the target expr. We only
@@ -231,88 +239,48 @@ void PostsolveLinMax(const ConstraintProto& ct, std::vector<Domain>* domains) {
   (*domains)[target.vars(0)] = Domain(max_value);
 }
 
-// We only support 3 cases in the presolve currently.
+// We only support 2 cases:  either the index was removed, of the target, not
+// both.
 void PostsolveElement(const ConstraintProto& ct, std::vector<Domain>* domains) {
-  const int index_ref = ct.element().index();
-  const int index_var = PositiveRef(index_ref);
-  const int target_ref = ct.element().target();
-  const int target_var = PositiveRef(target_ref);
-
-  // Deal with non-fixed target and non-fixed index. This only happen if
-  // whatever the value of the index and selected variable, we can choose a
-  // valid target, so we just fix the index to its min value in this case.
-  if (!(*domains)[target_var].IsFixed() && !(*domains)[index_var].IsFixed()) {
-    const int64_t index_var_value = (*domains)[index_var].Min();
-    (*domains)[index_var] = Domain(index_var_value);
-
-    // If the selected variable is not fixed, we also need to fix it.
-    const int selected_ref = ct.element().vars(
-        RefIsPositive(index_ref) ? index_var_value : -index_var_value);
-    const int selected_var = PositiveRef(selected_ref);
-    if (!(*domains)[selected_var].IsFixed()) {
-      (*domains)[selected_var] = Domain((*domains)[selected_var].Min());
-    }
-  }
+  const LinearExpressionProto& index = ct.element().linear_index();
+  const LinearExpressionProto& target = ct.element().linear_target();
+
+  DCHECK(LinearExpressionIsFixed(index, *domains) ||
+         LinearExpressionIsFixed(target, *domains));
 
   // Deal with fixed index.
-  if ((*domains)[index_var].IsFixed()) {
-    const int64_t index_var_value = (*domains)[index_var].FixedValue();
-    const int selected_ref = ct.element().vars(
-        RefIsPositive(index_ref) ? index_var_value : -index_var_value);
-    const int selected_var = PositiveRef(selected_ref);
-    if ((*domains)[selected_var].IsFixed()) {
-      const int64_t selected_value = (*domains)[selected_var].FixedValue();
-      (*domains)[target_var] = (*domains)[target_var].IntersectionWith(
-          Domain(RefIsPositive(target_ref) == RefIsPositive(selected_ref)
-                     ? selected_value
-                     : -selected_value));
-      DCHECK(!(*domains)[target_var].IsEmpty());
+  if (LinearExpressionIsFixed(index, *domains)) {
+    const int64_t index_value = EvaluateLinearExpression(index, *domains);
+    const LinearExpressionProto& expr = ct.element().exprs(index_value);
+    DCHECK(LinearExpressionIsFixed(expr, *domains));
+    const int64_t expr_value = EvaluateLinearExpression(expr, *domains);
+    if (target.vars().empty()) {
+      DCHECK_EQ(expr_value, target.offset());
     } else {
-      const bool same_sign =
-          (selected_var == selected_ref) == (target_var == target_ref);
-      const Domain target_domain = (*domains)[target_var];
-      const Domain selected_domain = same_sign
-                                         ? (*domains)[selected_var]
-                                         : (*domains)[selected_var].Negation();
-      const Domain final = target_domain.IntersectionWith(selected_domain);
-      const int64_t value = final.SmallestValue();
-      (*domains)[target_var] =
-          (*domains)[target_var].IntersectionWith(Domain(value));
-      (*domains)[selected_var] = (*domains)[selected_var].IntersectionWith(
-          Domain(same_sign ? value : -value));
-      DCHECK(!(*domains)[target_var].IsEmpty());
-      DCHECK(!(*domains)[selected_var].IsEmpty());
+      (*domains)[target.vars(0)] = Domain(GetInnerVarValue(target, expr_value));
     }
     return;
   }
 
   // Deal with fixed target (and constant vars).
-  const int64_t target_value = (*domains)[target_var].FixedValue();
+  const int64_t target_value = EvaluateLinearExpression(target, *domains);
   int selected_index_value = -1;
-  for (const int64_t v : (*domains)[index_var].Values()) {
-    const int64_t i = index_var == index_ref ? v : -v;
-    if (i < 0 || i >= ct.element().vars_size()) continue;
-
-    const int ref = ct.element().vars(i);
-    const int var = PositiveRef(ref);
-    const int64_t value = (*domains)[var].FixedValue();
-    if (RefIsPositive(target_ref) == RefIsPositive(ref)) {
-      if (value == target_value) {
-        selected_index_value = i;
-        break;
-      }
-    } else {
-      if (value == -target_value) {
-        selected_index_value = i;
-        break;
-      }
+  for (const int64_t v : (*domains)[index.vars(0)].Values()) {
+    const int64_t index_value = index.offset() + v * index.coeffs(0);
+    DCHECK_GE(index_value, 0);
+    DCHECK_LT(index_value, ct.element().exprs_size());
+
+    const LinearExpressionProto& expr = ct.element().exprs(index_value);
+    const int64_t value = EvaluateLinearExpression(expr, *domains);
+    if (value == target_value) {
+      selected_index_value = index_value;
+      break;
     }
   }
 
   CHECK_NE(selected_index_value, -1);
-  (*domains)[index_var] = (*domains)[index_var].IntersectionWith(Domain(
-      RefIsPositive(index_ref) ? selected_index_value : -selected_index_value));
-  DCHECK(!(*domains)[index_var].IsEmpty());
+  (*domains)[index.vars(0)] =
+      Domain(GetInnerVarValue(index, selected_index_value));
 }
 
 // We only support assigning to an affine target.
diff --git a/ortools/sat/cp_model_postsolve_test.cc b/ortools/sat/cp_model_postsolve_test.cc
index dee1d25448..2034147100 100644
--- a/ortools/sat/cp_model_postsolve_test.cc
+++ b/ortools/sat/cp_model_postsolve_test.cc
@@ -96,7 +96,7 @@ TEST(PostsolveResponseTest, ExactlyOneExample2) {
   EXPECT_THAT(solution, ::testing::ElementsAre(0, 1, 0));
 }
 
-TEST(PostsolveResponseTest, Element) {
+TEST(PostsolveResponseTest, FixedTarget) {
   // Fixing z will allow the postsolve code to reconstruct all values.
   const CpModelProto mapping_proto = ParseTestProto(R"pb(
     variables {
@@ -105,21 +105,22 @@ TEST(PostsolveResponseTest, Element) {
     }
     variables {
       name: 'a'
-      domain: [ 1, 10 ]
+      domain: [ 5, 5 ]
     }
     variables {
       name: 'b'
-      domain: [ 0, 10 ]
+      domain: [ 7, 7 ]
     }
     variables {
       name: 'target'
-      domain: [ 0, 10 ]
+      domain: [ 7, 7 ]
     }
     constraints {
       element {
-        index: 0
-        vars: [ 1, 2 ]
-        target: 3
+        linear_index { vars: 0 coeffs: 1 }
+        linear_target { vars: 3 coeffs: 1 }
+        exprs { vars: 1, coeffs: 1 }
+        exprs { vars: 2, coeffs: 1 }
       }
     }
   )pb");
@@ -128,23 +129,31 @@ TEST(PostsolveResponseTest, Element) {
   std::vector<int> postsolve_mapping = {};
   PostsolveResponse(/*num_variables_in_original_model=*/4, mapping_proto,
                     postsolve_mapping, &solution);
-  EXPECT_THAT(solution, ::testing::ElementsAre(0, 1, 0, 1));
+  EXPECT_THAT(solution, ::testing::ElementsAre(1, 5, 7, 7));
 }
 
-TEST(PostsolveResponseTest, VariableElement) {
+TEST(PostsolveResponseTest, FixedIndex) {
   const CpModelProto mapping_proto = ParseTestProto(R"pb(
-    variables { domain: [ 0, 129 ] }
-    variables { domain: [ 1, 5 ] }
+    variables { domain: [ 7, 7 ] }
+    variables { domain: [ 3, 3 ] }
     variables { domain: [ 0, 129 ] }
     variables { domain: [ 2, 2 ] }
-    constraints { element { index: 3 target: 2 vars: 0 vars: 1 vars: 0 } }
+    constraints {
+      element {
+        linear_index { vars: 3 coeffs: 1 }
+        linear_target { vars: 2 coeffs: 1 }
+        exprs { vars: 0, coeffs: 1 }
+        exprs { vars: 1, coeffs: 1 }
+        exprs { vars: 0, coeffs: 1 }
+      }
+    }
   )pb");
 
   std::vector<int64_t> solution;
   std::vector<int> postsolve_mapping = {};
   PostsolveResponse(/*num_variables_in_original_model=*/4, mapping_proto,
                     postsolve_mapping, &solution);
-  EXPECT_THAT(solution, ::testing::ElementsAre(0, 1, 0, 2));
+  EXPECT_THAT(solution, ::testing::ElementsAre(7, 3, 7, 2));
 }
 
 // Note that our postolve code is "limited" when it come to solving a single
@@ -158,7 +167,7 @@ TEST(PostsolveResponseTest, TrickyLinearCase) {
   // there is a possible z, but the reverse is not true, since y = 1, z = 0 is
   // not feasible.
   //
-  // The preosolve should deal with that by putting z first so that the
+  // The presolve should deal with that by putting z first so that the
   // postsolve code do not fail.
   const CpModelProto mapping_proto = ParseTestProto(R"pb(
     variables {
diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc
index a283115e39..7911fd59b8 100644
--- a/ortools/sat/cp_model_presolve.cc
+++ b/ortools/sat/cp_model_presolve.cc
@@ -1143,7 +1143,7 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
   }
 
   // Checks if the affine target domain is constraining.
-  bool affine_target_domain_contains_max_domain = false;
+  bool linear_target_domain_contains_max_domain = false;
   if (ExpressionContainsSingleRef(target)) {  // target = +/- var.
     int64_t infered_min = std::numeric_limits<int64_t>::min();
     int64_t infered_max = std::numeric_limits<int64_t>::min();
@@ -1165,7 +1165,7 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
     const Domain target_domain =
         target.coeffs(0) == 1 ? context_->DomainOf(target.vars(0))
                               : context_->DomainOf(target.vars(0)).Negation();
-    affine_target_domain_contains_max_domain =
+    linear_target_domain_contains_max_domain =
         rhs_domain.IsIncludedIn(target_domain);
   }
 
@@ -1181,7 +1181,7 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
   // Actually, I think for the expr=affine' case, it reduces to:
   // affine(x) >= affine'(x)
   // affine(x) = max(...);
-  if (affine_target_domain_contains_max_domain) {
+  if (linear_target_domain_contains_max_domain) {
     const int target_var = target.vars(0);
     bool abort = false;
     for (const LinearExpressionProto& expr : ct->lin_max().exprs()) {
@@ -1199,12 +1199,12 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
       // We only know that the target is affine here.
       context_->UpdateRuleStats(
           "TODO lin_max: affine(x) = max(affine'(x), ...) !!");
-      affine_target_domain_contains_max_domain = false;
+      linear_target_domain_contains_max_domain = false;
     }
   }
 
   // If the target is not used, and safe, we can remove the constraint.
-  if (affine_target_domain_contains_max_domain &&
+  if (linear_target_domain_contains_max_domain &&
       context_->VariableIsUniqueAndRemovable(target.vars(0))) {
     context_->UpdateRuleStats("lin_max: unused affine target");
     context_->MarkVariableAsRemoved(target.vars(0));
@@ -1214,7 +1214,7 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
 
   // If the target is only used in the objective, and safe, we can simplify the
   // constraint.
-  if (affine_target_domain_contains_max_domain &&
+  if (linear_target_domain_contains_max_domain &&
       context_->VariableWithCostIsUniqueAndRemovable(target.vars(0)) &&
       (target.coeffs(0) > 0) ==
           (context_->ObjectiveCoeff(target.vars(0)) > 0)) {
@@ -4124,7 +4124,7 @@ void CpModelPresolver::LowerThanCoeffStrengthening(bool from_lower_bound,
     for (int i = 0; i < num_vars; ++i) {
       const int64_t magnitude = std::abs(arg.coeffs(i));
       if (magnitude <= second_threshold) {
-        gcd = MathUtil::GCD64(gcd, magnitude);
+        gcd = std::gcd(gcd, magnitude);
         max_magnitude_left = std::max(max_magnitude_left, magnitude);
         const int64_t bound_diff =
             context_->MaxOf(arg.vars(i)) - context_->MinOf(arg.vars(i));
@@ -4856,194 +4856,222 @@ bool CpModelPresolver::PresolveInverse(ConstraintProto* ct) {
   return false;
 }
 
-bool CpModelPresolver::PresolveElement(ConstraintProto* ct) {
+bool CpModelPresolver::PresolveElement(int c, ConstraintProto* ct) {
   if (context_->ModelIsUnsat()) return false;
 
-  if (ct->element().vars().empty()) {
+  if (ct->element().exprs().empty()) {
     context_->UpdateRuleStats("element: empty array");
     return context_->NotifyThatModelIsUnsat();
   }
 
-  const int index_ref = ct->element().index();
-  const int target_ref = ct->element().target();
+  bool changed = false;
+  changed |= CanonicalizeLinearExpression(
+      *ct, ct->mutable_element()->mutable_linear_index());
+  changed |= CanonicalizeLinearExpression(
+      *ct, ct->mutable_element()->mutable_linear_target());
+  for (int i = 0; i < ct->element().exprs_size(); ++i) {
+    changed |= CanonicalizeLinearExpression(
+        *ct, ct->mutable_element()->mutable_exprs(i));
+  }
+
+  const LinearExpressionProto& index = ct->element().linear_index();
+  const LinearExpressionProto& target = ct->element().linear_target();
 
   // TODO(user): think about this once we do have such constraint.
   if (HasEnforcementLiteral(*ct)) return false;
 
-  bool all_constants = true;
-  std::vector<int64_t> constants;
-  bool all_included_in_target_domain = true;
-
+  // Reduce index domain from the array size.
   {
+    bool index_modified = false;
     if (!context_->IntersectDomainWith(
-            index_ref, Domain(0, ct->element().vars_size() - 1))) {
+            index, Domain(0, ct->element().exprs_size() - 1),
+            &index_modified)) {
       return false;
     }
+    if (index_modified) {
+      context_->UpdateRuleStats(
+          "element: reduced index domain from array size");
+    }
+  }
 
-    // Filter impossible index values if index == +/- target
-    //
-    // Note that this must be done before the unique_index/target rule.
-    if (PositiveRef(target_ref) == PositiveRef(index_ref)) {
-      std::vector<int64_t> possible_indices;
-      const Domain& index_domain = context_->DomainOf(index_ref);
-      for (const int64_t index_value : index_domain.Values()) {
-        const int ref = ct->element().vars(index_value);
-        const int64_t target_value =
-            target_ref == index_ref ? index_value : -index_value;
-        if (context_->DomainContains(ref, target_value)) {
-          possible_indices.push_back(target_value);
-        }
-      }
-      if (possible_indices.size() < index_domain.Size()) {
-        if (!context_->IntersectDomainWith(
-                index_ref, Domain::FromValues(possible_indices))) {
-          return true;
-        }
+  // Special case if the index is fixed.
+  if (context_->IsFixed(index)) {
+    const int64_t index_value = context_->FixedValue(index);
+    ConstraintProto* new_ct = context_->working_model->add_constraints();
+    new_ct->mutable_linear()->add_domain(0);
+    new_ct->mutable_linear()->add_domain(0);
+    AddLinearExpressionToLinearConstraint(target, 1, new_ct->mutable_linear());
+    AddLinearExpressionToLinearConstraint(ct->element().exprs(index_value), -1,
+                                          new_ct->mutable_linear());
+    context_->CanonicalizeLinearConstraint(new_ct);
+    context_->UpdateNewConstraintsVariableUsage();
+    context_->UpdateRuleStats("element: fixed index");
+    return RemoveConstraint(ct);
+  }
+
+  // We know index is not fixed.
+  const int index_var = index.vars(0);
+
+  {
+    // Cleanup the array: if exprs[i] contains index_var, fix its value.
+    const Domain& index_var_domain = context_->DomainOf(index_var);
+    std::vector<int64_t> reached_indices(ct->element().exprs_size(), false);
+    for (const int64_t index_var_value : index_var_domain.Values()) {
+      const int64_t index_value =
+          AffineExpressionValueAt(index, index_var_value);
+      reached_indices[index_value] = true;
+      const LinearExpressionProto& expr = ct->element().exprs(index_value);
+      if (expr.vars_size() == 1 && expr.vars(0) == index_var) {
+        const int64_t expr_value =
+            AffineExpressionValueAt(expr, index_var_value);
+        ct->mutable_element()->mutable_exprs(index_value)->clear_vars();
+        ct->mutable_element()->mutable_exprs(index_value)->clear_coeffs();
+        ct->mutable_element()
+            ->mutable_exprs(index_value)
+            ->set_offset(expr_value);
+        changed = true;
         context_->UpdateRuleStats(
-            "element: reduced index domain when target equals index");
+            "element: fix expression depending on the index");
       }
     }
 
-    // Filter possible index values. Accumulate variable domains to build
-    // a possible target domain.
-    Domain infered_domain;
-    const Domain& initial_index_domain = context_->DomainOf(index_ref);
-    const Domain& target_domain = context_->DomainOf(target_ref);
-    std::vector<int64_t> possible_indices;
-    for (const int64_t value : initial_index_domain.Values()) {
-      CHECK_GE(value, 0);
-      CHECK_LT(value, ct->element().vars_size());
-      const int ref = ct->element().vars(value);
-
-      // We cover the corner cases where the possible domain is actually fixed.
-      Domain domain = context_->DomainOf(ref);
-      if (ref == index_ref) {
-        domain = Domain(value);
-      } else if (ref == NegatedRef(index_ref)) {
-        domain = Domain(-value);
-      } else if (ref == NegatedRef(target_ref)) {
-        domain = Domain(0);
-      }
-
-      if (domain.IntersectionWith(target_domain).IsEmpty()) continue;
-      possible_indices.push_back(value);
-      if (domain.IsFixed()) {
-        constants.push_back(domain.Min());
-      } else {
-        all_constants = false;
+    // Cleanup the array: clear unreached expressions.
+    for (int i = 0; i < ct->element().exprs_size(); ++i) {
+      if (!reached_indices[i]) {
+        ct->mutable_element()->mutable_exprs(i)->Clear();
+        changed = true;
       }
-      if (!domain.IsIncludedIn(target_domain)) {
-        all_included_in_target_domain = false;
+    }
+  }
+
+  // Canonicalization and cleanups of the expressions could have messed up the
+  // var-constraint graph.
+  if (changed) context_->UpdateConstraintVariableUsage(c);
+
+  // Reduces the domain of the index.
+  {
+    const Domain& index_var_domain = context_->DomainOf(index_var);
+    const Domain& target_domain = context_->DomainSuperSetOf(target);
+    std::vector<int64_t> possible_index_var_values;
+    for (const int64_t index_var_value : index_var_domain.Values()) {
+      const int64_t index_value =
+          AffineExpressionValueAt(index, index_var_value);
+      const LinearExpressionProto& expr = ct->element().exprs(index_value);
+
+      // The target domain can be reduced if it shares its variable with the
+      // index.
+      Domain reduced_target_domain = target_domain;
+      if (target.vars_size() == 1 && target.vars(0) == index_var) {
+        reduced_target_domain =
+            Domain(AffineExpressionValueAt(target, index_var_value));
+      }
+
+      // TODO(user): Implement a more precise test here.
+      if (reduced_target_domain
+              .IntersectionWith(context_->DomainSuperSetOf(expr))
+              .IsEmpty()) {
+        ct->mutable_element()->mutable_exprs(index_value)->Clear();
+        changed = true;
+      } else {
+        possible_index_var_values.push_back(index_var_value);
       }
-      infered_domain = infered_domain.UnionWith(domain);
     }
-    if (possible_indices.size() < initial_index_domain.Size()) {
+    if (possible_index_var_values.size() < index_var_domain.Size()) {
       if (!context_->IntersectDomainWith(
-              index_ref, Domain::FromValues(possible_indices))) {
+              index_var, Domain::FromValues(possible_index_var_values))) {
         return true;
       }
-      context_->UpdateRuleStats("element: reduced index domain");
+      context_->UpdateRuleStats("element: reduced index domain ");
+      // If the index is fixed, this is a equality constraint.
+      if (context_->IsFixed(index)) {
+        ConstraintProto* const eq = context_->working_model->add_constraints();
+        eq->mutable_linear()->add_domain(0);
+        eq->mutable_linear()->add_domain(0);
+        AddLinearExpressionToLinearConstraint(target, 1, eq->mutable_linear());
+        AddLinearExpressionToLinearConstraint(
+            ct->element().exprs(context_->FixedValue(index)), -1,
+            eq->mutable_linear());
+        context_->CanonicalizeLinearConstraint(eq);
+        context_->UpdateNewConstraintsVariableUsage();
+        context_->UpdateRuleStats("element: fixed index");
+        return RemoveConstraint(ct);
+      }
+    }
+  }
+
+  bool all_included_in_target_domain = true;
+  {
+    // Accumulate expressions domains to build a superset of the target domain.
+    Domain infered_domain;
+    const Domain& index_var_domain = context_->DomainOf(index_var);
+    const Domain& target_domain = context_->DomainSuperSetOf(target);
+    for (const int64_t index_var_value : index_var_domain.Values()) {
+      const int64_t index_value =
+          AffineExpressionValueAt(index, index_var_value);
+      CHECK_GE(index_value, 0);
+      CHECK_LT(index_value, ct->element().exprs_size());
+      const LinearExpressionProto& expr = ct->element().exprs(index_value);
+      const Domain expr_domain = context_->DomainSuperSetOf(expr);
+      if (!expr_domain.IsIncludedIn(target_domain)) {
+        all_included_in_target_domain = false;
+      }
+      infered_domain = infered_domain.UnionWith(expr_domain);
     }
+
     bool domain_modified = false;
-    if (!context_->IntersectDomainWith(target_ref, infered_domain,
+    if (!context_->IntersectDomainWith(target, infered_domain,
                                        &domain_modified)) {
       return true;
     }
     if (domain_modified) {
-      context_->UpdateRuleStats("element: reduced target domain");
+      context_->UpdateRuleStats("element: reduce target domain");
     }
   }
 
-  // If the index is fixed, this is a equality constraint.
-  if (context_->IsFixed(index_ref)) {
-    const int var = ct->element().vars(context_->MinOf(index_ref));
-    if (var != target_ref) {
-      LinearConstraintProto* const lin =
-          context_->working_model->add_constraints()->mutable_linear();
-      lin->add_vars(var);
-      lin->add_coeffs(-1);
-      lin->add_vars(target_ref);
-      lin->add_coeffs(1);
-      lin->add_domain(0);
-      lin->add_domain(0);
-      context_->UpdateNewConstraintsVariableUsage();
+  bool all_constants = true;
+  {
+    const Domain& index_var_domain = context_->DomainOf(index_var);
+    std::vector<int64_t> expr_constants;
+
+    for (const int64_t index_var_value : index_var_domain.Values()) {
+      const int64_t index_value =
+          AffineExpressionValueAt(index, index_var_value);
+      const LinearExpressionProto& expr = ct->element().exprs(index_value);
+      if (context_->IsFixed(expr)) {
+        expr_constants.push_back(context_->FixedValue(expr));
+      } else {
+        all_constants = false;
+        break;
+      }
     }
-    context_->UpdateRuleStats("element: fixed index");
-    return RemoveConstraint(ct);
   }
 
   // If the accessible part of the array is made of a single constant value,
   // then we do not care about the index. And, because of the previous target
   // domain reduction, the target is also fixed.
-  if (all_constants && context_->IsFixed(target_ref)) {
+  if (all_constants && context_->IsFixed(target)) {
     context_->UpdateRuleStats("element: one value array");
     return RemoveConstraint(ct);
   }
 
   // Special case when the index is boolean, and the array does not contain
   // variables.
-  if (context_->MinOf(index_ref) == 0 && context_->MaxOf(index_ref) == 1 &&
+  if (context_->MinOf(index) == 0 && context_->MaxOf(index) == 1 &&
       all_constants) {
-    const int64_t v0 = constants[0];
-    const int64_t v1 = constants[1];
-
-    LinearConstraintProto* const lin =
-        context_->working_model->add_constraints()->mutable_linear();
-    lin->add_vars(target_ref);
-    lin->add_coeffs(1);
-    lin->add_vars(index_ref);
-    lin->add_coeffs(v0 - v1);
-    lin->add_domain(v0);
-    lin->add_domain(v0);
+    const int64_t v0 = context_->FixedValue(ct->element().exprs(0));
+    const int64_t v1 = context_->FixedValue(ct->element().exprs(1));
+
+    ConstraintProto* const eq = context_->working_model->add_constraints();
+    eq->mutable_linear()->add_domain(v0);
+    eq->mutable_linear()->add_domain(v0);
+    AddLinearExpressionToLinearConstraint(target, 1, eq->mutable_linear());
+    AddLinearExpressionToLinearConstraint(index, v0 - v1, eq->mutable_linear());
+    context_->CanonicalizeLinearConstraint(eq);
     context_->UpdateNewConstraintsVariableUsage();
     context_->UpdateRuleStats("element: linearize constant element of size 2");
     return RemoveConstraint(ct);
   }
 
-  // If the index has a canonical affine representative, rewrite the element.
-  const AffineRelation::Relation r_index =
-      context_->GetAffineRelation(index_ref);
-  if (r_index.representative != index_ref) {
-    // Checks the domains are synchronized.
-    if (context_->DomainOf(r_index.representative).Size() >
-        context_->DomainOf(index_ref).Size()) {
-      // Postpone, we will come back later when domains are synchronized.
-      return true;
-    }
-    const int r_ref = r_index.representative;
-    const int64_t r_min = context_->MinOf(r_ref);
-    const int64_t r_max = context_->MaxOf(r_ref);
-    const int array_size = ct->element().vars_size();
-    if (r_min != 0) {
-      context_->UpdateRuleStats("TODO element: representative has bad domain");
-    } else if (r_index.offset >= 0 && r_index.offset < array_size &&
-               r_index.offset + r_max * r_index.coeff >= 0 &&
-               r_index.offset + r_max * r_index.coeff < array_size) {
-      // This will happen eventually when domains are synchronized.
-      ElementConstraintProto* const element =
-          context_->working_model->add_constraints()->mutable_element();
-      for (int64_t v = 0; v <= r_max; ++v) {
-        const int64_t scaled_index = v * r_index.coeff + r_index.offset;
-        CHECK_GE(scaled_index, 0);
-        CHECK_LT(scaled_index, array_size);
-        element->add_vars(ct->element().vars(scaled_index));
-      }
-      element->set_index(r_ref);
-      element->set_target(target_ref);
-
-      if (r_index.coeff == 1) {
-        context_->UpdateRuleStats("element: shifed index ");
-      } else {
-        context_->UpdateRuleStats("element: scaled index");
-      }
-      context_->UpdateNewConstraintsVariableUsage();
-      return RemoveConstraint(ct);
-    }
-  }
-
-  // Should have been taken care of earlier.
-  DCHECK(!context_->IsFixed(index_ref));
-
   // If a variable (target or index) appears only in this constraint, it does
   // not necessarily mean that we can remove the constraint, as the variable
   // can be used multiple times in the element. So let's count the local
@@ -5052,26 +5080,31 @@ bool CpModelPresolver::PresolveElement(ConstraintProto* ct) {
   // TODO(user): now that we used fixed values for these case, this is no longer
   // needed I think.
   absl::flat_hash_map<int, int> local_var_occurrence_counter;
-  local_var_occurrence_counter[PositiveRef(index_ref)]++;
-  local_var_occurrence_counter[PositiveRef(target_ref)]++;
-
-  for (const ClosedInterval interval : context_->DomainOf(index_ref)) {
-    for (int64_t value = interval.start; value <= interval.end; ++value) {
-      DCHECK_GE(value, 0);
-      DCHECK_LT(value, ct->element().vars_size());
-      const int ref = ct->element().vars(value);
-      local_var_occurrence_counter[PositiveRef(ref)]++;
+  {
+    auto count = [&local_var_occurrence_counter](
+                     const LinearExpressionProto& expr) mutable {
+      for (const int var : expr.vars()) {
+        local_var_occurrence_counter[var]++;
+      }
+    };
+    count(index);
+    count(target);
+    for (const int64_t index_var_value :
+         context_->DomainOf(index_var).Values()) {
+      count(
+          ct->element().exprs(AffineExpressionValueAt(index, index_var_value)));
     }
   }
 
-  if (context_->VariableIsUniqueAndRemovable(index_ref) &&
-      local_var_occurrence_counter.at(PositiveRef(index_ref)) == 1) {
+  if (context_->VariableIsUniqueAndRemovable(index_var) &&
+      local_var_occurrence_counter.at(index_var) == 1) {
     if (all_constants) {
       // This constraint is just here to reduce the domain of the target! We can
       // add it to the mapping_model to reconstruct the index value during
       // postsolve and get rid of it now.
-      context_->UpdateRuleStats("element: trivial target domain reduction");
-      context_->MarkVariableAsRemoved(index_ref);
+      context_->UpdateRuleStats(
+          "element: removed  as the index is not used elsewhere");
+      context_->MarkVariableAsRemoved(index_var);
       context_->NewMappingConstraint(*ct, __FILE__, __LINE__);
       return RemoveConstraint(ct);
     } else {
@@ -5079,12 +5112,13 @@ bool CpModelPresolver::PresolveElement(ConstraintProto* ct) {
     }
   }
 
-  if (!context_->IsFixed(target_ref) &&
-      context_->VariableIsUniqueAndRemovable(target_ref) &&
-      local_var_occurrence_counter.at(PositiveRef(target_ref)) == 1) {
-    if (all_included_in_target_domain) {
-      context_->UpdateRuleStats("element: trivial index domain reduction");
-      context_->MarkVariableAsRemoved(target_ref);
+  if (target.vars_size() == 1 && !context_->IsFixed(target.vars(0)) &&
+      context_->VariableIsUniqueAndRemovable(target.vars(0)) &&
+      local_var_occurrence_counter.at(target.vars(0)) == 1) {
+    if (all_included_in_target_domain && std::abs(target.coeffs(0)) == 1) {
+      context_->UpdateRuleStats(
+          "element: removed as the target is not used elsewhere");
+      context_->MarkVariableAsRemoved(target.vars(0));
       context_->NewMappingConstraint(*ct, __FILE__, __LINE__);
       return RemoveConstraint(ct);
     } else {
@@ -5092,7 +5126,7 @@ bool CpModelPresolver::PresolveElement(ConstraintProto* ct) {
     }
   }
 
-  return false;
+  return changed;
 }
 
 bool CpModelPresolver::PresolveTable(ConstraintProto* ct) {
@@ -8271,7 +8305,7 @@ bool CpModelPresolver::PresolveOneConstraint(int c) {
     case ConstraintProto::kInverse:
       return PresolveInverse(ct);
     case ConstraintProto::kElement:
-      return PresolveElement(ct);
+      return PresolveElement(c, ct);
     case ConstraintProto::kTable:
       return PresolveTable(ct);
     case ConstraintProto::kAllDiff:
@@ -12013,6 +12047,9 @@ bool ModelCopy::ImportAndSimplifyConstraints(
       case ConstraintProto::kIntProd:
         if (!CopyIntProd(ct, ignore_names)) return CreateUnsatModel(c, ct);
         break;
+      case ConstraintProto::kElement:
+        if (!CopyElement(ct)) return CreateUnsatModel(c, ct);
+        break;
       case ConstraintProto::kAtMostOne:
         if (!CopyAtMostOne(ct)) return CreateUnsatModel(c, ct);
         break;
@@ -12363,6 +12400,8 @@ bool ModelCopy::CopyLinear(const ConstraintProto& ct) {
     return !temp_literals_.empty();
   }
 
+  DCHECK(!non_fixed_variables_.empty());
+
   ConstraintProto* new_ct = context_->working_model->add_constraints();
   FinishEnforcementCopy(new_ct);
   LinearConstraintProto* linear = new_ct->mutable_linear();
@@ -12374,6 +12413,36 @@ bool ModelCopy::CopyLinear(const ConstraintProto& ct) {
   return true;
 }
 
+bool ModelCopy::CopyElement(const ConstraintProto& ct) {
+  ConstraintProto* new_ct = context_->working_model->add_constraints();
+  if (ct.element().vars().empty() && !ct.element().exprs().empty()) {
+    // New format, just copy.
+    *new_ct = ct;
+    return true;
+  }
+
+  auto fill_expr = [this](int var, LinearExpressionProto* expr) mutable {
+    if (context_->IsFixed(var)) {
+      expr->set_offset(context_->FixedValue(var));
+    } else {
+      DCHECK(RefIsPositive(var));
+      expr->mutable_vars()->Reserve(1);
+      expr->mutable_coeffs()->Reserve(1);
+      expr->add_vars(var);
+      expr->add_coeffs(1);
+    }
+  };
+
+  fill_expr(ct.element().index(),
+            new_ct->mutable_element()->mutable_linear_index());
+  fill_expr(ct.element().target(),
+            new_ct->mutable_element()->mutable_linear_target());
+  for (const int var : ct.element().vars()) {
+    fill_expr(var, new_ct->mutable_element()->add_exprs());
+  }
+  return true;
+}
+
 bool ModelCopy::CopyAtMostOne(const ConstraintProto& ct) {
   int num_true = 0;
   temp_literals_.clear();
diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h
index 6759df2d0e..47926f7acf 100644
--- a/ortools/sat/cp_model_presolve.h
+++ b/ortools/sat/cp_model_presolve.h
@@ -122,7 +122,7 @@ class CpModelPresolver {
   // TODO(user): Make these public and unit test.
   bool PresolveAllDiff(ConstraintProto* ct);
   bool PresolveAutomaton(ConstraintProto* ct);
-  bool PresolveElement(ConstraintProto* ct);
+  bool PresolveElement(int c, ConstraintProto* ct);
   bool PresolveIntDiv(int c, ConstraintProto* ct);
   bool PresolveIntMod(int c, ConstraintProto* ct);
   bool PresolveIntProd(ConstraintProto* ct);
@@ -432,6 +432,7 @@ class ModelCopy {
                             LinearExpressionProto* dst);
   bool CopyIntProd(const ConstraintProto& ct, bool ignore_names);
   bool CopyLinear(const ConstraintProto& ct);
+  bool CopyElement(const ConstraintProto& ct);
   bool CopyAtMostOne(const ConstraintProto& ct);
   bool CopyExactlyOne(const ConstraintProto& ct);
 
diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc
index 9420aa4c6d..de7b61ef84 100644
--- a/ortools/sat/cp_model_solver.cc
+++ b/ortools/sat/cp_model_solver.cc
@@ -1628,7 +1628,7 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) {
         name_filter.Keep("lb_relax_lns")) {
       reentrant_interleaved_subsolvers.push_back(std::make_unique<LnsSolver>(
           std::make_unique<LocalBranchingLpBasedNeighborhoodGenerator>(
-              helper, name_filter.LastName(), shared->time_limit),
+              helper, name_filter.LastName(), shared->time_limit, shared),
           lns_params, helper, shared));
     }
 
@@ -1718,12 +1718,6 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) {
                 helper, name_filter.LastName()),
             lns_params, helper, shared));
       }
-      if (name_filter.Keep("routing_starts_lns")) {
-        reentrant_interleaved_subsolvers.push_back(std::make_unique<LnsSolver>(
-            std::make_unique<RoutingStartsNeighborhoodGenerator>(
-                helper, name_filter.LastName()),
-            lns_params, helper, shared));
-      }
     }
     if (num_routes > 0 || num_circuit > 1) {
       if (name_filter.Keep("routing_full_path_lns")) {
@@ -2399,6 +2393,16 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
   SOLVER_LOG(logger, "");
   SOLVER_LOG(logger, "Presolved ", CpModelStats(*new_cp_model_proto));
 
+  // Detect the symmetry of the presolved model.
+  // Note that this needs to be done before SharedClasses are created.
+  //
+  // TODO(user): We could actually report a complete feasible hint before this
+  // point. But the proper fix is to report it even before the presolve.
+  if (params.symmetry_level() > 1 && !params.stop_after_presolve() &&
+      !shared_time_limit->LimitReached()) {
+    DetectAndAddSymmetryToProto(params, new_cp_model_proto, logger);
+  }
+
   // TODO(user): reduce this function size and find a better place for this?
   SharedClasses shared(new_cp_model_proto, model);
 
@@ -2584,10 +2588,6 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
     TestSolutionHintForFeasibility(*new_cp_model_proto, logger, nullptr);
   }
 
-  if (params.symmetry_level() > 1) {
-    DetectAndAddSymmetryToProto(params, new_cp_model_proto, logger);
-  }
-
   LoadDebugSolution(*new_cp_model_proto, model);
 
   if (!model->GetOrCreate<TimeLimit>()->LimitReached()) {
diff --git a/ortools/sat/cp_model_solver_helpers.cc b/ortools/sat/cp_model_solver_helpers.cc
index b961bcf0f7..a7e50758d2 100644
--- a/ortools/sat/cp_model_solver_helpers.cc
+++ b/ortools/sat/cp_model_solver_helpers.cc
@@ -517,8 +517,11 @@ IntegerVariable AddLPConstraints(bool objective_need_to_be_tight,
       m->TakeOwnership(lp_constraints[c]);
     }
     // Load the constraint.
-    lp_constraints[c]->AddLinearConstraint(
-        std::move(relaxation.linear_constraints[i]));
+    if (!lp_constraints[c]->AddLinearConstraint(
+            std::move(relaxation.linear_constraints[i]))) {
+      m->GetOrCreate<SatSolver>()->NotifyThatModelIsUnsat();
+      return kNoIntegerVariable;
+    }
   }
 
   // Dispatch every cut generator to its LinearProgrammingConstraint.
@@ -532,6 +535,16 @@ IntegerVariable AddLPConstraints(bool objective_need_to_be_tight,
     lp_constraints[c]->AddCutGenerator(std::move(relaxation.cut_generators[i]));
   }
 
+  // We deal with the clique cut generator here now that the component have
+  // been computed. As we don't want to merge independent component with it.
+  if (params->linearization_level() > 1 && params->add_clique_cuts() &&
+      params->cut_level() > 0) {
+    for (LinearProgrammingConstraint* lp : lp_constraints) {
+      if (lp == nullptr) continue;
+      lp->AddCutGenerator(CreateCliqueCutGenerator(lp->integer_variables(), m));
+    }
+  }
+
   // Add the objective.
   std::vector<std::vector<std::pair<IntegerVariable, int64_t>>>
       component_to_cp_terms(num_components);
diff --git a/ortools/sat/cp_model_solver_helpers.h b/ortools/sat/cp_model_solver_helpers.h
index a2220f3b95..915bd91c17 100644
--- a/ortools/sat/cp_model_solver_helpers.h
+++ b/ortools/sat/cp_model_solver_helpers.h
@@ -17,6 +17,7 @@
 #include <cstdint>
 #include <memory>
 #include <string>
+#include <tuple>
 #include <vector>
 
 #include "absl/flags/declare.h"
diff --git a/ortools/sat/cp_model_utils.cc b/ortools/sat/cp_model_utils.cc
index 2af03ee0e2..6cd32ca6f7 100644
--- a/ortools/sat/cp_model_utils.cc
+++ b/ortools/sat/cp_model_utils.cc
@@ -142,9 +142,20 @@ void GetReferencesUsedByConstraint(const ConstraintProto& ct,
       AddIndices(ct.dummy_constraint().vars(), variables);
       break;
     case ConstraintProto::ConstraintCase::kElement:
-      variables->push_back(ct.element().index());
-      variables->push_back(ct.element().target());
-      AddIndices(ct.element().vars(), variables);
+      if (ct.element().index() != 0 || ct.element().target() != 0 ||
+          !ct.element().vars().empty()) {
+        variables->push_back(ct.element().index());
+        variables->push_back(ct.element().target());
+        AddIndices(ct.element().vars(), variables);
+      } else if (ct.element().has_linear_index() ||
+                 ct.element().has_linear_target() ||
+                 !ct.element().exprs().empty()) {
+        AddIndices(ct.element().linear_index().vars(), variables);
+        AddIndices(ct.element().linear_target().vars(), variables);
+        for (const LinearExpressionProto& expr : ct.element().exprs()) {
+          AddIndices(expr.vars(), variables);
+        }
+      }
       break;
     case ConstraintProto::ConstraintCase::kCircuit:
       AddIndices(ct.circuit().literals(), literals);
@@ -316,9 +327,20 @@ void ApplyToAllVariableIndices(const std::function<void(int*)>& f,
       APPLY_TO_REPEATED_FIELD(dummy_constraint, vars);
       break;
     case ConstraintProto::ConstraintCase::kElement:
-      APPLY_TO_SINGULAR_FIELD(element, index);
-      APPLY_TO_SINGULAR_FIELD(element, target);
-      APPLY_TO_REPEATED_FIELD(element, vars);
+      if (ct->element().index() != 0 || ct->element().target() != 0 ||
+          !ct->element().vars().empty()) {
+        APPLY_TO_SINGULAR_FIELD(element, index);
+        APPLY_TO_SINGULAR_FIELD(element, target);
+        APPLY_TO_REPEATED_FIELD(element, vars);
+      } else if (ct->element().has_linear_index() ||
+                 ct->element().has_linear_target() ||
+                 !ct->element().exprs().empty()) {
+        APPLY_TO_REPEATED_FIELD(element, linear_index()->mutable_vars);
+        APPLY_TO_REPEATED_FIELD(element, linear_target()->mutable_vars);
+        for (int i = 0; i < ct->element().exprs_size(); ++i) {
+          APPLY_TO_REPEATED_FIELD(element, exprs(i)->mutable_vars);
+        }
+      }
       break;
     case ConstraintProto::ConstraintCase::kCircuit:
       break;
@@ -710,6 +732,11 @@ uint64_t FingerprintModel(const CpModelProto& model, uint64_t seed) {
         fp = FingerprintSingleField(ct.element().index(), fp);
         fp = FingerprintSingleField(ct.element().target(), fp);
         fp = FingerprintRepeatedField(ct.element().vars(), fp);
+        fp = FingerprintExpression(ct.element().linear_index(), fp);
+        fp = FingerprintExpression(ct.element().linear_target(), fp);
+        for (const LinearExpressionProto& expr : ct.element().exprs()) {
+          fp = FingerprintExpression(expr, fp);
+        }
         break;
       case ConstraintProto::ConstraintCase::kCircuit:
         fp = FingerprintRepeatedField(ct.circuit().heads(), fp);
diff --git a/ortools/sat/cp_model_utils.h b/ortools/sat/cp_model_utils.h
index 2bb793e2ac..ad2d652ffb 100644
--- a/ortools/sat/cp_model_utils.h
+++ b/ortools/sat/cp_model_utils.h
@@ -205,6 +205,29 @@ bool ExpressionIsAffine(const LinearExpressionProto& expr);
 // ExpressionContainsSingleRef(expr) is true.
 int GetSingleRefFromExpression(const LinearExpressionProto& expr);
 
+// Evaluates an affine expression at the given value.
+inline int64_t AffineExpressionValueAt(const LinearExpressionProto& expr,
+                                       int64_t value) {
+  CHECK_EQ(1, expr.vars_size());
+  return expr.offset() + value * expr.coeffs(0);
+}
+
+// Returns the value of the inner variable of an affine expression from the
+// value of the expression. It will DCHECK that the result is valid.
+inline int64_t GetInnerVarValue(const LinearExpressionProto& expr,
+                                int64_t value) {
+  DCHECK_EQ(expr.vars_size(), 1);
+  const int64_t var_value = (value - expr.offset()) / expr.coeffs(0);
+  DCHECK_EQ(value, var_value * expr.coeffs(0) + expr.offset());
+  return var_value;
+}
+
+// Returns true if the expression is a * var + b.
+inline bool AffineExpressionContainsVar(const LinearExpressionProto& expr,
+                                        int var) {
+  return expr.vars_size() == 1 && expr.vars(0) == var;
+}
+
 // Adds a linear expression proto to a linear constraint in place.
 //
 // Important: The domain must already be set, otherwise the offset will be lost.
diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc
index f9025e06dc..296ff8e21c 100644
--- a/ortools/sat/linear_constraint_manager.cc
+++ b/ortools/sat/linear_constraint_manager.cc
@@ -101,7 +101,8 @@ void LinearConstraintSymmetrizer::AddSymmetryOrbit(
 //
 // We will remap & scale the constraint.
 // If not possible, we will drop it for now.
-bool LinearConstraintSymmetrizer::FoldLinearConstraint(LinearConstraint* ct) {
+bool LinearConstraintSymmetrizer::FoldLinearConstraint(LinearConstraint* ct,
+                                                       bool* folded) {
   if (!has_symmetry_) return true;
 
   // We assume the constraint had basic preprocessing with tight lb/ub for
@@ -133,6 +134,8 @@ bool LinearConstraintSymmetrizer::FoldLinearConstraint(LinearConstraint* ct) {
     return true;
   }
 
+  if (folded != nullptr) *folded = true;
+
   // We need to multiply each term by scaling_factor / orbit_size.
   //
   // TODO(user): Now that we know the actual coefficient we could scale less.
@@ -302,7 +305,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
 // to detect duplicate constraints and merge bounds. This is also relevant if
 // we regenerate identical cuts for some reason.
 LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
-    LinearConstraint ct, bool* added) {
+    LinearConstraint ct, bool* added, bool* folded) {
   DCHECK_GT(ct.num_terms, 0);
   DCHECK(!PossibleOverflow(integer_trail_, ct)) << ct.DebugString();
   DCHECK(NoDuplicateVariable(ct));
@@ -313,7 +316,8 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
 
   // If configured, store instead the folded version of this constraint.
   // TODO(user): Shall we simplify again?
-  if (symmetrizer_->HasSymmetry() && !symmetrizer_->FoldLinearConstraint(&ct)) {
+  if (symmetrizer_->HasSymmetry() &&
+      !symmetrizer_->FoldLinearConstraint(&ct, folded)) {
     return kInvalidConstraintIndex;
   }
   CHECK(std::is_sorted(ct.VarsAsSpan().begin(), ct.VarsAsSpan().end()));
diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h
index a58d09dc59..746064ddc3 100644
--- a/ortools/sat/linear_constraint_manager.h
+++ b/ortools/sat/linear_constraint_manager.h
@@ -66,11 +66,14 @@ struct ModelReducedCosts
 // constraint and summing them, we can derive a new constraint using the sum
 // of the variable on each orbit instead of the individual variables.
 //
-// For the integration in a MIP solver, I couldn't find any reference. The way
+// For the integration in a MIP solver, I couldn't find many reference. The way
 // I did it here is to introduce for each orbit a variable representing the
 // sum of the orbit variable. This allows to represent the folded LP in terms
 // of these variables (that are connected to the rest of the solver) and just
 // reuse the full machinery.
+//
+// There are related info in "Orbital Shrinking", Matteo Fischetti & Leo
+// Liberti, https://link.springer.com/chapter/10.1007/978-3-642-32147-4_6
 class LinearConstraintSymmetrizer {
  public:
   explicit LinearConstraintSymmetrizer(Model* model)
@@ -113,7 +116,7 @@ class LinearConstraintSymmetrizer {
   //
   // Preconditions: All IntegerVariable must be positive. And the constraint
   // lb/ub must be tight and not +/- int64_t max.
-  bool FoldLinearConstraint(LinearConstraint* ct);
+  bool FoldLinearConstraint(LinearConstraint* ct, bool* folded = nullptr);
 
  private:
   SharedStatistics* shared_stats_;
@@ -192,7 +195,8 @@ class LinearConstraintManager {
   // constraint was actually a new one and to false if it was dominated by an
   // already existing one.
   DEFINE_STRONG_INDEX_TYPE(ConstraintIndex);
-  ConstraintIndex Add(LinearConstraint ct, bool* added = nullptr);
+  ConstraintIndex Add(LinearConstraint ct, bool* added = nullptr,
+                      bool* folded = nullptr);
 
   // Same as Add(), but logs some information about the newly added constraint.
   // Cuts are also handled slightly differently than normal constraints.
diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc
index c2382745cb..eb7771c504 100644
--- a/ortools/sat/linear_programming_constraint.cc
+++ b/ortools/sat/linear_programming_constraint.cc
@@ -48,7 +48,6 @@
 #include "ortools/lp_data/lp_types.h"
 #include "ortools/lp_data/scattered_vector.h"
 #include "ortools/lp_data/sparse.h"
-#include "ortools/lp_data/sparse_column.h"
 #include "ortools/sat/cp_model_mapping.h"
 #include "ortools/sat/cuts.h"
 #include "ortools/sat/implied_bounds.h"
@@ -56,6 +55,7 @@
 #include "ortools/sat/integer_expr.h"
 #include "ortools/sat/linear_constraint.h"
 #include "ortools/sat/linear_constraint_manager.h"
+#include "ortools/sat/linear_propagation.h"
 #include "ortools/sat/model.h"
 #include "ortools/sat/sat_base.h"
 #include "ortools/sat/sat_parameters.pb.h"
@@ -282,6 +282,7 @@ LinearProgrammingConstraint::LinearProgrammingConstraint(
       shared_response_manager_(model->GetOrCreate<SharedResponseManager>()),
       random_(model->GetOrCreate<ModelRandomGenerator>()),
       symmetrizer_(model->GetOrCreate<LinearConstraintSymmetrizer>()),
+      linear_propagator_(model->GetOrCreate<LinearPropagator>()),
       rlt_cut_helper_(model),
       implied_bounds_processor_({}, integer_trail_,
                                 model->GetOrCreate<ImpliedBounds>()),
@@ -317,18 +318,36 @@ LinearProgrammingConstraint::LinearProgrammingConstraint(
   // presence of symmetry. However they can still appear in cut, so it is a
   // bit tricky and require some refactoring to be tried.
   ColIndex col{0};
-  integer_variables_.assign(vars.begin(), vars.end());
   for (const IntegerVariable positive_variable : vars) {
     CHECK(VariableIsPositive(positive_variable));
     implied_bounds_processor_.AddLpVariable(positive_variable);
     (*dispatcher_)[positive_variable] = this;
+
+    if (!symmetrizer_->AppearInFoldedProblem(positive_variable)) continue;
+
+    integer_variables_.push_back(positive_variable);
+    extended_integer_variables_.push_back(positive_variable);
     mirror_lp_variable_[positive_variable] = col;
     ++col;
   }
 
-  lp_solution_.assign(integer_variables_.size(),
-                      std::numeric_limits<double>::infinity());
-  lp_reduced_cost_.assign(integer_variables_.size(), 0.0);
+  // Complete the extended variables with the orbit afterwards.
+  if (symmetrizer_->HasSymmetry()) {
+    for (const IntegerVariable var : integer_variables_) {
+      if (!symmetrizer_->IsOrbitSumVar(var)) continue;
+      const int orbit_index = symmetrizer_->OrbitIndex(var);
+      for (const IntegerVariable var : symmetrizer_->Orbit(orbit_index)) {
+        extended_integer_variables_.push_back(var);
+        mirror_lp_variable_[var] = col;
+        ++col;
+      }
+    }
+  }
+
+  // We use the "extended" vector for lp_solution_/lp_reduced_cost_.
+  CHECK_EQ(vars.size(), extended_integer_variables_.size());
+  lp_solution_.assign(vars.size(), std::numeric_limits<double>::infinity());
+  lp_reduced_cost_.assign(vars.size(), 0.0);
 
   if (!vars.empty()) {
     const int max_index = NegationOf(vars.back()).value();
@@ -341,9 +360,34 @@ LinearProgrammingConstraint::LinearProgrammingConstraint(
   }
 }
 
-void LinearProgrammingConstraint::AddLinearConstraint(LinearConstraint ct) {
+bool LinearProgrammingConstraint::AddLinearConstraint(LinearConstraint ct) {
   DCHECK(!lp_constraint_is_registered_);
-  constraint_manager_.Add(std::move(ct));
+  bool added = false;
+  bool folded = false;
+  const auto index = constraint_manager_.Add(std::move(ct), &added, &folded);
+
+  // If we added a folded constraints, lets add it to the CP propagators.
+  // This is important as this should tighten the bounds of the orbit sum vars.
+  if (added && folded) {
+    const auto& info = constraint_manager_.AllConstraints()[index];
+    const LinearConstraint& new_ct = info.constraint;
+    const absl::Span<const IntegerVariable> vars = new_ct.VarsAsSpan();
+    if (!info.ub_is_trivial) {
+      if (!linear_propagator_->AddConstraint({}, vars, new_ct.CoeffsAsSpan(),
+                                             new_ct.ub)) {
+        return false;
+      }
+    }
+    if (!info.lb_is_trivial) {
+      tmp_vars_.assign(vars.begin(), vars.end());
+      for (IntegerVariable& var : tmp_vars_) var = NegationOf(var);
+      if (!linear_propagator_->AddConstraint(
+              {}, tmp_vars_, new_ct.CoeffsAsSpan(), -new_ct.lb)) {
+        return false;
+      }
+    }
+  }
+  return true;
 }
 
 void LinearProgrammingConstraint::RegisterWith(Model* model) {
@@ -383,10 +427,8 @@ void LinearProgrammingConstraint::RegisterWith(Model* model) {
   orbit_indices_.clear();
   for (int i = 0; i < num_vars; i++) {
     const IntegerVariable pos_var = integer_variables_[i];
-    if (symmetrizer_->AppearInFoldedProblem(pos_var)) {
-      watcher_->WatchIntegerVariable(pos_var, watcher_id_, i);
-    }
-
+    DCHECK(symmetrizer_->AppearInFoldedProblem(pos_var));
+    watcher_->WatchIntegerVariable(pos_var, watcher_id_, i);
     if (symmetrizer_->IsOrbitSumVar(pos_var)) {
       orbit_indices_.push_back(symmetrizer_->OrbitIndex(pos_var));
     }
@@ -755,9 +797,9 @@ void LinearProgrammingConstraint::SetLevel(int level) {
     lp_solution_ = level_zero_lp_solution_;
     lp_solution_level_ = 0;
     for (int i = 0; i < lp_solution_.size(); i++) {
-      expanded_lp_solution_[integer_variables_[i]] = lp_solution_[i];
-      expanded_lp_solution_[NegationOf(integer_variables_[i])] =
-          -lp_solution_[i];
+      const IntegerVariable var = extended_integer_variables_[i];
+      expanded_lp_solution_[var] = lp_solution_[i];
+      expanded_lp_solution_[NegationOf(var)] = -lp_solution_[i];
     }
   }
 }
@@ -972,9 +1014,8 @@ bool LinearProgrammingConstraint::SolveLp() {
 
     // Compute integrality.
     lp_solution_is_integer_ = true;
-    for (int i = 0; i < num_vars; i++) {
-      if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) >
-          kCpEpsilon) {
+    for (const double value : lp_solution_) {
+      if (std::abs(value - std::round(value)) > kCpEpsilon) {
         lp_solution_is_integer_ = false;
         break;
       }
@@ -1244,9 +1285,9 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
   // in principle be removed. Easy for cuts, but not so much for
   // implied_bounds_processor_. Note that in theory this could allow us to
   // use Literal directly without the need to have an IntegerVariable for them.
-  tmp_scattered_vector_.ConvertToCutData(cut_ub.value(), integer_variables_,
-                                         lp_solution_, integer_trail_,
-                                         &base_ct_);
+  tmp_scattered_vector_.ConvertToCutData(
+      cut_ub.value(), extended_integer_variables_, lp_solution_, integer_trail_,
+      &base_ct_);
 
   // If there are no integer (all Booleans), no need to try implied bounds
   // heurititics. By setting this to nullptr, we are a bit faster.
@@ -1425,7 +1466,7 @@ bool LinearProgrammingConstraint::PostprocessAndAddCut(
   }
 
   // Substitute any slack left.
-  tmp_scattered_vector_.ClearAndResize(integer_variables_.size());
+  tmp_scattered_vector_.ClearAndResize(extended_integer_variables_.size());
   IntegerValue cut_ub = static_cast<int64_t>(cut.rhs);
   for (const CutTerm& term : cut.terms) {
     if (term.coeff == 0) continue;
@@ -1474,8 +1515,8 @@ bool LinearProgrammingConstraint::PostprocessAndAddCut(
   // something we already added. This tends to happen if the duplicate was
   // already a generated cut which is currently not part of the LP.
   LinearConstraint converted_cut =
-      tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_,
-                                                      cut_ub);
+      tmp_scattered_vector_.ConvertToLinearConstraint(
+          extended_integer_variables_, cut_ub);
 
   // TODO(user): We probably generate too many cuts, keep best one only?
   // Note that we do need duplicate removal and maybe some orthogonality here?
@@ -2257,7 +2298,7 @@ bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
     ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound) const {
   // Initialize the new constraint.
   *upper_bound = 0;
-  scattered_vector->ClearAndResize(integer_variables_.size());
+  scattered_vector->ClearAndResize(extended_integer_variables_.size());
 
   // Compute the new constraint by taking the linear combination given by
   // integer_multipliers of the integer constraints in integer_lp_.
@@ -2565,14 +2606,19 @@ bool LinearProgrammingConstraint::PropagateExactLpReason() {
     extra_term = {objective_cp_, -obj_scale};
   }
 
-  AdjustNewLinearConstraint(&tmp_integer_multipliers_, &tmp_scattered_vector_,
-                            &rc_ub);
+  // TODO(user): It seems when the LP as a single variable and the equation is
+  // obj >= lower_bound, this removes it ? For now we disable this if the
+  // objective variable is part of the LP (i.e. single objective).
+  if (take_objective_into_account) {
+    AdjustNewLinearConstraint(&tmp_integer_multipliers_, &tmp_scattered_vector_,
+                              &rc_ub);
+  }
 
   // Create the IntegerSumLE that will allow to propagate the objective and more
   // generally do the reduced cost fixing.
   LinearConstraint explanation =
-      tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, rc_ub,
-                                                      extra_term);
+      tmp_scattered_vector_.ConvertToLinearConstraint(
+          extended_integer_variables_, rc_ub, extra_term);
 
   // Corner case where prevent overflow removed all terms.
   if (explanation.num_terms == 0) {
@@ -2612,8 +2658,8 @@ bool LinearProgrammingConstraint::PropagateExactDualRay() {
                             &new_constraint_ub);
 
   LinearConstraint explanation =
-      tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_,
-                                                      new_constraint_ub);
+      tmp_scattered_vector_.ConvertToLinearConstraint(
+          extended_integer_variables_, new_constraint_ub);
 
   std::string message;
   if (VLOG_IS_ON(1)) {
diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h
index 715de7ba38..093e68de9a 100644
--- a/ortools/sat/linear_programming_constraint.h
+++ b/ortools/sat/linear_programming_constraint.h
@@ -145,7 +145,8 @@ class LinearProgrammingConstraint : public PropagatorInterface,
                               absl::Span<const IntegerVariable> vars);
 
   // Add a new linear constraint to this LP.
-  void AddLinearConstraint(LinearConstraint ct);
+  // Return false if we prove infeasibility of the global model.
+  bool AddLinearConstraint(LinearConstraint ct);
 
   // Set the coefficient of the variable in the objective. Calling it twice will
   // overwrite the previous value.
@@ -183,12 +184,14 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // ReversibleInterface API.
   void SetLevel(int level) override;
 
+  // From outside, the lp should be seen as containing all extended variables.
   int NumVariables() const {
-    return static_cast<int>(integer_variables_.size());
+    return static_cast<int>(extended_integer_variables_.size());
   }
   const std::vector<IntegerVariable>& integer_variables() const {
-    return integer_variables_;
+    return extended_integer_variables_;
   }
+
   std::string DimensionString() const;
 
   // Returns a IntegerLiteral guided by the underlying LP constraints.
@@ -453,6 +456,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
 
   std::vector<glop::ColIndex> tmp_cols_;
   std::vector<IntegerValue> tmp_coeffs_;
+  std::vector<IntegerVariable> tmp_vars_;
 
   LinearExpression integer_objective_;
   IntegerValue integer_objective_offset_ = IntegerValue(0);
@@ -501,6 +505,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // TODO(user): This should be util_intops::StrongVector<glop::ColIndex,
   // IntegerVariable>.
   std::vector<IntegerVariable> integer_variables_;
+  std::vector<IntegerVariable> extended_integer_variables_;
   absl::flat_hash_map<IntegerVariable, glop::ColIndex> mirror_lp_variable_;
 
   // This is only used if we use symmetry folding.
@@ -530,6 +535,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   SharedResponseManager* shared_response_manager_;
   ModelRandomGenerator* random_;
   LinearConstraintSymmetrizer* symmetrizer_;
+  LinearPropagator* linear_propagator_;
 
   int watcher_id_;
 
diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc
index 5175b818c7..8d306f3b59 100644
--- a/ortools/sat/linear_relaxation.cc
+++ b/ortools/sat/linear_relaxation.cc
@@ -1968,30 +1968,6 @@ LinearRelaxation ComputeLinearRelaxation(const CpModelProto& model_proto,
           [](const LinearConstraint& lc) { return lc.num_terms <= 1; }),
       relaxation.linear_constraints.end());
 
-  // We add a clique cut generation over all Booleans of the problem.
-  // Note that in practice this might regroup independent LP together.
-  //
-  // TODO(user): compute connected components of the original problem and
-  // split these cuts accordingly.
-  if (params.linearization_level() > 1 && params.add_clique_cuts()) {
-    LinearConstraintBuilder builder(m);
-    for (int i = 0; i < model_proto.variables_size(); ++i) {
-      if (!mapping->IsBoolean(i)) continue;
-
-      // Note that it is okay to simply ignore the literal if it has no
-      // integer view.
-      const bool unused ABSL_ATTRIBUTE_UNUSED =
-          builder.AddLiteralTerm(mapping->Literal(i), IntegerValue(1));
-    }
-
-    // We add a generator touching all the variable in the builder.
-    const LinearExpression& expr = builder.BuildExpression();
-    if (!expr.vars.empty()) {
-      relaxation.cut_generators.push_back(
-          CreateCliqueCutGenerator(expr.vars, m));
-    }
-  }
-
   if (!relaxation.linear_constraints.empty() ||
       !relaxation.cut_generators.empty()) {
     SOLVER_LOG(logger,
diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc
index d844d2f08a..738297267d 100644
--- a/ortools/sat/presolve_context.cc
+++ b/ortools/sat/presolve_context.cc
@@ -745,6 +745,10 @@ void PresolveContext::UpdateNewConstraintsVariableUsage() {
 bool PresolveContext::HasUnusedAffineVariable() const {
   if (is_unsat_) return false;  // We do not care in this case.
   if (keep_all_feasible_solutions) return false;
+
+  // We can leave non-optimal stuff around if we reach the time limit.
+  if (time_limit_->LimitReached()) return false;
+
   for (int var = 0; var < working_model->variables_size(); ++var) {
     if (VariableIsNotUsedAnymore(var)) continue;
     if (IsFixed(var)) continue;
@@ -1584,12 +1588,9 @@ bool PresolveContext::InsertVarValueEncoding(int literal, int var,
                                       /*add_constraints=*/true)) {
     return false;
   }
-  if (!StoreLiteralImpliesVarEqValue(literal, var, value)) {
-    return false;
-  }
-  if (!StoreLiteralImpliesVarNEqValue(NegatedRef(literal), var, value)) {
-    return false;
-  }
+
+  eq_half_encoding_.insert({{literal, var}, value});
+  neq_half_encoding_.insert({{NegatedRef(literal), var}, value});
 
   if (hint_is_loaded_) {
     const int bool_var = PositiveRef(literal);
diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc
index fee8a20463..edbec59f10 100644
--- a/ortools/sat/synchronization.cc
+++ b/ortools/sat/synchronization.cc
@@ -25,6 +25,7 @@
 #include <deque>
 #include <functional>
 #include <limits>
+#include <memory>
 #include <string>
 #include <utility>
 #include <vector>
@@ -48,6 +49,7 @@
 #include "absl/strings/string_view.h"
 #include "absl/synchronization/mutex.h"
 #include "absl/types/span.h"
+#include "ortools/algorithms/sparse_permutation.h"
 #include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/cp_model_utils.h"
 #include "ortools/sat/integer.h"
@@ -55,6 +57,7 @@
 #include "ortools/sat/sat_base.h"
 #include "ortools/sat/sat_parameters.pb.h"
 #include "ortools/sat/sat_solver.h"
+#include "ortools/sat/symmetry_util.h"
 #include "ortools/sat/util.h"
 #include "ortools/util/bitset.h"
 #include "ortools/util/logging.h"
@@ -842,6 +845,45 @@ SharedBoundsManager::SharedBoundsManager(const CpModelProto& model_proto)
     synchronized_lower_bounds_[i] = lower_bounds_[i];
     synchronized_upper_bounds_[i] = upper_bounds_[i];
   }
+
+  // Fill symmetry data.
+  if (model_proto.has_symmetry()) {
+    const int num_vars = model_proto.variables().size();
+    std::vector<std::unique_ptr<SparsePermutation>> generators;
+    for (const SparsePermutationProto& perm :
+         model_proto.symmetry().permutations()) {
+      generators.emplace_back(CreateSparsePermutationFromProto(num_vars, perm));
+    }
+    if (generators.empty()) return;
+
+    // Get orbits in term of IntegerVariable.
+    var_to_orbit_index_ = GetOrbits(num_vars, generators);
+
+    // Fill orbits_.
+    std::vector<int> keys;
+    std::vector<int> values;
+    for (int var = 0; var < num_vars; ++var) {
+      const int orbit_index = var_to_orbit_index_[var];
+      if (orbit_index == -1) continue;
+      keys.push_back(orbit_index);
+      values.push_back(var);
+    }
+    if (keys.empty()) return;
+
+    has_symmetry_ = true;
+    orbits_.ResetFromFlatMapping(keys, values);
+
+    // Fill representative.
+    var_to_representative_.resize(num_vars);
+    for (int var = 0; var < num_vars; ++var) {
+      const int orbit_index = var_to_orbit_index_[var];
+      if (orbit_index == -1) {
+        var_to_representative_[var] = var;
+      } else {
+        var_to_representative_[var] = orbits_[orbit_index][0];
+      }
+    }
+  }
 }
 
 void SharedBoundsManager::ReportPotentialNewBounds(
@@ -851,11 +893,17 @@ void SharedBoundsManager::ReportPotentialNewBounds(
   CHECK_EQ(variables.size(), new_lower_bounds.size());
   CHECK_EQ(variables.size(), new_upper_bounds.size());
   int num_improvements = 0;
+  int num_symmetric_improvements = 0;
 
   absl::MutexLock mutex_lock(&mutex_);
   for (int i = 0; i < variables.size(); ++i) {
-    const int var = variables[i];
+    int var = variables[i];
     if (var >= num_variables_) continue;
+
+    // In the presence of symmetry we only update the representative.
+    if (has_symmetry_) {
+      var = var_to_representative_[var];
+    }
     const int64_t old_lb = lower_bounds_[var];
     const int64_t old_ub = upper_bounds_[var];
     const int64_t new_lb = new_lower_bounds[i];
@@ -881,18 +929,29 @@ void SharedBoundsManager::ReportPotentialNewBounds(
     }
     changed_variables_since_last_synchronize_.Set(var);
     num_improvements++;
+
+    if (has_symmetry_ && variables[i] != var) {
+      // We count -1 so that num_improvements + num_symmetric_improvements
+      // corresponds to the number of actual bound improvement.
+      num_symmetric_improvements +=
+          orbits_[var_to_orbit_index_[var]].size() - 1;
+    }
   }
   if (num_improvements > 0) {
     total_num_improvements_ += num_improvements;
     VLOG(3) << total_num_improvements_ << "/" << num_variables_;
-    bounds_exported_[worker_name] += num_improvements;
+    bounds_exported_[worker_name].num_exported += num_improvements;
+    bounds_exported_[worker_name].num_symmetric += num_symmetric_improvements;
     if (absl::GetFlag(FLAGS_cp_model_dump_tightened_models)) {
       CpModelProto tight_model = model_proto_;
       for (int i = 0; i < num_variables_; ++i) {
         IntegerVariableProto* var_proto = tight_model.mutable_variables(i);
-        const Domain domain =
-            ReadDomainFromProto(*var_proto)
-                .IntersectionWith(Domain(lower_bounds_[i], upper_bounds_[i]));
+
+        int rep = i;
+        if (has_symmetry_) rep = var_to_representative_[i];
+        const Domain domain = ReadDomainFromProto(*var_proto)
+                                  .IntersectionWith(Domain(lower_bounds_[rep],
+                                                           upper_bounds_[rep]));
         FillDomainInProto(domain, var_proto);
       }
       const std::string filename = absl::StrCat(dump_prefix_, "tighened_model_",
@@ -910,6 +969,8 @@ void SharedBoundsManager::ReportPotentialNewBounds(
 void SharedBoundsManager::FixVariablesFromPartialSolution(
     const std::vector<int64_t>& solution,
     const std::vector<int>& variables_to_fix) {
+  // This function shouldn't be called if we has symmetry.
+  CHECK(!has_symmetry_);
   absl::MutexLock mutex_lock(&mutex_);
 
   // Abort if incompatible. Note that we only check the position that we are
@@ -957,6 +1018,7 @@ void SharedBoundsManager::Synchronize() {
   absl::MutexLock mutex_lock(&mutex_);
   for (const int var :
        changed_variables_since_last_synchronize_.PositionsSetAtLeastOnce()) {
+    DCHECK(!has_symmetry_ || var_to_representative_[var] == var);
     synchronized_lower_bounds_[var] = lower_bounds_[var];
     synchronized_upper_bounds_[var] = upper_bounds_[var];
     for (int j = 0; j < id_to_changed_variables_.size(); ++j) {
@@ -977,6 +1039,7 @@ int SharedBoundsManager::RegisterNewId() {
     const int64_t ub = model_proto_.variables(var).domain(domain_size - 1);
     if (lb != synchronized_lower_bounds_[var] ||
         ub != synchronized_upper_bounds_[var]) {
+      DCHECK(!has_symmetry_ || var_to_representative_[var] == var);
       id_to_changed_variables_[id].Set(var);
     }
   }
@@ -990,19 +1053,46 @@ void SharedBoundsManager::GetChangedBounds(
   new_lower_bounds->clear();
   new_upper_bounds->clear();
 
-  absl::MutexLock mutex_lock(&mutex_);
-  for (const int var : id_to_changed_variables_[id].PositionsSetAtLeastOnce()) {
-    variables->push_back(var);
+  {
+    absl::MutexLock mutex_lock(&mutex_);
+    for (const int var :
+         id_to_changed_variables_[id].PositionsSetAtLeastOnce()) {
+      DCHECK(!has_symmetry_ || var_to_representative_[var] == var);
+      variables->push_back(var);
+    }
+    id_to_changed_variables_[id].ClearAll();
+
+    // We need to report the bounds in a deterministic order as it is difficult
+    // to guarantee that nothing depend on the order in which the new bounds are
+    // processed.
+    absl::c_sort(*variables);
+    for (const int var : *variables) {
+      new_lower_bounds->push_back(synchronized_lower_bounds_[var]);
+      new_upper_bounds->push_back(synchronized_upper_bounds_[var]);
+    }
   }
-  id_to_changed_variables_[id].ClearAll();
 
-  // We need to report the bounds in a deterministic order as it is difficult to
-  // guarantee that nothing depend on the order in which the new bounds are
-  // processed.
-  absl::c_sort(*variables);
-  for (const int var : *variables) {
-    new_lower_bounds->push_back(synchronized_lower_bounds_[var]);
-    new_upper_bounds->push_back(synchronized_upper_bounds_[var]);
+  // Now that the mutex is released, we can add all symmetric version if any.
+  // Note that alternatively we could do that in the client side, but the
+  // complexity will be the same, we will just save some memory that is usually
+  // just reused.
+  if (has_symmetry_) {
+    const int old_size = variables->size();
+    for (int i = 0; i < old_size; ++i) {
+      const int var = (*variables)[i];
+      const int orbit_index = var_to_orbit_index_[var];
+      if (orbit_index == -1) continue;
+
+      const int64_t lb = (*new_lower_bounds)[i];
+      const int64_t ub = (*new_upper_bounds)[i];
+      const auto orbit = orbits_[orbit_index];
+      CHECK_EQ(var, orbit[0]);
+      for (const int other : orbit.subspan(1)) {
+        variables->push_back(other);
+        new_lower_bounds->push_back(lb);
+        new_upper_bounds->push_back(ub);
+      }
+    }
   }
 }
 
@@ -1019,9 +1109,11 @@ void SharedBoundsManager::LogStatistics(SolverLogger* logger) {
   absl::MutexLock mutex_lock(&mutex_);
   if (!bounds_exported_.empty()) {
     std::vector<std::vector<std::string>> table;
-    table.push_back({"Improving bounds shared", "Num"});
+    table.push_back({"Improving bounds shared", "Num", "Sym"});
     for (const auto& entry : bounds_exported_) {
-      table.push_back({FormatName(entry.first), FormatCounter(entry.second)});
+      table.push_back({FormatName(entry.first),
+                       FormatCounter(entry.second.num_exported),
+                       FormatCounter(entry.second.num_symmetric)});
     }
     SOLVER_LOG(logger, FormatTable(table));
   }
@@ -1031,7 +1123,7 @@ int SharedBoundsManager::NumBoundsExported(const std::string& worker_name) {
   absl::MutexLock mutex_lock(&mutex_);
   const auto it = bounds_exported_.find(worker_name);
   if (it == bounds_exported_.end()) return 0;
-  return it->second;
+  return it->second.num_exported;
 }
 
 UniqueClauseStream::UniqueClauseStream() {
diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h
index c5d5acc116..2b07f65448 100644
--- a/ortools/sat/synchronization.h
+++ b/ortools/sat/synchronization.h
@@ -505,6 +505,9 @@ class SharedBoundsManager {
   // Note that because there can be more than one optimal solution on an
   // independent subproblem, it is important to do that in a locked fashion, and
   // reject future incompatible fixing.
+  //
+  // Note that this do not work with symmetries. And for now we don't call it
+  // when this is the case.
   void FixVariablesFromPartialSolution(
       const std::vector<int64_t>& solution,
       const std::vector<int>& variables_to_fix);
@@ -561,7 +564,21 @@ class SharedBoundsManager {
   std::vector<int64_t> synchronized_upper_bounds_ ABSL_GUARDED_BY(mutex_);
   std::deque<SparseBitset<int>> id_to_changed_variables_
       ABSL_GUARDED_BY(mutex_);
-  absl::btree_map<std::string, int> bounds_exported_ ABSL_GUARDED_BY(mutex_);
+
+  // We track the number of bounds exported by each solver, and the "extra"
+  // bounds pushed due to symmetries.
+  struct Counters {
+    int64_t num_exported = 0;
+    int64_t num_symmetric = 0;
+  };
+  absl::btree_map<std::string, Counters> bounds_exported_
+      ABSL_GUARDED_BY(mutex_);
+
+  // Symmetry info.
+  bool has_symmetry_ = false;
+  std::vector<int> var_to_representative_;  // Identity if not touched.
+  std::vector<int> var_to_orbit_index_;
+  CompactVectorVector<int, int> orbits_;
 
   std::vector<int64_t> debug_solution_;
   std::string dump_prefix_;