diff --git a/ortools/constraint_solver/constraint_solveri.h b/ortools/constraint_solver/constraint_solveri.h
index 4bf288af5c..6243c1793f 100644
--- a/ortools/constraint_solver/constraint_solveri.h
+++ b/ortools/constraint_solver/constraint_solveri.h
@@ -1862,7 +1862,8 @@ class LocalSearchState {
   DEFINE_STRONG_INT_TYPE(ConstraintId, int);
   // Adds a variable domain to this state, returns a handler to the new domain.
   VariableDomainId AddVariableDomain(int64_t relaxed_min, int64_t relaxed_max);
-  void RelaxVariableDomain(VariableDomainId domain_id);
+  // Relaxes the domain, returns false iff the domain was already relaxed.
+  bool RelaxVariableDomain(VariableDomainId domain_id);
   bool TightenVariableDomainMin(VariableDomainId domain_id, int64_t value);
   bool TightenVariableDomainMax(VariableDomainId domain_id, int64_t value);
   int64_t VariableDomainMin(VariableDomainId domain_id) const;
@@ -2105,8 +2106,9 @@ class LocalSearchState::Variable {
   }
   void Relax() const {
     if (state_ == nullptr) return;
-    state_->RelaxVariableDomain(domain_id_);
-    state_->PropagateRelax(domain_id_);
+    if (state_->RelaxVariableDomain(domain_id_)) {
+      state_->PropagateRelax(domain_id_);
+    }
   }
   bool Exists() const { return state_ != nullptr; }
 
diff --git a/ortools/constraint_solver/local_search.cc b/ortools/constraint_solver/local_search.cc
index e15bce5b43..ec33c28127 100644
--- a/ortools/constraint_solver/local_search.cc
+++ b/ortools/constraint_solver/local_search.cc
@@ -3282,7 +3282,7 @@ LocalSearchState::Variable LocalSearchState::DummyVariable() {
   return {nullptr, VariableDomainId(-1)};
 }
 
-void LocalSearchState::RelaxVariableDomain(VariableDomainId domain_id) {
+bool LocalSearchState::RelaxVariableDomain(VariableDomainId domain_id) {
   DCHECK(state_domains_are_all_nonempty_);
   if (!state_has_relaxed_domains_) {
     trailed_num_committed_empty_domains_ = num_committed_empty_domains_;
@@ -3297,7 +3297,9 @@ void LocalSearchState::RelaxVariableDomain(VariableDomainId domain_id) {
       --num_committed_empty_domains_;
     }
     current_domains_[domain_id] = relaxed_domains_[domain_id];
+    return true;
   }
+  return false;
 }
 
 int64_t LocalSearchState::VariableDomainMin(VariableDomainId domain_id) const {
diff --git a/ortools/routing/decision_builders.cc b/ortools/routing/decision_builders.cc
index e6c3b01de7..a8ea4f2f6c 100644
--- a/ortools/routing/decision_builders.cc
+++ b/ortools/routing/decision_builders.cc
@@ -17,6 +17,7 @@
 #include <cstdint>
 #include <functional>
 #include <limits>
+#include <new>
 #include <string>
 #include <utility>
 #include <vector>
@@ -356,10 +357,10 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder {
       std::vector<int64_t>* break_start_end_values) {
     cumul_values->clear();
     break_start_end_values->clear();
-    const RouteDimensionTravelInfo& dimension_travel_info =
+    const RouteDimensionTravelInfo* const dimension_travel_info =
         dimension_travel_info_per_route_.empty()
-            ? RouteDimensionTravelInfo()
-            : dimension_travel_info_per_route_[vehicle];
+            ? nullptr
+            : &dimension_travel_info_per_route_[vehicle];
     const Resource* resource = nullptr;
     if (rg_index_ >= 0 && model_.ResourceVar(vehicle, rg_index_)->Bound()) {
       const int resource_index =
diff --git a/ortools/routing/filters.cc b/ortools/routing/filters.cc
index 03f0fb1c56..9eedfea741 100644
--- a/ortools/routing/filters.cc
+++ b/ortools/routing/filters.cc
@@ -43,7 +43,6 @@
 #include "ortools/base/map_util.h"
 #include "ortools/base/small_map.h"
 #include "ortools/base/strong_vector.h"
-#include "ortools/base/types.h"
 #include "ortools/constraint_solver/constraint_solver.h"
 #include "ortools/constraint_solver/constraint_solveri.h"
 #include "ortools/routing/lp_scheduling.h"
@@ -230,38 +229,43 @@ class ActiveNodeGroupFilter : public IntVarLocalSearchFilter {
   explicit ActiveNodeGroupFilter(const RoutingModel& routing_model)
       : IntVarLocalSearchFilter(routing_model.Nexts()),
         routing_model_(routing_model),
-        group_is_active_(routing_model.GetSameActivityGroupsCount(), false) {}
+        active_count_per_group_(routing_model.GetSameActivityGroupsCount(),
+                                {.active = 0, .unknown = 0}),
+        node_is_active_(routing_model.Nexts().size(), false),
+        node_is_unknown_(routing_model.Nexts().size(), false) {}
+
   bool Accept(const Assignment* delta, const Assignment* /*deltadelta*/,
               int64_t /*objective_min*/, int64_t /*objective_max*/) override {
-    absl::flat_hash_map<int, int> active_count_per_group;
-    absl::flat_hash_map<int, int> lns_count_per_group;
+    active_count_per_group_.Revert();
     const Assignment::IntContainer& container = delta->IntVarContainer();
     for (const IntVarElement& new_element : container.elements()) {
       IntVar* const var = new_element.Var();
       int64_t index = -1;
       if (!FindIndex(var, &index)) continue;
+      const int group = routing_model_.GetSameActivityGroupOfIndex(index);
+      ActivityCounts counts = active_count_per_group_.Get(group);
+      // Change contribution to counts: remove old state, add new state.
+      if (node_is_unknown_[index]) --counts.unknown;
+      if (node_is_active_[index]) --counts.active;
       if (new_element.Min() != new_element.Max()) {
-        // LNS detected.
-        lns_count_per_group[routing_model_.GetSameActivityGroupOfIndex(
-            index)]++;
-        continue;
-      }
-      int& active_count = gtl::LookupOrInsert(
-          &active_count_per_group,
-          routing_model_.GetSameActivityGroupOfIndex(index), 0);
-      if (new_element.Min() != index) {
-        if (active_count < 0) return false;
-        ++active_count;
-      } else {
-        if (active_count > 0) return false;
-        --active_count;
+        ++counts.unknown;
+      } else if (new_element.Min() != index) {
+        ++counts.active;
       }
+      active_count_per_group_.Set(group, counts);
     }
-    for (const auto& [group, count] : active_count_per_group) {
+    for (const int group : active_count_per_group_.ChangedIndices()) {
+      const ActivityCounts counts = active_count_per_group_.Get(group);
       const int group_size =
           routing_model_.GetSameActivityIndicesOfGroup(group).size();
-      if (std::abs(count) + lns_count_per_group[group] == group_size) continue;
-      if ((count > 0) != group_is_active_[group]) return false;
+      // The group constraint is respected iff either 0 or group size is inside
+      // interval [num_active, num_active + num_unknown],
+      if (counts.active == 0) continue;
+      if (counts.active <= group_size &&
+          group_size <= counts.active + counts.unknown) {
+        continue;
+      }
+      return false;
     }
     return true;
   }
@@ -269,28 +273,38 @@ class ActiveNodeGroupFilter : public IntVarLocalSearchFilter {
 
  private:
   void OnSynchronize(const Assignment* /*delta*/) override {
-    for (int group = 0; group < routing_model_.GetSameActivityGroupsCount();
-         ++group) {
-      bool is_group_active = false;
+    const int num_groups = routing_model_.GetSameActivityGroupsCount();
+    for (int group = 0; group < num_groups; ++group) {
+      ActivityCounts counts = {.active = 0, .unknown = 0};
       for (int node : routing_model_.GetSameActivityIndicesOfGroup(group)) {
         if (IsVarSynced(node)) {
-          is_group_active = (Value(node) != node);
-          break;
+          const bool is_active = (Value(node) != node);
+          node_is_active_[node] = is_active;
+          node_is_unknown_[node] = false;
+          counts.active += is_active ? 1 : 0;
+        } else {
+          ++counts.unknown;
+          node_is_unknown_[node] = true;
+          node_is_active_[node] = false;
         }
       }
-#ifndef NDEBUG
-      for (int node : routing_model_.GetSameActivityIndicesOfGroup(group)) {
-        if (IsVarSynced(node)) {
-          DCHECK_EQ((Value(node) != node), is_group_active);
-        }
-      }
-#endif  // NDEBUG
-      group_is_active_[group] = is_group_active;
+      active_count_per_group_.Set(group, counts);
     }
+    active_count_per_group_.Commit();
   }
 
   const RoutingModel& routing_model_;
-  std::vector<bool> group_is_active_;
+  struct ActivityCounts {
+    int active;
+    int unknown;
+  };
+  CommittableVector<ActivityCounts> active_count_per_group_;
+  // node_is_active_[node] is true iff node was synced and active at last
+  // Synchronize().
+  std::vector<bool> node_is_active_;
+  // node_is_unknown_[node] is true iff node was not synced at last
+  // Synchronize().
+  std::vector<bool> node_is_unknown_;
 };
 
 }  // namespace
@@ -310,54 +324,53 @@ class NodeDisjunctionFilter : public IntVarLocalSearchFilter {
                                  bool filter_cost)
       : IntVarLocalSearchFilter(routing_model.Nexts()),
         routing_model_(routing_model),
-        active_per_disjunction_(routing_model.GetNumberOfDisjunctions(), 0),
-        inactive_per_disjunction_(routing_model.GetNumberOfDisjunctions(), 0),
+        count_per_disjunction_(routing_model.GetNumberOfDisjunctions(),
+                               {.active = 0, .inactive = 0}),
         synchronized_objective_value_(std::numeric_limits<int64_t>::min()),
         accepted_objective_value_(std::numeric_limits<int64_t>::min()),
         filter_cost_(filter_cost),
         has_mandatory_disjunctions_(routing_model.HasMandatoryDisjunctions()) {}
 
+  using Disjunction = RoutingModel::DisjunctionIndex;
+
   bool Accept(const Assignment* delta, const Assignment* /*deltadelta*/,
               int64_t /*objective_min*/, int64_t objective_max) override {
-    const int64_t kUnassigned = -1;
-    const Assignment::IntContainer& container = delta->IntVarContainer();
-    gtl::small_map<absl::flat_hash_map<RoutingModel::DisjunctionIndex, int>>
-        disjunction_active_deltas;
-    gtl::small_map<absl::flat_hash_map<RoutingModel::DisjunctionIndex, int>>
-        disjunction_inactive_deltas;
+    count_per_disjunction_.Revert();
     bool lns_detected = false;
-    // Update active/inactive count per disjunction for each element of delta.
-    for (const IntVarElement& new_element : container.elements()) {
-      IntVar* const var = new_element.Var();
-      int64_t index = kUnassigned;
-      if (FindIndex(var, &index)) {
-        const bool is_inactive =
-            (new_element.Min() <= index && new_element.Max() >= index);
-        if (new_element.Min() != new_element.Max()) {
-          lns_detected = true;
-        }
-        for (const RoutingModel::DisjunctionIndex disjunction_index :
-             routing_model_.GetDisjunctionIndices(index)) {
-          const bool is_var_synced = IsVarSynced(index);
-          if (!is_var_synced || (Value(index) == index) != is_inactive) {
-            ++gtl::LookupOrInsert(is_inactive ? &disjunction_inactive_deltas
-                                              : &disjunction_active_deltas,
-                                  disjunction_index, 0);
-            if (is_var_synced) {
-              --gtl::LookupOrInsert(is_inactive ? &disjunction_active_deltas
-                                                : &disjunction_inactive_deltas,
-                                    disjunction_index, 0);
-            }
-          }
-        }
+    // Update the active/inactive counts of each modified disjunction.
+    for (const IntVarElement& element : delta->IntVarContainer().elements()) {
+      int64_t node = -1;
+      if (!FindIndex(element.Var(), &node)) continue;
+      lns_detected |= element.Min() != element.Max();
+      // Compute difference in how this node contributes to activity counts.
+      const bool is_var_synced = IsVarSynced(node);
+      const bool was_active = is_var_synced && Value(node) != node;
+      const bool is_active = node < element.Min() || element.Max() < node;
+      ActivityCount contribution_delta = {.active = 0, .inactive = 0};
+      if (is_var_synced) {
+        contribution_delta.active -= was_active;
+        contribution_delta.inactive -= !was_active;
+      }
+      contribution_delta.active += is_active;
+      contribution_delta.inactive += !is_active;
+      // Common shortcut: if the change is neutral, counts stay the same.
+      if (contribution_delta.active == 0 && contribution_delta.inactive == 0) {
+        continue;
+      }
+      // Change counts of all disjunctions affected by this node.
+      for (const Disjunction disjunction :
+           routing_model_.GetDisjunctionIndices(node)) {
+        ActivityCount new_count =
+            count_per_disjunction_.Get(disjunction.value());
+        new_count.active += contribution_delta.active;
+        new_count.inactive += contribution_delta.inactive;
+        count_per_disjunction_.Set(disjunction.value(), new_count);
       }
     }
     // Check if any disjunction has too many active nodes.
-    for (const auto [disjunction_index, active_nodes_delta] :
-         disjunction_active_deltas) {
-      // Too many active nodes.
-      if (active_per_disjunction_[disjunction_index] + active_nodes_delta >
-          routing_model_.GetDisjunctionMaxCardinality(disjunction_index)) {
+    for (const int index : count_per_disjunction_.ChangedIndices()) {
+      if (count_per_disjunction_.Get(index).active >
+          routing_model_.GetDisjunctionMaxCardinality(Disjunction(index))) {
         return false;
       }
     }
@@ -367,70 +380,32 @@ class NodeDisjunctionFilter : public IntVarLocalSearchFilter {
     }
     // Update penalty costs for disjunctions.
     accepted_objective_value_ = synchronized_objective_value_;
-    for (const auto [disjunction_index, inactive_nodes_delta] :
-         disjunction_inactive_deltas) {
-      const int64_t penalty =
-          routing_model_.GetDisjunctionPenalty(disjunction_index);
+    for (const int index : count_per_disjunction_.ChangedIndices()) {
+      // If num inactives did not change, skip. Common shortcut.
+      const int old_inactives =
+          count_per_disjunction_.GetCommitted(index).inactive;
+      const int new_inactives = count_per_disjunction_.Get(index).inactive;
+      if (old_inactives == new_inactives) continue;
+      // If this disjunction has no penalty for inactive nodes, skip.
+      const Disjunction disjunction(index);
+      const int64_t penalty = routing_model_.GetDisjunctionPenalty(disjunction);
       if (penalty == 0) continue;
-      const int current_inactive_nodes =
-          inactive_per_disjunction_[disjunction_index];
-      const int max_inactive_cardinality =
-          routing_model_.GetDisjunctionNodeIndices(disjunction_index).size() -
-          routing_model_.GetDisjunctionMaxCardinality(disjunction_index);
-      // Too many inactive nodes.
-      const int inactive_nodes_above_limit =
-          (current_inactive_nodes + inactive_nodes_delta) -
-          max_inactive_cardinality;
-      if (inactive_nodes_above_limit > 0 && penalty < 0) {
-        // Nodes are mandatory, i.e. exactly max_cardinality nodes must be
-        // performed, so the move is not acceptable.
-        return false;
-      }
 
-      const RoutingModel::PenaltyCostBehavior penalty_cost_behavior =
-          routing_model_.GetDisjunctionPenaltyCostBehavior(disjunction_index);
-      switch (penalty_cost_behavior) {
-        case RoutingModel::PenaltyCostBehavior::PENALIZE_ONCE:
-          if (inactive_nodes_above_limit > 0) {  // penalty cost to update
-            if (current_inactive_nodes <= max_inactive_cardinality) {
-              // Add penalty if there were not too many inactive nodes before
-              // the move.
-              CapAddTo(penalty, &accepted_objective_value_);
-            }
-          } else {  // no more penalty cost
-            if (current_inactive_nodes > max_inactive_cardinality) {
-              // Remove penalty if there were too many inactive nodes before the
-              // move.
-              accepted_objective_value_ =
-                  CapSub(accepted_objective_value_, penalty);
-            }
-          }
-          break;
-        case RoutingModel::PenaltyCostBehavior::PENALIZE_PER_INACTIVE:
-          if (inactive_nodes_above_limit > 0) {  // penalty cost to update
-            if (current_inactive_nodes <= max_inactive_cardinality) {
-              // Add penalty if there were not too many inactive nodes before
-              // the move.
-              CapAddTo(penalty * inactive_nodes_above_limit,
-                       &accepted_objective_value_);
-            } else if (inactive_nodes_delta != 0) {
-              // Update penalty cost if there are new or fewer inactive nodes.
-              CapAddTo(penalty * inactive_nodes_delta,
-                       &accepted_objective_value_);
-            }
-          } else {  // no more penalty cost
-            if (current_inactive_nodes > max_inactive_cardinality) {
-              // Remove penalty if there were too many inactive nodes before the
-              // move.
-              const int current_inactive_nodes_above_limit =
-                  current_inactive_nodes - max_inactive_cardinality;
-              accepted_objective_value_ =
-                  CapSub(accepted_objective_value_,
-                         penalty * current_inactive_nodes_above_limit);
-            }
-          }
-          break;
+      // Compute the new cost of activity bound violations.
+      const int max_inactives =
+          routing_model_.GetDisjunctionNodeIndices(disjunction).size() -
+          routing_model_.GetDisjunctionMaxCardinality(disjunction);
+      int new_violation = std::max(0, new_inactives - max_inactives);
+      int old_violation = std::max(0, old_inactives - max_inactives);
+      // If nodes are mandatory, there can be no violation.
+      if (penalty < 0 && new_violation > 0) return false;
+      if (routing_model_.GetDisjunctionPenaltyCostBehavior(disjunction) ==
+          RoutingModel::PenaltyCostBehavior::PENALIZE_ONCE) {
+        new_violation = std::min(1, new_violation);
+        old_violation = std::min(1, old_violation);
       }
+      CapAddTo(CapProd(penalty, (new_violation - old_violation)),
+               &accepted_objective_value_);
     }
     // Only compare to max as a cost lower bound is computed.
     return accepted_objective_value_ <= objective_max;
@@ -446,51 +421,44 @@ class NodeDisjunctionFilter : public IntVarLocalSearchFilter {
  private:
   void OnSynchronize(const Assignment* /*delta*/) override {
     synchronized_objective_value_ = 0;
-    for (RoutingModel::DisjunctionIndex i(0);
-         i < active_per_disjunction_.size(); ++i) {
-      active_per_disjunction_[i] = 0;
-      inactive_per_disjunction_[i] = 0;
-      const std::vector<int64_t>& disjunction_indices =
-          routing_model_.GetDisjunctionNodeIndices(i);
-      for (const int64_t index : disjunction_indices) {
-        if (IsVarSynced(index)) {
-          if (Value(index) != index) {
-            ++active_per_disjunction_[i];
-          } else {
-            ++inactive_per_disjunction_[i];
-          }
-        }
+    count_per_disjunction_.Revert();
+    const int num_disjunctions = routing_model_.GetNumberOfDisjunctions();
+    for (Disjunction disjunction(0); disjunction < num_disjunctions;
+         ++disjunction) {
+      // Count number of active/inactive nodes of this disjunction.
+      ActivityCount count = {.active = 0, .inactive = 0};
+      const auto& nodes = routing_model_.GetDisjunctionNodeIndices(disjunction);
+      for (const int64_t node : nodes) {
+        if (!IsVarSynced(node)) continue;
+        const int is_active = Value(node) != node;
+        count.active += is_active;
+        count.inactive += !is_active;
       }
+      count_per_disjunction_.Set(disjunction.value(), count);
+      // Add penalty of this disjunction to total cost.
       if (!filter_cost_) continue;
-      const int64_t penalty = routing_model_.GetDisjunctionPenalty(i);
-      const int max_cardinality =
-          routing_model_.GetDisjunctionMaxCardinality(i);
-      const RoutingModel::PenaltyCostBehavior penalty_cost_behavior =
-          routing_model_.GetDisjunctionPenaltyCostBehavior(i);
-
-      const int inactive_nodes_above_limit =
-          inactive_per_disjunction_[i] -
-          (disjunction_indices.size() - max_cardinality);
-      if (inactive_nodes_above_limit > 0 && penalty > 0) {
-        switch (penalty_cost_behavior) {
-          case RoutingModel::PenaltyCostBehavior::PENALIZE_ONCE:
-            CapAddTo(penalty, &synchronized_objective_value_);
-            break;
-          case RoutingModel::PenaltyCostBehavior::PENALIZE_PER_INACTIVE:
-            CapAddTo(penalty * inactive_nodes_above_limit,
-                     &synchronized_objective_value_);
-            break;
+      const int64_t penalty = routing_model_.GetDisjunctionPenalty(disjunction);
+      const int max_actives =
+          routing_model_.GetDisjunctionMaxCardinality(disjunction);
+      int violation = count.inactive - (nodes.size() - max_actives);
+      if (violation > 0 && penalty > 0) {
+        if (routing_model_.GetDisjunctionPenaltyCostBehavior(disjunction) ==
+            RoutingModel::PenaltyCostBehavior::PENALIZE_ONCE) {
+          violation = std::min(1, violation);
         }
+        CapAddTo(CapProd(penalty, violation), &synchronized_objective_value_);
       }
     }
+    count_per_disjunction_.Commit();
+    accepted_objective_value_ = synchronized_objective_value_;
   }
 
   const RoutingModel& routing_model_;
-
-  util_intops::StrongVector<RoutingModel::DisjunctionIndex, int>
-      active_per_disjunction_;
-  util_intops::StrongVector<RoutingModel::DisjunctionIndex, int>
-      inactive_per_disjunction_;
+  struct ActivityCount {
+    int active = 0;
+    int inactive = 0;
+  };
+  CommittableVector<ActivityCount> count_per_disjunction_;
   int64_t synchronized_objective_value_;
   int64_t accepted_objective_value_;
   const bool filter_cost_;
@@ -1149,7 +1117,7 @@ class PathCumulFilter : public BasePathFilter {
   PathCumulFilter(const RoutingModel& routing_model,
                   const RoutingDimension& dimension,
                   bool propagate_own_objective_value,
-                  bool filter_objective_cost, bool can_use_lp);
+                  bool filter_objective_cost, bool may_use_optimizers);
   ~PathCumulFilter() override {}
   std::string DebugString() const override {
     return "PathCumulFilter(" + name_ + ")";
@@ -1163,7 +1131,7 @@ class PathCumulFilter : public BasePathFilter {
                : accepted_objective_value_;
   }
   bool UsesDimensionOptimizers() {
-    if (!can_use_lp_) return false;
+    if (!may_use_optimizers_) return false;
     for (int vehicle = 0; vehicle < routing_model_.vehicles(); ++vehicle) {
       if (FilterWithDimensionCumulOptimizerForVehicle(vehicle)) return true;
     }
@@ -1271,7 +1239,7 @@ class PathCumulFilter : public BasePathFilter {
   }
 
   bool FilterWithDimensionCumulOptimizerForVehicle(int vehicle) const {
-    if (!can_use_lp_ || FilterCumulPiecewiseLinearCosts()) {
+    if (!may_use_optimizers_ || FilterCumulPiecewiseLinearCosts()) {
       return false;
     }
 
@@ -1414,13 +1382,13 @@ class PathCumulFilter : public BasePathFilter {
   absl::btree_set<int> delta_paths_;
   const std::string name_;
 
-  LocalDimensionCumulOptimizer* optimizer_;
+  LocalDimensionCumulOptimizer* lp_optimizer_;
   LocalDimensionCumulOptimizer* mp_optimizer_;
   const std::function<int64_t(int64_t)> path_accessor_;
   const bool filter_objective_cost_;
   // This boolean indicates if the LP optimizer can be used if necessary to
   // optimize the dimension cumuls.
-  const bool can_use_lp_;
+  const bool may_use_optimizers_;
   const bool propagate_own_objective_value_;
 
   std::vector<int64_t> min_path_cumuls_;
@@ -1440,7 +1408,8 @@ std::vector<T> SumOfVectors(const std::vector<T>& v1,
 PathCumulFilter::PathCumulFilter(const RoutingModel& routing_model,
                                  const RoutingDimension& dimension,
                                  bool propagate_own_objective_value,
-                                 bool filter_objective_cost, bool can_use_lp)
+                                 bool filter_objective_cost,
+                                 bool may_use_optimizers)
     : BasePathFilter(routing_model.Nexts(), dimension.cumuls().size(),
                      routing_model.GetPathsMetadata()),
       routing_model_(routing_model),
@@ -1472,11 +1441,11 @@ PathCumulFilter::PathCumulFilter(const RoutingModel& routing_model,
       delta_max_end_cumul_(0),
       delta_nodes_with_precedences_and_changed_cumul_(routing_model.Size()),
       name_(dimension.name()),
-      optimizer_(routing_model.GetMutableLocalCumulLPOptimizer(dimension)),
+      lp_optimizer_(routing_model.GetMutableLocalCumulLPOptimizer(dimension)),
       mp_optimizer_(routing_model.GetMutableLocalCumulMPOptimizer(dimension)),
       path_accessor_([this](int64_t node) { return GetNext(node); }),
       filter_objective_cost_(filter_objective_cost),
-      can_use_lp_(can_use_lp),
+      may_use_optimizers_(may_use_optimizers),
       propagate_own_objective_value_(propagate_own_objective_value) {
   cumul_soft_bounds_.resize(cumuls_.size());
   cumul_soft_lower_bounds_.resize(cumuls_.size());
@@ -1555,7 +1524,7 @@ PathCumulFilter::PathCumulFilter(const RoutingModel& routing_model,
 #ifndef NDEBUG
   for (int vehicle = 0; vehicle < routing_model.vehicles(); vehicle++) {
     if (FilterWithDimensionCumulOptimizerForVehicle(vehicle)) {
-      DCHECK_NE(optimizer_, nullptr);
+      DCHECK_NE(lp_optimizer_, nullptr);
       DCHECK_NE(mp_optimizer_, nullptr);
     }
   }
@@ -1723,7 +1692,7 @@ void PathCumulFilter::OnBeforeSynchronizePaths() {
         LocalDimensionCumulOptimizer* const optimizer =
             (FilterSoftSpanQuadraticCost(vehicle) || FilterBreakCost(vehicle))
                 ? mp_optimizer_
-                : optimizer_;
+                : lp_optimizer_;
         DCHECK(optimizer != nullptr);
         const DimensionSchedulingStatus status =
             optimizer->ComputeRouteCumulCostWithoutFixedTransits(
@@ -2053,7 +2022,7 @@ bool PathCumulFilter::FinalizeAcceptPath(int64_t /*objective_min*/,
       CapAdd(cumul_cost_delta_, CapProd(global_span_cost_coefficient_,
                                         CapSub(new_max_end, new_min_start)));
 
-  if (can_use_lp_ && optimizer_ != nullptr &&
+  if (may_use_optimizers_ && lp_optimizer_ != nullptr &&
       accepted_objective_value_ <= objective_max) {
     const size_t num_touched_paths = GetTouchedPathStarts().size();
     std::vector<int64_t> path_delta_cost_values(num_touched_paths, 0);
@@ -2066,7 +2035,7 @@ bool PathCumulFilter::FinalizeAcceptPath(int64_t /*objective_min*/,
       }
       int64_t path_delta_cost_with_lp = 0;
       const DimensionSchedulingStatus status =
-          optimizer_->ComputeRouteCumulCostWithoutFixedTransits(
+          lp_optimizer_->ComputeRouteCumulCostWithoutFixedTransits(
               vehicle, path_accessor_, /*resource=*/nullptr,
               filter_objective_cost_ ? &path_delta_cost_with_lp : nullptr);
       if (status == DimensionSchedulingStatus::INFEASIBLE) {
@@ -2242,11 +2211,11 @@ int64_t PathCumulFilter::ComputePathMaxStartFromEndCumul(
 IntVarLocalSearchFilter* MakePathCumulFilter(const RoutingDimension& dimension,
                                              bool propagate_own_objective_value,
                                              bool filter_objective_cost,
-                                             bool can_use_lp) {
+                                             bool may_use_optimizers) {
   RoutingModel& model = *dimension.model();
   return model.solver()->RevAlloc(
       new PathCumulFilter(model, dimension, propagate_own_objective_value,
-                          filter_objective_cost, can_use_lp));
+                          filter_objective_cost, may_use_optimizers));
 }
 
 namespace {
@@ -2812,11 +2781,11 @@ class LPCumulFilter : public IntVarLocalSearchFilter {
   void OnSynchronize(const Assignment* delta) override;
   int64_t GetSynchronizedObjectiveValue() const override;
   std::string DebugString() const override {
-    return "LPCumulFilter(" + optimizer_.dimension()->name() + ")";
+    return "LPCumulFilter(" + lp_optimizer_.dimension()->name() + ")";
   }
 
  private:
-  GlobalDimensionCumulOptimizer& optimizer_;
+  GlobalDimensionCumulOptimizer& lp_optimizer_;
   GlobalDimensionCumulOptimizer& mp_optimizer_;
   const bool filter_objective_cost_;
   int64_t synchronized_cost_without_transit_;
@@ -2826,11 +2795,11 @@ class LPCumulFilter : public IntVarLocalSearchFilter {
 };
 
 LPCumulFilter::LPCumulFilter(const std::vector<IntVar*>& nexts,
-                             GlobalDimensionCumulOptimizer* optimizer,
+                             GlobalDimensionCumulOptimizer* lp_optimizer,
                              GlobalDimensionCumulOptimizer* mp_optimizer,
                              bool filter_objective_cost)
     : IntVarLocalSearchFilter(nexts),
-      optimizer_(*optimizer),
+      lp_optimizer_(*lp_optimizer),
       mp_optimizer_(*mp_optimizer),
       filter_objective_cost_(filter_objective_cost),
       synchronized_cost_without_transit_(-1),
@@ -2861,8 +2830,8 @@ bool LPCumulFilter::Accept(const Assignment* delta,
   if (!filter_objective_cost_) {
     // No need to compute the cost of the LP, only verify its feasibility.
     delta_cost_without_transit_ = 0;
-    const DimensionSchedulingStatus status =
-        optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr, nullptr);
+    const DimensionSchedulingStatus status = lp_optimizer_.ComputeCumuls(
+        next_accessor, {}, nullptr, nullptr, nullptr);
     if (status == DimensionSchedulingStatus::OPTIMAL) return true;
     if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY &&
         mp_optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr,
@@ -2874,7 +2843,7 @@ bool LPCumulFilter::Accept(const Assignment* delta,
   }
 
   const DimensionSchedulingStatus status =
-      optimizer_.ComputeCumulCostWithoutFixedTransits(
+      lp_optimizer_.ComputeCumulCostWithoutFixedTransits(
           next_accessor, &delta_cost_without_transit_);
   if (status == DimensionSchedulingStatus::INFEASIBLE) {
     delta_cost_without_transit_ = std::numeric_limits<int64_t>::max();
@@ -2899,7 +2868,7 @@ int64_t LPCumulFilter::GetAcceptedObjectiveValue() const {
 void LPCumulFilter::OnSynchronize(const Assignment* /*delta*/) {
   // TODO(user): Try to optimize this so the LP is not called when the last
   // computed delta cost corresponds to the solution being synchronized.
-  const RoutingModel& model = *optimizer_.dimension()->model();
+  const RoutingModel& model = *lp_optimizer_.dimension()->model();
   const auto& next_accessor = [this, &model](int64_t index) {
     return IsVarSynced(index)     ? Value(index)
            : model.IsStart(index) ? model.End(model.VehicleIndex(index))
@@ -2911,10 +2880,10 @@ void LPCumulFilter::OnSynchronize(const Assignment* /*delta*/) {
   }
   DimensionSchedulingStatus status =
       filter_objective_cost_
-          ? optimizer_.ComputeCumulCostWithoutFixedTransits(
+          ? lp_optimizer_.ComputeCumulCostWithoutFixedTransits(
                 next_accessor, &synchronized_cost_without_transit_)
-          : optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr,
-                                     nullptr);
+          : lp_optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr,
+                                        nullptr);
   if (status == DimensionSchedulingStatus::INFEASIBLE) {
     // TODO(user): This should only happen if the LP solver times out.
     // DCHECK the fail wasn't due to an infeasible model.
@@ -2941,13 +2910,13 @@ int64_t LPCumulFilter::GetSynchronizedObjectiveValue() const {
 }  // namespace
 
 IntVarLocalSearchFilter* MakeGlobalLPCumulFilter(
-    GlobalDimensionCumulOptimizer* optimizer,
+    GlobalDimensionCumulOptimizer* lp_optimizer,
     GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost) {
-  DCHECK_NE(optimizer, nullptr);
+  DCHECK_NE(lp_optimizer, nullptr);
   DCHECK_NE(mp_optimizer, nullptr);
-  const RoutingModel& model = *optimizer->dimension()->model();
+  const RoutingModel& model = *lp_optimizer->dimension()->model();
   return model.solver()->RevAlloc(new LPCumulFilter(
-      model.Nexts(), optimizer, mp_optimizer, filter_objective_cost));
+      model.Nexts(), lp_optimizer, mp_optimizer, filter_objective_cost));
 }
 
 namespace {
@@ -3175,11 +3144,11 @@ void ResourceGroupAssignmentFilter::OnSynchronizePathFromStart(int64_t start) {
     vehicle_to_resource_class_assignment_costs_[v] = {route_cost};
     return;
   }
-  // NOTE(user): Even if filter_objective_cost_ is false, we still need to
-  // call ComputeVehicleToResourceClassAssignmentCosts() for every vehicle
-  // requiring resource assignment to keep track of whether or not a given
-  // vehicle-to-resource-class assignment is possible by storing 0 or -1 in
-  // vehicle_to_resource_class_assignment_costs_.
+  // NOTE(user): Even if filter_objective_cost_ is false, we
+  // still need to call ComputeVehicleToResourceClassAssignmentCosts() for every
+  // vehicle requiring resource assignment to keep track of whether or not a
+  // given vehicle-to-resource-class assignment is possible by storing 0 or -1
+  // in vehicle_to_resource_class_assignment_costs_.
   if (!ComputeVehicleToResourceClassAssignmentCosts(
           v, resource_group_, ignored_resources_per_class_, next_accessor,
           dimension_.transit_evaluator(v), filter_objective_cost_,
@@ -3315,16 +3284,17 @@ class ResourceAssignmentFilter : public LocalSearchFilter {
 };
 
 ResourceAssignmentFilter::ResourceAssignmentFilter(
-    const std::vector<IntVar*>& nexts, LocalDimensionCumulOptimizer* optimizer,
+    const std::vector<IntVar*>& nexts,
+    LocalDimensionCumulOptimizer* lp_optimizer,
     LocalDimensionCumulOptimizer* mp_optimizer,
     bool propagate_own_objective_value, bool filter_objective_cost)
     : propagate_own_objective_value_(propagate_own_objective_value),
-      dimension_name_(optimizer->dimension()->name()) {
-  const RoutingModel& model = *optimizer->dimension()->model();
+      dimension_name_(lp_optimizer->dimension()->name()) {
+  const RoutingModel& model = *lp_optimizer->dimension()->model();
   for (const auto& resource_group : model.GetResourceGroups()) {
     resource_group_assignment_filters_.push_back(
         model.solver()->RevAlloc(new ResourceGroupAssignmentFilter(
-            nexts, resource_group.get(), optimizer, mp_optimizer,
+            nexts, resource_group.get(), lp_optimizer, mp_optimizer,
             filter_objective_cost)));
   }
 }
@@ -3361,14 +3331,14 @@ void ResourceAssignmentFilter::Synchronize(const Assignment* assignment,
 }  // namespace
 
 LocalSearchFilter* MakeResourceAssignmentFilter(
-    LocalDimensionCumulOptimizer* optimizer,
+    LocalDimensionCumulOptimizer* lp_optimizer,
     LocalDimensionCumulOptimizer* mp_optimizer,
     bool propagate_own_objective_value, bool filter_objective_cost) {
-  const RoutingModel& model = *optimizer->dimension()->model();
-  DCHECK_NE(optimizer, nullptr);
+  const RoutingModel& model = *lp_optimizer->dimension()->model();
+  DCHECK_NE(lp_optimizer, nullptr);
   DCHECK_NE(mp_optimizer, nullptr);
   return model.solver()->RevAlloc(new ResourceAssignmentFilter(
-      model.Nexts(), optimizer, mp_optimizer, propagate_own_objective_value,
+      model.Nexts(), lp_optimizer, mp_optimizer, propagate_own_objective_value,
       filter_objective_cost));
 }
 
@@ -4243,16 +4213,23 @@ bool LightVehicleBreaksChecker::Check() const {
     if (!path_data_[path].span.Exists()) continue;
     const int64_t total_transit = path_data_[path].total_transit.Min();
     // Compute lower bound of path span from break and path time windows.
+    const PathData& data = path_data_[path];
     int64_t lb_span_tw = total_transit;
-    const int64_t start_max = path_data_[path].start_cumul.Max();
-    const int64_t end_min = path_data_[path].end_cumul.Min();
-    for (const auto& br : path_data_[path].vehicle_breaks) {
+    const int64_t start_max = data.start_cumul.Max();
+    const int64_t end_min = data.end_cumul.Min();
+    for (const auto& br : data.vehicle_breaks) {
       if (!br.is_performed_min) continue;
       if (br.start_max < end_min && start_max < br.end_min) {
         CapAddTo(br.duration_min, &lb_span_tw);
       }
     }
-    if (!path_data_[path].span.SetMin(lb_span_tw)) return false;
+    if (!data.span.SetMin(lb_span_tw)) return false;
+    if (!data.start_cumul.SetMax(CapSub(data.end_cumul.Max(), lb_span_tw))) {
+      return false;
+    }
+    if (!data.end_cumul.SetMin(CapAdd(data.start_cumul.Min(), lb_span_tw))) {
+      return false;
+    }
   }
   return true;
 }
diff --git a/ortools/routing/filters.h b/ortools/routing/filters.h
index c88578bef5..ee332c8f81 100644
--- a/ortools/routing/filters.h
+++ b/ortools/routing/filters.h
@@ -14,6 +14,7 @@
 #ifndef OR_TOOLS_ROUTING_FILTERS_H_
 #define OR_TOOLS_ROUTING_FILTERS_H_
 
+#include <cstddef>
 #include <cstdint>
 #include <functional>
 #include <initializer_list>
@@ -22,6 +23,7 @@
 #include <utility>
 #include <vector>
 
+#include "absl/log/check.h"
 #include "absl/strings/string_view.h"
 #include "absl/types/span.h"
 #include "ortools/base/types.h"
@@ -36,6 +38,84 @@
 
 namespace operations_research::routing {
 
+// A vector that allows to revert back to a previously committed state,
+// get the set of changed indices, and get current and committed values.
+template <typename T>
+class CommittableVector {
+ public:
+  // Makes a vector with initial elements all committed to value.
+  CommittableVector<T>(size_t num_elements, const T& value)
+      : elements_(num_elements, {value, value}), changed_(num_elements) {}
+
+  // Return the size of the vector.
+  size_t Size() const { return elements_.size(); }
+
+  // Returns a copy of the value stored at index in the current state.
+  // Does not return a reference, because the class needs to know when elements
+  // are modified.
+  T Get(size_t index) const {
+    DCHECK_LT(index, elements_.size());
+    return elements_[index].current;
+  }
+
+  // Set the value stored at index in the current state to given value.
+  void Set(size_t index, const T& value) {
+    DCHECK_GE(index, 0);
+    DCHECK_LT(index, elements_.size());
+    changed_.Set(index);
+    elements_[index].current = value;
+  }
+
+  // Changes the values of the vector to those in the last Commit().
+  void Revert() {
+    for (const size_t index : changed_.PositionsSetAtLeastOnce()) {
+      elements_[index].current = elements_[index].committed;
+    }
+    changed_.ClearAll();
+  }
+
+  // Makes the current state committed, clearing all changes.
+  void Commit() {
+    for (const size_t index : changed_.PositionsSetAtLeastOnce()) {
+      elements_[index].committed = elements_[index].current;
+    }
+    changed_.ClearAll();
+  }
+
+  // Sets all elements of this vector to given value, and commits to this state.
+  // Supposes that there are no changes since the last Commit() or Revert().
+  void SetAllAndCommit(const T& value) {
+    DCHECK_EQ(0, changed_.NumberOfSetCallsWithDifferentArguments());
+    elements_.assign(elements_.size(), {value, value});
+  }
+
+  // Returns a copy of the value stored at index in the last committed state.
+  T GetCommitted(size_t index) const {
+    DCHECK_LT(index, elements_.size());
+    return elements_[index].committed;
+  }
+
+  // Return true iff the value at index has been Set() since the last Commit()
+  // or Revert(), even if the current value is the same as the committed value.
+  bool HasChanged(size_t index) const { return changed_[index]; }
+
+  // Returns the set of indices that have been Set() since the last Commit() or
+  // Revert().
+  const std::vector<size_t>& ChangedIndices() const {
+    return changed_.PositionsSetAtLeastOnce();
+  }
+
+ private:
+  struct VersionedElement {
+    T current;
+    T committed;
+  };
+  // Holds current and committed versions of values of this vector.
+  std::vector<VersionedElement> elements_;
+  // Holds indices that were Set() since the last Commit() or Revert().
+  SparseBitset<size_t> changed_;
+};
+
 /// Returns a filter tracking route constraints.
 IntVarLocalSearchFilter* MakeRouteConstraintFilter(
     const RoutingModel& routing_model);
@@ -76,7 +156,7 @@ IntVarLocalSearchFilter* MakeVehicleVarFilter(
 IntVarLocalSearchFilter* MakePathCumulFilter(const RoutingDimension& dimension,
                                              bool propagate_own_objective_value,
                                              bool filter_objective_cost,
-                                             bool can_use_lp);
+                                             bool may_use_optimizers);
 
 /// Returns a filter handling dimension cumul bounds.
 IntVarLocalSearchFilter* MakeCumulBoundsPropagatorFilter(
@@ -84,7 +164,7 @@ IntVarLocalSearchFilter* MakeCumulBoundsPropagatorFilter(
 
 /// Returns a filter checking global linear constraints and costs.
 IntVarLocalSearchFilter* MakeGlobalLPCumulFilter(
-    GlobalDimensionCumulOptimizer* optimizer,
+    GlobalDimensionCumulOptimizer* lp_optimizer,
     GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost);
 
 /// Returns a filter checking the feasibility and cost of the resource
diff --git a/ortools/routing/ils.cc b/ortools/routing/ils.cc
index d9463fa6d9..8af8cbb820 100644
--- a/ortools/routing/ils.cc
+++ b/ortools/routing/ils.cc
@@ -93,6 +93,11 @@ std::unique_ptr<RuinProcedure> MakeRuinProcedure(
       return std::make_unique<RandomWalkRemovalRuinProcedure>(
           model, rnd, ruin.random_walk().num_removed_visits(),
           num_neighbors_for_route_selection);
+    case RuinStrategy::kSisr:
+      return std::make_unique<SISRRuinProcedure>(
+          model, rnd, ruin.sisr().max_removed_sequence_size(),
+          ruin.sisr().avg_num_removed_visits(), ruin.sisr().bypass_factor(),
+          num_neighbors_for_route_selection);
     default:
       LOG(DFATAL) << "Unsupported ruin procedure.";
       return nullptr;
@@ -895,9 +900,14 @@ int64_t RandomWalkRemovalRuinProcedure::GetNextNodeToRemove(
 }
 
 SISRRuinProcedure::SISRRuinProcedure(RoutingModel* model, std::mt19937* rnd,
-                                     int num_neighbors)
+                                     int max_removed_sequence_size,
+                                     int avg_num_removed_visits,
+                                     double bypass_factor, int num_neighbors)
     : model_(*model),
       rnd_(*rnd),
+      max_removed_sequence_size_(max_removed_sequence_size),
+      avg_num_removed_visits_(avg_num_removed_visits),
+      bypass_factor_(bypass_factor),
       neighbors_manager_(model->GetOrCreateNodeNeighborsByCostClass(
           {num_neighbors,
            /*add_vehicle_starts_to_neighbors=*/false,
@@ -921,23 +931,17 @@ std::function<int64_t(int64_t)> SISRRuinProcedure::Ruin(
   routing_solution_.Reset(assignment);
   ruined_routes_.SparseClearAll();
 
-  // TODO(user): add to proto.
-  const int max_cardinality_removed_sequences = 10;
-
-  // TODO(user): add to proto.
-  const int avg_num_removed_visits = 10;
-
   const double max_sequence_size =
-      std::min<double>(max_cardinality_removed_sequences,
+      std::min<double>(max_removed_sequence_size_,
                        ComputeAverageNonEmptyRouteSize(model_, *assignment));
 
   const double max_num_removed_sequences =
-      (4 * avg_num_removed_visits) / (1 + max_sequence_size) - 1;
+      (4 * avg_num_removed_visits_) / (1 + max_sequence_size) - 1;
   DCHECK_GE(max_num_removed_sequences, 1);
 
   const int num_sequences_to_remove =
       std::floor(std::uniform_real_distribution<double>(
-          1.0, max_num_removed_sequences)(rnd_));
+          1.0, max_num_removed_sequences + 1)(rnd_));
 
   // We start by disrupting the route where the seed visit is served.
   const int seed_route = RuinRoute(*assignment, seed_node, max_sequence_size);
@@ -979,7 +983,7 @@ int SISRRuinProcedure::RuinRoute(const Assignment& assignment,
       routing_solution_.GetRouteSize(route), global_max_sequence_size);
 
   int sequence_size = std::floor(
-      std::uniform_real_distribution<double>(1.0, max_sequence_size)(rnd_));
+      std::uniform_real_distribution<double>(1.0, max_sequence_size + 1)(rnd_));
 
   if (sequence_size == 1 || sequence_size == max_sequence_size ||
       boolean_dist_(rnd_)) {
@@ -1011,14 +1015,11 @@ void SISRRuinProcedure::RuinRouteWithSequenceProcedure(int64_t seed_visit,
 void SISRRuinProcedure::RuinRouteWithSplitSequenceProcedure(int64_t route,
                                                             int64_t seed_visit,
                                                             int sequence_size) {
-  // TODO(user): add to proto.
-  const double alpha = 0.01;
-
   const int max_num_bypassed_visits =
       routing_solution_.GetRouteSize(route) - sequence_size;
   int num_bypassed_visits = 1;
   while (num_bypassed_visits < max_num_bypassed_visits &&
-         probability_dist_(rnd_) >= alpha * probability_dist_(rnd_)) {
+         probability_dist_(rnd_) >= bypass_factor_ * probability_dist_(rnd_)) {
     ++num_bypassed_visits;
   }
 
diff --git a/ortools/routing/ils.h b/ortools/routing/ils.h
index 4b9b35568c..cd8dd82472 100644
--- a/ortools/routing/ils.h
+++ b/ortools/routing/ils.h
@@ -206,9 +206,12 @@ class CompositeRuinProcedure : public RuinProcedure {
 // combination of user-defined parameters and solution and instance properties.
 // Every selected route is then disrupted by removing a contiguous sequence of
 // visits, possibly bypassing a contiguous subsequence.
+// See also SISRRuinStrategy in ils.proto.
 class SISRRuinProcedure : public RuinProcedure {
  public:
-  SISRRuinProcedure(RoutingModel* model, std::mt19937* rnd, int num_neighbors);
+  SISRRuinProcedure(RoutingModel* model, std::mt19937* rnd,
+                    int max_removed_sequence_size, int avg_num_removed_visits,
+                    double bypass_factor, int num_neighbors);
 
   std::function<int64_t(int64_t)> Ruin(const Assignment* assignment) override;
 
@@ -226,6 +229,9 @@ class SISRRuinProcedure : public RuinProcedure {
 
   const RoutingModel& model_;
   std::mt19937& rnd_;
+  int max_removed_sequence_size_;
+  int avg_num_removed_visits_;
+  double bypass_factor_;
   const RoutingModel::NodeNeighborsByCostClass* const neighbors_manager_;
   std::uniform_int_distribution<int64_t> customer_dist_;
   std::bernoulli_distribution boolean_dist_;
diff --git a/ortools/routing/ils.proto b/ortools/routing/ils.proto
index 8df144bedb..77e0328f1b 100644
--- a/ortools/routing/ils.proto
+++ b/ortools/routing/ils.proto
@@ -50,11 +50,99 @@ message RandomWalkRuinStrategy {
   optional uint32 num_removed_visits = 7;
 }
 
+// Ruin strategy based on the "Slack Induction by String Removals for Vehicle
+// Routing Problems" by Jan Christiaens and Greet Vanden Berghe, Transportation
+// Science 2020.
+// Link to paper:
+// https://kuleuven.limo.libis.be/discovery/fulldisplay?docid=lirias1988666&context=SearchWebhook&vid=32KUL_KUL:Lirias&lang=en&search_scope=lirias_profile&adaptor=SearchWebhook&tab=LIRIAS&query=any,contains,LIRIAS1988666&offset=0
+//
+// Note that, in this implementation, the notion of "string" is replaced by
+// "sequence".
+//
+// The main idea of this ruin is to remove a number of geographically close
+// sequences of nodes. In particular, at every ruin application, a maximum
+// number max_ruined_routes of routes are disrupted. The value for
+// max_ruined_routes is defined as
+//       (4 * avg_num_removed_visits) / (1 + max_sequence_size) + 1
+// with
+// - avg_num_removed_visits: user-defined parameter ruling the average number of
+//   visits that are removed in face of several ruin applications (see also the
+//   proto message below)
+// - max_sequence_size is defined as
+//        min{max_removed_sequence_size, average_route_size}
+//   with
+//   - max_removed_sequence_size: user-defined parameter that specifies
+//     the maximum number of visits removed from a single route (see also the
+//     proto message below)
+//   - average_route_size: the average size of a non-empty route in the current
+//     solution
+//
+// The actual number of ruined routes is then obtained as
+//                  floor(U(1, max_ruined_routes + 1))
+// where U is a continuous uniform distribution of real values in the given
+// interval.
+//
+// The routes affected by the ruin procedure are selected as follows.
+// First, a non start/end seed node is randomly selected. The route serving this
+// node is the first ruined route. Then, until the required number of routes has
+// been ruined, neighbor nodes of the initial seed node are scanned and the
+// associated not yet ruined routes are disrupted. Nodes defining the selected
+// routes are designated as seed nodes for the "sequence" and "split sequence"
+// removal procedures described below.
+//
+// For every selected route, a maximum number route_max_sequence_size of nodes
+// are removed. In particular, route_max_sequence_size is defined as
+//                min{route_size, max_sequence_size}
+// with route_size being the size of the current route.
+//
+// Then, the actual number of removed nodes num_removed_nodes is defined as
+//              floor(U(1, route_max_sequence_size + 1))
+// where U is a continuous uniform distribution of real values in the given
+// interval.
+//
+// As mentioned above, the selected num_removed_nodes number of nodes is removed
+// either via the "sequence" removal or "split sequence" removal procedures. The
+// two removal procedures are executed with equal probabilities.
+//
+// The "sequence" removal procedure removes a randomly selected sequence of size
+// num_removed_nodes that includes the seed node.
+//
+// The "split sequence" removal procedure also removes a randomly selected
+// sequence of size num_removed_nodes that includes the seed node, but it can
+// possibly preserve a subsequence of contiguous nodes.
+// In particular, the procedure first selects a sequence of size
+// num_removed_nodes + num_bypassed, then num_bypassed contiguous nodes in the
+// selected sequence are preserved while the others removed.
+//
+// The definition of num_bypassed is as follows. First num_bypassed = 1. The
+// current value of num_bypassed is maintaned if
+//          U(0, 1) < bypass_factor * U(0, 1)
+// or the maximum value for num_bypassed, equal to
+//           route_size - num_removed_nodes
+// is reached. The value is incremented of a unit otherwise,
+// and the process is repeated. The value assigned to bypass_factor affects the
+// number of preserved visits (see also the proto message below).
+message SISRRuinStrategy {
+  // Maximum number of removed visits per sequence. The parameter name in the
+  // paper is L^{max} and the suggested value is 10.
+  optional uint32 max_removed_sequence_size = 1;
+
+  // Number of visits that are removed on average. The parameter name in the
+  // paper is \bar{c} and the suggested value is 10.
+  optional uint32 avg_num_removed_visits = 2;
+
+  // Value in [0, 1] ruling the number of preserved customers in the split
+  // sequence removal. The parameter name in the paper is \alpha and the
+  // suggested value is 0.01.
+  optional double bypass_factor = 3;
+}
+
 // Ruin strategies, used in perturbation based on ruin and recreate approaches.
 message RuinStrategy {
   oneof strategy {
     SpatiallyCloseRoutesRuinStrategy spatially_close_routes = 1;
     RandomWalkRuinStrategy random_walk = 2;
+    SISRRuinStrategy sisr = 3;
   }
 }
 
diff --git a/ortools/routing/lp_scheduling.cc b/ortools/routing/lp_scheduling.cc
index 6f37bfebb0..4a0482dc53 100644
--- a/ortools/routing/lp_scheduling.cc
+++ b/ortools/routing/lp_scheduling.cc
@@ -17,11 +17,8 @@
 #include <cmath>
 #include <cstdint>
 #include <cstdlib>
-#include <deque>
 #include <functional>
-#include <iterator>
 #include <limits>
-#include <memory>
 #include <numeric>
 #include <optional>
 #include <string>
@@ -29,8 +26,6 @@
 #include <vector>
 
 #include "absl/algorithm/container.h"
-#include "absl/container/flat_hash_map.h"
-#include "absl/container/flat_hash_set.h"
 #include "absl/log/check.h"
 #include "absl/strings/str_format.h"
 #include "absl/strings/string_view.h"
@@ -40,7 +35,6 @@
 #include "ortools/base/logging.h"
 #include "ortools/base/map_util.h"
 #include "ortools/base/mathutil.h"
-#include "ortools/base/strong_vector.h"
 #include "ortools/base/types.h"
 #include "ortools/constraint_solver/constraint_solver.h"
 #include "ortools/glop/parameters.pb.h"
@@ -52,6 +46,7 @@
 #include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/lp_utils.h"
 #include "ortools/util/flat_matrix.h"
+#include "ortools/util/piecewise_linear_function.h"
 #include "ortools/util/saturated_arithmetic.h"
 #include "ortools/util/sorted_interval_list.h"
 
@@ -219,7 +214,7 @@ DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumulCost(
   const DimensionSchedulingStatus status =
       optimizer_core_.OptimizeSingleRouteWithResource(
           vehicle, next_accessor,
-          /*dimension_travel_info=*/{},
+          /*dimension_travel_info=*/nullptr,
           /*resource=*/nullptr,
           /*optimize_vehicle_costs=*/optimal_cost != nullptr,
           solver_[vehicle].get(), /*cumul_values=*/nullptr,
@@ -238,7 +233,7 @@ LocalDimensionCumulOptimizer::ComputeRouteCumulCostWithoutFixedTransits(
     const RoutingModel::ResourceGroup::Resource* resource,
     int64_t* optimal_cost_without_transits) {
   return optimizer_core_.OptimizeSingleRouteWithResource(
-      vehicle, next_accessor, /*dimension_travel_info=*/{}, resource,
+      vehicle, next_accessor, /*dimension_travel_info=*/nullptr, resource,
       /*optimize_vehicle_costs=*/optimal_cost_without_transits != nullptr,
       solver_[vehicle].get(), /*cumul_values=*/nullptr,
       /*break_values=*/nullptr, optimal_cost_without_transits, nullptr);
@@ -254,14 +249,14 @@ std::vector<DimensionSchedulingStatus> LocalDimensionCumulOptimizer::
         std::vector<std::vector<int64_t>>* optimal_cumuls,
         std::vector<std::vector<int64_t>>* optimal_breaks) {
   return optimizer_core_.OptimizeSingleRouteWithResources(
-      vehicle, next_accessor, transit_accessor, {}, resources, resource_indices,
-      optimize_vehicle_costs, solver_[vehicle].get(), optimal_cumuls,
-      optimal_breaks, optimal_costs_without_transits, nullptr);
+      vehicle, next_accessor, transit_accessor, nullptr, resources,
+      resource_indices, optimize_vehicle_costs, solver_[vehicle].get(),
+      optimal_cumuls, optimal_breaks, optimal_costs_without_transits, nullptr);
 }
 
 DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumuls(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+    const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
     const RoutingModel::ResourceGroup::Resource* resource,
     std::vector<int64_t>* optimal_cumuls,
     std::vector<int64_t>* optimal_breaks) {
@@ -275,7 +270,7 @@ DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumuls(
 DimensionSchedulingStatus
 LocalDimensionCumulOptimizer::ComputeRouteCumulsAndCostWithoutFixedTransits(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+    const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
     std::vector<int64_t>* optimal_cumuls, std::vector<int64_t>* optimal_breaks,
     int64_t* optimal_cost_without_transits) {
   return optimizer_core_.OptimizeSingleRouteWithResource(
@@ -287,7 +282,7 @@ LocalDimensionCumulOptimizer::ComputeRouteCumulsAndCostWithoutFixedTransits(
 DimensionSchedulingStatus
 LocalDimensionCumulOptimizer::ComputeRouteSolutionCostWithoutFixedTransits(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+    const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
     absl::Span<const int64_t> solution_cumul_values,
     absl::Span<const int64_t> solution_break_values, int64_t* solution_cost,
     int64_t* cost_offset, bool reuse_previous_model_if_possible, bool clear_lp,
@@ -303,7 +298,7 @@ LocalDimensionCumulOptimizer::ComputeRouteSolutionCostWithoutFixedTransits(
 DimensionSchedulingStatus
 LocalDimensionCumulOptimizer::ComputePackedRouteCumuls(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+    const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
     const RoutingModel::ResourceGroup::Resource* resource,
     std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks) {
   return optimizer_core_.OptimizeAndPackSingleRoute(
@@ -580,7 +575,7 @@ DimensionCumulOptimizerCore::DimensionCumulOptimizerCore(
 DimensionSchedulingStatus
 DimensionCumulOptimizerCore::ComputeSingleRouteSolutionCostWithoutFixedTransits(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RouteDimensionTravelInfo& dimension_travel_info,
+    const RouteDimensionTravelInfo* dimension_travel_info,
     RoutingLinearSolverWrapper* solver,
     absl::Span<const int64_t> solution_cumul_values,
     absl::Span<const int64_t> solution_break_values,
@@ -691,7 +686,7 @@ void ClearIfNonNull(std::vector<T>* v) {
 DimensionSchedulingStatus
 DimensionCumulOptimizerCore::OptimizeSingleRouteWithResource(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RouteDimensionTravelInfo& dimension_travel_info,
+    const RouteDimensionTravelInfo* dimension_travel_info,
     const RoutingModel::ResourceGroup::Resource* resource,
     bool optimize_vehicle_costs, RoutingLinearSolverWrapper* solver,
     std::vector<int64_t>* cumul_values, std::vector<int64_t>* break_values,
@@ -830,7 +825,7 @@ std::vector<DimensionSchedulingStatus>
 DimensionCumulOptimizerCore::OptimizeSingleRouteWithResources(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
     const std::function<int64_t(int64_t, int64_t)>& transit_accessor,
-    const RouteDimensionTravelInfo& dimension_travel_info,
+    const RouteDimensionTravelInfo* dimension_travel_info,
     absl::Span<const RoutingModel::ResourceGroup::Resource> resources,
     absl::Span<const int> resource_indices, bool optimize_vehicle_costs,
     RoutingLinearSolverWrapper* solver,
@@ -969,10 +964,10 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize(
         !model->IsEnd(next_accessor(model->Start(vehicle))) ||
         model->IsVehicleUsedWhenEmpty(vehicle);
     const bool optimize_vehicle_costs = optimize_costs && vehicle_is_used;
-    const RouteDimensionTravelInfo& dimension_travel_info =
+    const RouteDimensionTravelInfo* const dimension_travel_info =
         dimension_travel_info_per_route.empty()
-            ? RouteDimensionTravelInfo()
-            : dimension_travel_info_per_route[vehicle];
+            ? nullptr
+            : &dimension_travel_info_per_route[vehicle];
     if (!SetRouteCumulConstraints(
             vehicle, next_accessor, dimension_->transit_evaluator(vehicle),
             dimension_travel_info, cumul_offset, optimize_vehicle_costs, solver,
@@ -1078,7 +1073,7 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack(
 DimensionSchedulingStatus
 DimensionCumulOptimizerCore::OptimizeAndPackSingleRoute(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-    const RouteDimensionTravelInfo& dimension_travel_info,
+    const RouteDimensionTravelInfo* dimension_travel_info,
     const RoutingModel::ResourceGroup::Resource* resource,
     RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
     std::vector<int64_t>* break_values) {
@@ -1271,6 +1266,24 @@ bool DimensionCumulOptimizerCore::TightenRouteCumulBounds(
   return true;
 }
 
+std::vector<SlopeAndYIntercept> PiecewiseLinearFunctionToSlopeAndYIntercept(
+    const FloatSlopePiecewiseLinearFunction& pwl_function, int index_start,
+    int index_end) {
+  const auto& x_anchors = pwl_function.x_anchors();
+  const auto& y_anchors = pwl_function.y_anchors();
+  if (index_end < 0) index_end = x_anchors.size() - 1;
+  const int num_segments = index_end - index_start;
+  DCHECK_GE(num_segments, 1);
+  std::vector<SlopeAndYIntercept> slope_and_y_intercept(num_segments);
+  for (int seg = index_start; seg < index_end; ++seg) {
+    auto& [slope, y_intercept] = slope_and_y_intercept[seg - index_start];
+    slope = (y_anchors[seg + 1] - y_anchors[seg]) /
+            static_cast<double>(x_anchors[seg + 1] - x_anchors[seg]);
+    y_intercept = y_anchors[seg] - slope * x_anchors[seg];
+  }
+  return slope_and_y_intercept;
+}
+
 std::vector<bool> SlopeAndYInterceptToConvexityRegions(
     absl::Span<const SlopeAndYIntercept> slope_and_y_intercept) {
   CHECK(!slope_and_y_intercept.empty());
@@ -1286,27 +1299,7 @@ std::vector<bool> SlopeAndYInterceptToConvexityRegions(
   return convex;
 }
 
-std::vector<SlopeAndYIntercept> PiecewiseLinearFormulationToSlopeAndYIntercept(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl_function,
-    int index_start, int index_end) {
-  if (index_end < 0) index_end = pwl_function.x_anchors.size() - 1;
-  const int num_segments = index_end - index_start;
-  DCHECK_GE(num_segments, 1);
-  std::vector<SlopeAndYIntercept> slope_and_y_intercept(num_segments);
-  for (int seg = index_start; seg < index_end; ++seg) {
-    auto& [slope, y_intercept] = slope_and_y_intercept[seg - index_start];
-    slope = (pwl_function.y_anchors[seg + 1] - pwl_function.y_anchors[seg]) /
-            static_cast<double>(pwl_function.x_anchors[seg + 1] -
-                                pwl_function.x_anchors[seg]);
-    y_intercept =
-        pwl_function.y_anchors[seg] - slope * pwl_function.x_anchors[seg];
-  }
-  return slope_and_y_intercept;
-}
-
 namespace {
-
 // Find a "good" scaling factor for constraints with non-integers coefficients.
 // See sat::FindBestScalingAndComputeErrors() for more infos.
 double FindBestScaling(const std::vector<double>& coefficients,
@@ -1321,88 +1314,17 @@ double FindBestScaling(const std::vector<double>& coefficients,
       wanted_absolute_activity_precision, &unused_relative_coeff_error,
       &unused_scaled_sum_error);
 }
-
-// Returns the value of pwl(x) with pwl a PiecewiseLinearFormulation, knowing
-// that x ∈ [pwl.x[upper_segment_index-1], pwl.x[upper_segment_index]].
-int64_t PieceWiseLinearFormulationValueKnownSegment(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl,
-    int64_t x, int upper_segment_index, double delta = 0) {
-  DCHECK_GE(upper_segment_index, 1);
-  DCHECK_LE(upper_segment_index, pwl.x_anchors.size() - 1);
-  const double alpha =
-      static_cast<double>(pwl.y_anchors[upper_segment_index] -
-                          pwl.y_anchors[upper_segment_index - 1]) /
-      (pwl.x_anchors[upper_segment_index] -
-       pwl.x_anchors[upper_segment_index - 1]);
-  const double beta = pwl.y_anchors[upper_segment_index] -
-                      pwl.x_anchors[upper_segment_index] * alpha;
-  return std::ceil(alpha * x + beta + delta);
-}
-
 }  // namespace
 
-PiecewiseEvaluationStatus ComputePiecewiseLinearFormulationValue(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl,
-    int64_t x, int64_t* value, double delta) {
-  // Search for first element xi such that xi < x.
-  const auto upper_segment =
-      std::upper_bound(pwl.x_anchors.begin(), pwl.x_anchors.end(), x);
-  const int upper_segment_index =
-      std::distance(pwl.x_anchors.begin(), upper_segment);
-
-  // Checking bounds
-  if (upper_segment_index == 0) {
-    return PiecewiseEvaluationStatus::SMALLER_THAN_LOWER_BOUND;
-  } else if (upper_segment == pwl.x_anchors.end()) {
-    if (x == pwl.x_anchors.back()) {
-      *value = std::ceil(pwl.y_anchors.back() + delta);
-      return PiecewiseEvaluationStatus::WITHIN_BOUNDS;
-    }
-    return PiecewiseEvaluationStatus::LARGER_THAN_UPPER_BOUND;
-  }
-
-  *value = PieceWiseLinearFormulationValueKnownSegment(
-      pwl, x, upper_segment_index, delta);
-  return PiecewiseEvaluationStatus::WITHIN_BOUNDS;
-}
-
-int64_t ComputeConvexPiecewiseLinearFormulationValue(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl,
-    int64_t x, double delta) {
-  int64_t y_value;
-  const PiecewiseEvaluationStatus status =
-      ComputePiecewiseLinearFormulationValue(pwl, x, &y_value, delta);
-  switch (status) {
-    case PiecewiseEvaluationStatus::UNSPECIFIED:
-      // The status should be specified.
-      LOG(FATAL) << "Unspecified PiecewiseEvaluationStatus.";
-      break;
-    case PiecewiseEvaluationStatus::WITHIN_BOUNDS:
-      // x is in the bounds, therefore, simply return the computed value.
-      return y_value;
-    case PiecewiseEvaluationStatus::SMALLER_THAN_LOWER_BOUND:
-      // In the convex case, if x <= lower_bound, the most restrictive
-      // constraint will be the first one.
-      return PieceWiseLinearFormulationValueKnownSegment(pwl, x, 1, delta);
-    case PiecewiseEvaluationStatus::LARGER_THAN_UPPER_BOUND:
-      // In the convex case, if x >= upper_bound, the most restrictive
-      // constraint will be the last one.
-      return PieceWiseLinearFormulationValueKnownSegment(
-          pwl, x, pwl.x_anchors.size() - 1, delta);
-  }
-}
-
 bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
-    const RouteDimensionTravelInfo& dimension_travel_info,
+    const RouteDimensionTravelInfo* dimension_travel_info,
     absl::Span<const int> lp_slacks, absl::Span<const int64_t> fixed_transit,
     RoutingLinearSolverWrapper* solver) {
   const std::vector<int>& lp_cumuls = current_route_cumul_variables_;
   const int path_size = lp_cumuls.size();
 
-  if (dimension_travel_info.transition_info.empty()) {
+  if (dimension_travel_info == nullptr ||
+      dimension_travel_info->transition_info.empty()) {
     // Travel is not travel-start dependent.
     // Add all path constraints to LP:
     // cumul[i] + fixed_transit[i] + slack[i] == cumul[i+1]
@@ -1438,14 +1360,10 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
         absl::StrFormat("relative_compression_cost(%ld)", pos));
 
     const RoutingModel::RouteDimensionTravelInfo::TransitionInfo&
-        transition_info = dimension_travel_info.transition_info[pos];
-    const RouteDimensionTravelInfo::TransitionInfo::PiecewiseLinearFormulation&
-        travel_function = transition_info.travel_start_dependent_travel;
-    const int num_pwl_anchors = travel_function.x_anchors.size();
-    DCHECK_GE(num_pwl_anchors, 2)
-        << "Travel value PWL must have at least 2 points";
-    DCHECK_EQ(num_pwl_anchors, travel_function.y_anchors.size())
-        << "Travel value PWL must have as many x anchors than y.";
+        transition_info = dimension_travel_info->transition_info[pos];
+    const FloatSlopePiecewiseLinearFunction& travel_function =
+        transition_info.travel_start_dependent_travel;
+    const auto& travel_x_anchors = travel_function.x_anchors();
 
     // 1. Create the travel value variable and set its constraints.
     // 1.a. Create Variables for the start and value of a travel
@@ -1455,7 +1373,7 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
     const int64_t compressed_travel_value_lower_bound =
         transition_info.compressed_travel_value_lower_bound;
     const int64_t travel_value_upper_bound =
-        dimension_travel_info.transition_info[pos].travel_value_upper_bound;
+        transition_info.travel_value_upper_bound;
     // The lower bound of travel_value is already implemented by constraints as
     // travel_value >= compressed_travel_value (defined below) and
     // compressed_travel_value has compressed_travel_value_lower_bound as a
@@ -1479,15 +1397,16 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
     // Find segments that are in bounds
     // Only the segments in [index_anchor_start, index_anchor_end[ are in
     // bounds, the others can therefore be discarded.
+    const int num_pwl_anchors = travel_x_anchors.size();
     int index_anchor_start = 0;
     while (index_anchor_start < num_pwl_anchors - 1 &&
-           travel_function.x_anchors[index_anchor_start + 1] <=
+           travel_x_anchors[index_anchor_start + 1] <=
                current_route_min_cumuls_[pos] + pre_travel_transit) {
       ++index_anchor_start;
     }
     int index_anchor_end = num_pwl_anchors - 1;
     while (index_anchor_end > 0 &&
-           travel_function.x_anchors[index_anchor_end - 1] >=
+           travel_x_anchors[index_anchor_end - 1] >=
                current_route_max_cumuls_[pos] + pre_travel_transit) {
       --index_anchor_end;
     }
@@ -1497,7 +1416,7 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
     // Precompute the slopes and y-intercept as they will be used to detect
     // convexities and in the constraints.
     const std::vector<SlopeAndYIntercept> slope_and_y_intercept =
-        PiecewiseLinearFormulationToSlopeAndYIntercept(
+        PiecewiseLinearFunctionToSlopeAndYIntercept(
             travel_function, index_anchor_start, index_anchor_end);
 
     // Optimize binary variables by detecting convexities
@@ -1531,7 +1450,7 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
         // If the travel_start value is outside the PWL, the closest segment
         // will be used. This is why some bounds are infinite.
         const int64_t lower_bound_interval =
-            seg > 0 ? travel_function.x_anchors[index_anchor_start + seg]
+            seg > 0 ? travel_x_anchors[index_anchor_start + seg]
                     : current_route_min_cumuls_[pos] + pre_travel_transit;
         int64_t end_of_seg = seg + 1;
         while (end_of_seg < num_pwl_anchors - 1 && !convexities[end_of_seg]) {
@@ -1539,7 +1458,7 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
         }
         const int64_t higher_bound_interval =
             end_of_seg < num_pwl_anchors - 1
-                ? travel_function.x_anchors[index_anchor_start + end_of_seg]
+                ? travel_x_anchors[index_anchor_start + end_of_seg]
                 : current_route_max_cumuls_[pos] + pre_travel_transit;
         const int travel_start_in_segment_ct = solver->AddLinearConstraint(
             lower_bound_interval, higher_bound_interval, {{travel_start, 1}});
@@ -1640,25 +1559,25 @@ bool DimensionCumulOptimizerCore::SetRouteTravelConstraints(
     // compressed_travel_value to not give the incentive to compress a little
     // bit in order to same some cost per travel.
     solver->SetObjectiveCoefficient(
-        travel_value, dimension_travel_info.travel_cost_coefficient);
+        travel_value, dimension_travel_info->travel_cost_coefficient);
 
     // 4. Adds a convex cost in epsilon
     // Here we DCHECK that the cost function is indeed convex
-    const RouteDimensionTravelInfo::TransitionInfo::PiecewiseLinearFormulation&
-        cost_function =
-            dimension_travel_info.transition_info[pos].travel_compression_cost;
+    const FloatSlopePiecewiseLinearFunction& cost_function =
+        transition_info.travel_compression_cost;
+    const auto& cost_x_anchors = cost_function.x_anchors();
+
     const std::vector<SlopeAndYIntercept> cost_slope_and_y_intercept =
-        PiecewiseLinearFormulationToSlopeAndYIntercept(cost_function);
-    const double cost_max = ComputeConvexPiecewiseLinearFormulationValue(
-        cost_function,
+        PiecewiseLinearFunctionToSlopeAndYIntercept(cost_function);
+    const double cost_max = cost_function.ComputeConvexValue(
         travel_value_upper_bound - compressed_travel_value_lower_bound);
     double previous_slope = 0;
-    for (int seg = 0; seg < cost_function.x_anchors.size() - 1; ++seg) {
+    for (int seg = 0; seg < cost_x_anchors.size() - 1; ++seg) {
       const auto [slope, y_intercept] = cost_slope_and_y_intercept[seg];
       // Check convexity
       DCHECK_GE(slope, previous_slope)
           << "Compression error is not convex. Segment " << (1 + seg)
-          << " out of " << (cost_function.x_anchors.size() - 1);
+          << " out of " << (cost_x_anchors.size() - 1);
       previous_slope = slope;
       const double factor = FindBestScaling(
           {1.0, -slope, y_intercept}, /*lower_bounds=*/
@@ -1723,7 +1642,7 @@ bool RouteIsValid(const RoutingModel& model, int vehicle,
 bool DimensionCumulOptimizerCore::SetRouteCumulConstraints(
     int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
     const std::function<int64_t(int64_t, int64_t)>& transit_accessor,
-    const RouteDimensionTravelInfo& dimension_travel_info, int64_t cumul_offset,
+    const RouteDimensionTravelInfo* dimension_travel_info, int64_t cumul_offset,
     bool optimize_costs, RoutingLinearSolverWrapper* solver,
     int64_t* route_transit_cost, int64_t* route_cost_offset) {
   RoutingModel* const model = dimension_->model();
@@ -1753,7 +1672,8 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints(
   if (!ExtractRouteCumulBounds(path, cumul_offset)) {
     return false;
   }
-  if (dimension_travel_info.transition_info.empty()) {
+  if (dimension_travel_info == nullptr ||
+      dimension_travel_info->transition_info.empty()) {
     if (!TightenRouteCumulBounds(path, fixed_transit, cumul_offset)) {
       return false;
     }
@@ -1762,7 +1682,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints(
     std::vector<int64_t> min_transit(path_size - 1);
     for (int pos = 0; pos < path_size - 1; ++pos) {
       const RouteDimensionTravelInfo::TransitionInfo& transition =
-          dimension_travel_info.transition_info[pos];
+          dimension_travel_info->transition_info[pos];
       min_transit[pos] = transition.pre_travel_transit_value +
                          transition.compressed_travel_value_lower_bound +
                          transition.post_travel_transit_value;
diff --git a/ortools/routing/lp_scheduling.h b/ortools/routing/lp_scheduling.h
index d2706ff995..7176298e08 100644
--- a/ortools/routing/lp_scheduling.h
+++ b/ortools/routing/lp_scheduling.h
@@ -45,6 +45,7 @@
 #include "ortools/sat/cp_model_solver.h"
 #include "ortools/sat/model.h"
 #include "ortools/sat/sat_parameters.pb.h"
+#include "ortools/util/piecewise_linear_function.h"
 #include "ortools/util/sorted_interval_list.h"
 
 namespace operations_research::routing {
@@ -648,7 +649,7 @@ class DimensionCumulOptimizerCore {
   // resource. If the resource is null, it is simply ignored.
   DimensionSchedulingStatus OptimizeSingleRouteWithResource(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       const Resource* resource, bool optimize_vehicle_costs,
       RoutingLinearSolverWrapper* solver, std::vector<int64_t>* cumul_values,
       std::vector<int64_t>* break_values, int64_t* cost_without_transit,
@@ -659,7 +660,7 @@ class DimensionCumulOptimizerCore {
   // constraints for cumuls and breaks.
   DimensionSchedulingStatus ComputeSingleRouteSolutionCostWithoutFixedTransits(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       RoutingLinearSolverWrapper* solver,
       absl::Span<const int64_t> solution_cumul_values,
       absl::Span<const int64_t> solution_break_values,
@@ -675,7 +676,7 @@ class DimensionCumulOptimizerCore {
   std::vector<DimensionSchedulingStatus> OptimizeSingleRouteWithResources(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
       const std::function<int64_t(int64_t, int64_t)>& transit_accessor,
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       absl::Span<const Resource> resources,
       absl::Span<const int> resource_indices, bool optimize_vehicle_costs,
       RoutingLinearSolverWrapper* solver,
@@ -708,7 +709,7 @@ class DimensionCumulOptimizerCore {
 
   DimensionSchedulingStatus OptimizeAndPackSingleRoute(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       const Resource* resource, RoutingLinearSolverWrapper* solver,
       std::vector<int64_t>* cumul_values, std::vector<int64_t>* break_values);
 
@@ -738,7 +739,7 @@ class DimensionCumulOptimizerCore {
   bool SetRouteCumulConstraints(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
       const std::function<int64_t(int64_t, int64_t)>& transit_accessor,
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       int64_t cumul_offset, bool optimize_costs,
       RoutingLinearSolverWrapper* solver, int64_t* route_transit_cost,
       int64_t* route_cost_offset);
@@ -747,7 +748,7 @@ class DimensionCumulOptimizerCore {
   // static or time-dependent travel values.
   // Returns false if some infeasibility was detected, true otherwise.
   bool SetRouteTravelConstraints(
-      const RouteDimensionTravelInfo& dimension_travel_info,
+      const RouteDimensionTravelInfo* dimension_travel_info,
       absl::Span<const int> lp_slacks, absl::Span<const int64_t> fixed_transit,
       RoutingLinearSolverWrapper* solver);
 
@@ -865,7 +866,7 @@ class LocalDimensionCumulOptimizer {
   // Returns false if the route is not feasible.
   DimensionSchedulingStatus ComputeRouteCumuls(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+      const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
       const RoutingModel::ResourceGroup::Resource* resource,
       std::vector<int64_t>* optimal_cumuls,
       std::vector<int64_t>* optimal_breaks);
@@ -874,7 +875,7 @@ class LocalDimensionCumulOptimizer {
   // ComputeRouteCumuls().
   DimensionSchedulingStatus ComputeRouteCumulsAndCostWithoutFixedTransits(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+      const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
       std::vector<int64_t>* optimal_cumuls,
       std::vector<int64_t>* optimal_breaks,
       int64_t* optimal_cost_without_transits);
@@ -883,7 +884,7 @@ class LocalDimensionCumulOptimizer {
   // defined by its cumuls and breaks.
   DimensionSchedulingStatus ComputeRouteSolutionCostWithoutFixedTransits(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+      const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
       absl::Span<const int64_t> solution_cumul_values,
       absl::Span<const int64_t> solution_break_values, int64_t* solution_cost,
       int64_t* cost_offset = nullptr,
@@ -897,7 +898,7 @@ class LocalDimensionCumulOptimizer {
   // time window.
   DimensionSchedulingStatus ComputePackedRouteCumuls(
       int vehicle, const std::function<int64_t(int64_t)>& next_accessor,
-      const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info,
+      const RoutingModel::RouteDimensionTravelInfo* dimension_travel_info,
       const RoutingModel::ResourceGroup::Resource* resource,
       std::vector<int64_t>* packed_cumuls, std::vector<int64_t>* packed_breaks);
 
@@ -1009,35 +1010,8 @@ bool ComputeVehicleToResourceClassAssignmentCosts(
     std::vector<std::vector<int64_t>>* cumul_values,
     std::vector<std::vector<int64_t>>* break_values);
 
-// Simple struct returned by ComputePiecewiseLinearFormulationValue() to
-// indicate if the value could be computed and if not, on what side the value
-// was from the definition interval.
-enum class PiecewiseEvaluationStatus {
-  UNSPECIFIED = 0,
-  WITHIN_BOUNDS,
-  SMALLER_THAN_LOWER_BOUND,
-  LARGER_THAN_UPPER_BOUND
-};
-
-// Computes pwl(x) for pwl a PieceWiseLinearFormulation.
-// Returns a PieceWiseEvaluationStatus to indicate if the value could be
-// computed (filled in value) and if not, why.
-PiecewiseEvaluationStatus ComputePiecewiseLinearFormulationValue(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl,
-    int64_t x, int64_t* value, double delta = 0);
-
-// Like ComputePiecewiseLinearFormulationValue(), computes pwl(x) for pwl a
-// PiecewiseLinearFormulation. For convex PiecewiseLinearFormulations, if x is
-// outside the bounds of the function, instead of returning an error like in
-// PiecewiseLinearFormulation, the function will still be defined by its outer
-// segments.
-int64_t ComputeConvexPiecewiseLinearFormulationValue(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl,
-    int64_t x, double delta = 0);
-
-// Structure to store the slope and y_intercept of a segment.
+// Structure to store the slope and y_intercept of a segment of a piecewise
+// linear function.
 struct SlopeAndYIntercept {
   double slope;
   double y_intercept;
@@ -1048,6 +1022,16 @@ struct SlopeAndYIntercept {
   }
 };
 
+// Given a FloatSlopePiecewiseLinearFunction, returns a vector of slope and
+// y-intercept corresponding to each segment. Only the segments in
+// [index_start, index_end[ will be considered.
+// TODO(user): Consider making the following two functions methods of
+// FloatSlopePiecewiseLinearFunction. They're only called in lp_scheduling.cc
+// and ../tour_optimization/model_test.cc, but they might come in handy.
+std::vector<SlopeAndYIntercept> PiecewiseLinearFunctionToSlopeAndYIntercept(
+    const FloatSlopePiecewiseLinearFunction& pwl_function, int index_start = 0,
+    int index_end = -1);
+
 // Converts a vector of SlopeAndYIntercept to a vector of convexity regions.
 // Convexity regions are defined such that, all segment in a convexity region
 // form a convex function. The boolean in the vector is set to true if the
@@ -1057,14 +1041,6 @@ struct SlopeAndYIntercept {
 std::vector<bool> SlopeAndYInterceptToConvexityRegions(
     absl::Span<const SlopeAndYIntercept> slope_and_y_intercept);
 
-// Given a PiecewiseLinearFormulation, returns a vector of slope and y-intercept
-// corresponding to each segment. Only the segments in [index_start, index_end[
-// will be considered.
-std::vector<SlopeAndYIntercept> PiecewiseLinearFormulationToSlopeAndYIntercept(
-    const RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-        PiecewiseLinearFormulation& pwl_function,
-    int index_start = 0, int index_end = -1);
-
 }  // namespace operations_research::routing
 
 #endif  // OR_TOOLS_ROUTING_LP_SCHEDULING_H_
diff --git a/ortools/routing/parameters.cc b/ortools/routing/parameters.cc
index a991cb07a9..57b93f3651 100644
--- a/ortools/routing/parameters.cc
+++ b/ortools/routing/parameters.cc
@@ -122,7 +122,7 @@ RoutingSearchParameters CreateDefaultRoutingSearchParameters() {
   o->set_use_make_inactive(BOOL_TRUE);
   o->set_use_make_chain_inactive(BOOL_TRUE);
   o->set_use_swap_active(BOOL_TRUE);
-  o->set_use_swap_active_chain(BOOL_FALSE);
+  o->set_use_swap_active_chain(BOOL_TRUE);
   o->set_use_extended_swap_active(BOOL_FALSE);
   o->set_use_shortest_path_swap_active(BOOL_TRUE);
   o->set_use_node_pair_swap_active(BOOL_FALSE);
@@ -143,6 +143,7 @@ RoutingSearchParameters CreateDefaultRoutingSearchParameters() {
   p.set_use_multi_armed_bandit_concatenate_operators(false);
   p.set_multi_armed_bandit_compound_operator_memory_coefficient(0.04);
   p.set_multi_armed_bandit_compound_operator_exploration_coefficient(1e12);
+  p.set_max_swap_active_chain_size(10);
   p.set_relocate_expensive_chain_num_arcs_to_consider(4);
   p.set_heuristic_expensive_chain_lns_num_arcs_to_consider(4);
   p.set_heuristic_close_nodes_lns_num_nodes(5);
@@ -317,6 +318,28 @@ void FindErrorsInIteratedLocalSearchParameters(
                    "ruin_strategy is set to RandomWalkRuinStrategy"
                    " but random_walk.num_removed_visits is 0 (should be "
                    "strictly positive)"));
+      } else if (ruin.strategy_case() == RuinStrategy::kSisr) {
+        if (ruin.sisr().avg_num_removed_visits() == 0) {
+          errors.emplace_back(
+              "iterated_local_search_parameters.ruin_recreate_parameters."
+              "ruin is set to SISRRuinStrategy"
+              " but sisr.avg_num_removed_visits is 0 (should be strictly "
+              "positive)");
+        }
+        if (ruin.sisr().max_removed_sequence_size() == 0) {
+          errors.emplace_back(
+              "iterated_local_search_parameters.ruin_recreate_parameters.ruin "
+              "is set to SISRRuinStrategy but "
+              "sisr.max_removed_sequence_size is 0 (should be strictly "
+              "positive)");
+        }
+        if (ruin.sisr().bypass_factor() < 0 ||
+            ruin.sisr().bypass_factor() > 1) {
+          errors.emplace_back(StrCat(
+              "iterated_local_search_parameters.ruin_recreate_parameters."
+              "ruin is set to SISRRuinStrategy"
+              " but sisr.bypass_factor is not in [0, 1]"));
+        }
       }
     }
 
@@ -687,6 +710,14 @@ std::vector<std::string> FindErrorsInRoutingSearchParameters(
         " search");
   }
 
+  if (search_parameters.max_swap_active_chain_size() < 1 &&
+      search_parameters.local_search_operators().use_swap_active_chain() ==
+          OptionalBoolean::BOOL_TRUE) {
+    errors.emplace_back(
+        "max_swap_active_chain_size must be greater than 1 if "
+        "local_search_operators.use_swap_active_chain is BOOL_TRUE");
+  }
+
   FindErrorsInIteratedLocalSearchParameters(search_parameters, errors);
 
   return errors;
diff --git a/ortools/routing/parameters.proto b/ortools/routing/parameters.proto
index 0863ccac4f..96a4fc9b17 100644
--- a/ortools/routing/parameters.proto
+++ b/ortools/routing/parameters.proto
@@ -36,7 +36,7 @@ package operations_research.routing;
 // then the routing library will pick its preferred value for that parameter
 // automatically: this should be the case for most parameters.
 // To see those "default" parameters, call GetDefaultRoutingSearchParameters().
-// Next ID: 66
+// Next ID: 67
 message RoutingSearchParameters {
   reserved 19;
 
@@ -459,6 +459,9 @@ message RoutingSearchParameters {
   // unsuccessful in the past operators
   double multi_armed_bandit_compound_operator_exploration_coefficient = 43;
 
+  // Maximum size of the chain to make inactive in SwapActiveChainOperator.
+  int32 max_swap_active_chain_size = 66;
+
   // Number of expensive arcs to consider cutting in the RelocateExpensiveChain
   // neighborhood operator (see
   // LocalSearchNeighborhoodOperators.use_relocate_expensive_chain()).
diff --git a/ortools/routing/routing.cc b/ortools/routing/routing.cc
index fd4098192a..bbe5d39088 100644
--- a/ortools/routing/routing.cc
+++ b/ortools/routing/routing.cc
@@ -50,7 +50,6 @@
 #include "absl/time/time.h"
 #include "absl/types/span.h"
 #include "google/protobuf/util/message_differencer.h"
-#include "ortools/base/dump_vars.h"
 #include "ortools/base/int_type.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/map_util.h"
@@ -138,16 +137,6 @@ std::string RoutingModel::RouteDimensionTravelInfo::TransitionInfo::DebugString(
       line_prefix, travel_compression_cost.DebugString(line_prefix + "\t"));
 }
 
-std::string RoutingModel::RouteDimensionTravelInfo::TransitionInfo::
-    PiecewiseLinearFormulation::DebugString(std::string line_prefix) const {
-  if (x_anchors.size() <= 10) {
-    return "{ " + DUMP_VARS(x_anchors, y_anchors).str() + "}";
-  }
-  return absl::StrFormat("{\n%s%s\n%s%s\n}", line_prefix,
-                         DUMP_VARS(x_anchors).str(), line_prefix,
-                         DUMP_VARS(y_anchors).str());
-}
-
 const Assignment* RoutingModel::PackCumulsOfOptimizerDimensionsFromAssignment(
     const Assignment* original_assignment, absl::Duration duration_limit,
     bool* time_limit_was_reached) {
@@ -4707,7 +4696,8 @@ void RoutingModel::CreateNeighborhoodOperators(
       CreateCPOperator<MakeChainInactiveOperator>();
   local_search_operators_[SWAP_ACTIVE] = CreateCPOperator<SwapActiveOperator>();
   local_search_operators_[SWAP_ACTIVE_CHAIN] =
-      CreateCPOperator<SwapActiveChainOperator>();
+      CreateCPOperatorWithArg<SwapActiveChainOperator, int>(
+          parameters.max_swap_active_chain_size());
   local_search_operators_[EXTENDED_SWAP_ACTIVE] =
       CreateCPOperator<ExtendedSwapActiveOperator>();
   std::vector<std::vector<int64_t>> alternative_sets(disjunctions_.size());
diff --git a/ortools/routing/routing.h b/ortools/routing/routing.h
index 3ce1ff1519..d80d093b3a 100644
--- a/ortools/routing/routing.h
+++ b/ortools/routing/routing.h
@@ -172,7 +172,6 @@
 
 #include "absl/container/flat_hash_map.h"
 #include "absl/container/flat_hash_set.h"
-#include "absl/container/inlined_vector.h"
 #include "absl/hash/hash.h"
 #include "absl/log/check.h"
 #include "absl/strings/string_view.h"
@@ -1521,36 +1520,21 @@ class RoutingModel {
   const Assignment* PackCumulsOfOptimizerDimensionsFromAssignment(
       const Assignment* original_assignment, absl::Duration duration_limit,
       bool* time_limit_was_reached = nullptr);
+
+#ifndef SWIG
   /// Contains the information needed by the solver to optimize a dimension's
   /// cumuls with travel-start dependent transit values.
   struct RouteDimensionTravelInfo {
     /// Contains the information for a single transition on the route.
     struct TransitionInfo {
-      /// The following struct defines a piecewise linear formulation, with
-      /// int64_t values for the "anchor" x and y values, and potential double
-      /// values for the slope of each linear function.
-      // TODO(user): Adjust the inlined vector sizes based on experiments.
-      struct PiecewiseLinearFormulation {
-        /// The set of *increasing* anchor cumul values for the interpolation.
-        absl::InlinedVector<int64_t, 8> x_anchors;
-        /// The y values used for the interpolation:
-        /// For any x anchor value, let i be an index such that
-        /// x_anchors[i] ≤ x < x_anchors[i+1], then the y value for x is
-        /// y_anchors[i] * (1-λ) + y_anchors[i+1] * λ, with
-        /// λ = (x - x_anchors[i]) / (x_anchors[i+1] - x_anchors[i]).
-        absl::InlinedVector<int64_t, 8> y_anchors;
-
-        std::string DebugString(std::string line_prefix = "") const;
-      };
-
       /// Models the (real) travel value Tᵣ, for this transition based on the
       /// departure value of the travel.
-      PiecewiseLinearFormulation travel_start_dependent_travel;
+      FloatSlopePiecewiseLinearFunction travel_start_dependent_travel;
 
       /// travel_compression_cost models the cost of the difference between the
       /// (real) travel value Tᵣ given by travel_start_dependent_travel and the
       /// compressed travel value considered in the scheduling.
-      PiecewiseLinearFormulation travel_compression_cost;
+      FloatSlopePiecewiseLinearFunction travel_compression_cost;
 
       /// The parts of the transit which occur pre/post travel between the
       /// nodes. The total transit between the two nodes i and j is
@@ -1580,6 +1564,8 @@ class RoutingModel {
     std::string DebugString(std::string line_prefix = "") const;
   };
 
+#endif  // SWIG
+
 #ifndef SWIG
   // TODO(user): Revisit if coordinates are added to the RoutingModel class.
   void SetSweepArranger(SweepArranger* sweep_arranger);
@@ -2345,6 +2331,11 @@ class RoutingModel {
   LocalSearchOperator* CreateCPOperator() {
     return CreateCPOperator(MakeLocalSearchOperator<T>);
   }
+  template <class T, typename ArgType>
+  LocalSearchOperator* CreateCPOperatorWithArg(ArgType arg) {
+    return CreateCPOperatorWithArg(MakeLocalSearchOperatorWithArg<T, ArgType>,
+                                   std::move(arg));
+  }
   using NeighborAccessor = std::function<const std::vector<int>&(int, int)>;
   template <class T>
   LocalSearchOperator* CreateCPOperatorWithNeighbors(
@@ -2368,6 +2359,15 @@ class RoutingModel {
                                             : vehicle_vars_,
         vehicle_start_class_callback_, std::move(get_neighbors));
   }
+  template <class T, typename ArgType>
+  LocalSearchOperator* CreateCPOperatorWithArg(const T& operator_factory,
+                                               ArgType arg) {
+    return operator_factory(solver_.get(), nexts_,
+                            CostsAreHomogeneousAcrossVehicles()
+                                ? std::vector<IntVar*>()
+                                : vehicle_vars_,
+                            vehicle_start_class_callback_, std::move(arg));
+  }
   template <class T, class Arg>
   LocalSearchOperator* CreateOperator(const Arg& arg) {
     return solver_->RevAlloc(new T(nexts_,
diff --git a/ortools/util/piecewise_linear_function.cc b/ortools/util/piecewise_linear_function.cc
index e691f9ae4a..23bc4eda23 100644
--- a/ortools/util/piecewise_linear_function.cc
+++ b/ortools/util/piecewise_linear_function.cc
@@ -14,15 +14,22 @@
 #include "ortools/util/piecewise_linear_function.h"
 
 #include <algorithm>
+#include <cstdint>
 #include <functional>
-#include <set>
 #include <string>
 #include <utility>
 #include <vector>
 
+#include "absl/algorithm/container.h"
 #include "absl/container/btree_set.h"
+#include "absl/container/inlined_vector.h"
+#include "absl/log/check.h"
 #include "absl/strings/str_format.h"
+#include "absl/strings/string_view.h"
+#include "ortools/base/dump_vars.h"
 #include "ortools/base/logging.h"
+#include "ortools/base/mathutil.h"
+#include "ortools/base/types.h"
 #include "ortools/util/saturated_arithmetic.h"
 
 namespace operations_research {
@@ -78,6 +85,7 @@ uint64_t UnsignedCapProd(uint64_t left, uint64_t right) {
 }
 }  // namespace
 
+// PiecewiseSegment
 PiecewiseSegment::PiecewiseSegment(int64_t point_x, int64_t point_y,
                                    int64_t slope, int64_t other_point_x)
     : slope_(slope), reference_x_(point_x), reference_y_(point_y) {
@@ -269,6 +277,7 @@ std::string PiecewiseSegment::DebugString() const {
   return result;
 }
 
+// PiecewiseLinearFunction
 const int PiecewiseLinearFunction::kNotFound = -1;
 
 PiecewiseLinearFunction::PiecewiseLinearFunction(
@@ -800,4 +809,60 @@ bool PiecewiseLinearFunction::IsNonIncreasingInternal() const {
   return true;
 }
 
+// FloatSlopePiecewiseLinearFunction
+const int FloatSlopePiecewiseLinearFunction::kNoValue = -1;
+
+FloatSlopePiecewiseLinearFunction::FloatSlopePiecewiseLinearFunction(
+    absl::InlinedVector<int64_t, 8> x_anchors,
+    absl::InlinedVector<int64_t, 8> y_anchors)
+    : x_anchors_(std::move(x_anchors)), y_anchors_(std::move(y_anchors)) {
+  DCHECK(absl::c_is_sorted(x_anchors_));
+  DCHECK_EQ(x_anchors_.size(), y_anchors_.size());
+  DCHECK_NE(x_anchors_.size(), 1);
+}
+
+std::string FloatSlopePiecewiseLinearFunction::DebugString(
+    absl::string_view line_prefix) const {
+  if (x_anchors_.size() <= 10) {
+    return "{ " + DUMP_VARS(x_anchors_, y_anchors_).str() + "}";
+  }
+  return absl::StrFormat("{\n%s%s\n%s%s\n}", line_prefix,
+                         DUMP_VARS(x_anchors_).str(), line_prefix,
+                         DUMP_VARS(y_anchors_).str());
+}
+
+int64_t FloatSlopePiecewiseLinearFunction::ComputeInBoundsValue(
+    int64_t x) const {
+  const int segment_index = GetSegmentIndex(x);
+  if (segment_index == kNoValue) return kNoValue;
+  return GetValueOnSegment(x, segment_index);
+}
+
+int64_t FloatSlopePiecewiseLinearFunction::ComputeConvexValue(int64_t x) const {
+  if (x_anchors_.empty()) return kNoValue;
+
+  int segment_index = kNoValue;
+  if (x <= x_anchors_[0]) {
+    segment_index = 0;
+  } else if (x >= x_anchors_.back()) {
+    segment_index = x_anchors_.size() - 2;
+  } else {
+    segment_index = GetSegmentIndex(x);
+  }
+
+  return GetValueOnSegment(x, segment_index);
+}
+
+int64_t FloatSlopePiecewiseLinearFunction::GetValueOnSegment(
+    int64_t x, int segment_index) const {
+  DCHECK_GE(segment_index, 0);
+  DCHECK_LE(segment_index, x_anchors_.size() - 2);
+  const double slope =
+      static_cast<double>(y_anchors_[segment_index + 1] -
+                          y_anchors_[segment_index]) /
+      (x_anchors_[segment_index + 1] - x_anchors_[segment_index]);
+  return MathUtil::Round<int64_t>(slope * (x - x_anchors_[segment_index]) +
+                                  y_anchors_[segment_index]);
+}
+
 }  // namespace operations_research
diff --git a/ortools/util/piecewise_linear_function.h b/ortools/util/piecewise_linear_function.h
index 26af5778bf..b076baddd8 100644
--- a/ortools/util/piecewise_linear_function.h
+++ b/ortools/util/piecewise_linear_function.h
@@ -22,10 +22,16 @@
 
 #include <cstdint>
 #include <functional>
+#include <iterator>
 #include <string>
 #include <utility>
 #include <vector>
 
+#include "absl/algorithm/container.h"
+#include "absl/container/inlined_vector.h"
+#include "absl/log/check.h"
+#include "absl/strings/string_view.h"
+
 namespace operations_research {
 // This structure stores one straight line. It contains the start point, the
 // end point and the slope.
@@ -268,5 +274,86 @@ class PiecewiseLinearFunction {
   bool is_non_decreasing_;
   bool is_non_increasing_;
 };
+
+// The following class defines a piecewise linear formulation with potential
+// double values for the slope of each linear function.
+// This formulation is meant to be used with a small number of segments (see
+// InlinedVector sizes below).
+// These segments are determined by int64_t values for the "anchor" x and y
+// values, such that (x_anchors_[i], y_anchors_[i]) and
+// (x_anchors_[i+1], y_anchors_[i+1]) are respectively the start and end point
+// of the i-th segment.
+// TODO(user): Adjust the inlined vector sizes based on experiments.
+class FloatSlopePiecewiseLinearFunction {
+ public:
+  static const int kNoValue;
+
+  FloatSlopePiecewiseLinearFunction() = default;
+  FloatSlopePiecewiseLinearFunction(absl::InlinedVector<int64_t, 8> x_anchors,
+                                    absl::InlinedVector<int64_t, 8> y_anchors);
+  FloatSlopePiecewiseLinearFunction(
+      FloatSlopePiecewiseLinearFunction&& other) noexcept {
+    *this = std::move(other);
+  }
+
+  FloatSlopePiecewiseLinearFunction& operator=(
+      FloatSlopePiecewiseLinearFunction&& other) noexcept {
+    x_anchors_ = std::move(other.x_anchors_);
+    y_anchors_ = std::move(other.y_anchors_);
+    return *this;
+  }
+
+  std::string DebugString(absl::string_view line_prefix = {}) const;
+
+  const absl::InlinedVector<int64_t, 8>& x_anchors() const {
+    return x_anchors_;
+  }
+  const absl::InlinedVector<int64_t, 8>& y_anchors() const {
+    return y_anchors_;
+  }
+
+  // Computes the y value associated to 'x'. Returns kNoValue if 'x' is out of
+  // bounds, i.e. lower than the first x_anchor and largest than the last.
+  int64_t ComputeInBoundsValue(int64_t x) const;
+
+  // Computes the y value associated to 'x'. Unlike ComputeInBoundsValue(), if
+  // 'x' is outside the bounds of the function, the function will still be
+  // defined by its outer segments.
+  int64_t ComputeConvexValue(int64_t x) const;
+
+ private:
+  // Returns the index of the segment x belongs to, i.e. the index i such that
+  // x_anchors_[i] ≤ x < x_anchors_[i+1]. For x = x_anchors_.back(), also
+  // returns the last segment (i.e. x_anchors_.size() - 2).
+  // Returns kNoValue if x is out of bounds for the function.
+  int GetSegmentIndex(int64_t x) const {
+    if (x_anchors_.empty() || x < x_anchors_[0] || x > x_anchors_.back()) {
+      return kNoValue;
+    }
+    if (x == x_anchors_.back()) return x_anchors_.size() - 2;
+
+    // Search for first element xi such that xi > x.
+    const auto upper_segment = absl::c_upper_bound(x_anchors_, x);
+    const int segment_index =
+        std::distance(x_anchors_.begin(), upper_segment) - 1;
+    DCHECK_GE(segment_index, 0);
+    DCHECK_LE(segment_index, x_anchors_.size() - 2);
+    return segment_index;
+  }
+
+  // Returns the value of 'x' on the linear segment determined by
+  // x_anchors_[segment_index] and x_anchors_[segment_index + 1].
+  int64_t GetValueOnSegment(int64_t x, int segment_index) const;
+
+  // The set of *increasing* anchor cumul values for the interpolation.
+  absl::InlinedVector<int64_t, 8> x_anchors_;
+  // The y values used for the interpolation:
+  // For any x anchor value, let i be an index such that
+  // x_anchors[i] ≤ x < x_anchors[i+1], then the y value for x is
+  // y_anchors[i] * (1-λ) + y_anchors[i+1] * λ, with
+  // λ = (x - x_anchors[i]) / (x_anchors[i+1] - x_anchors[i]).
+  absl::InlinedVector<int64_t, 8> y_anchors_;
+};
+
 }  // namespace operations_research
 #endif  // OR_TOOLS_UTIL_PIECEWISE_LINEAR_FUNCTION_H_