diff --git a/ortools/sat/cp_model.cc b/ortools/sat/cp_model.cc
index 8b206962fa..bcf4d1becf 100644
--- a/ortools/sat/cp_model.cc
+++ b/ortools/sat/cp_model.cc
@@ -612,6 +612,28 @@ Constraint CpModelBuilder::AddMinEquality(IntVar target,
   return Constraint(proto);
 }
 
+void CpModelBuilder::LinearExprToProto(const LinearExpr& expr,
+                                       LinearExpressionProto* expr_proto) {
+  for (const IntVar var : expr.variables()) {
+    expr_proto->add_vars(GetOrCreateIntegerIndex(var.index_));
+  }
+  for (const int64 coeff : expr.coefficients()) {
+    expr_proto->add_coeffs(coeff);
+  }
+  expr_proto->set_offset(expr.constant());
+}
+
+Constraint CpModelBuilder::AddLinMinEquality(
+    const LinearExpr& target, absl::Span<const LinearExpr> exprs) {
+  ConstraintProto* const proto = cp_model_.add_constraints();
+  LinearExprToProto(target, proto->mutable_lin_min()->mutable_target());
+  for (const LinearExpr& expr : exprs) {
+    LinearExpressionProto* expr_proto = proto->mutable_lin_min()->add_exprs();
+    LinearExprToProto(expr, expr_proto);
+  }
+  return Constraint(proto);
+}
+
 Constraint CpModelBuilder::AddMaxEquality(IntVar target,
                                           absl::Span<const IntVar> vars) {
   ConstraintProto* const proto = cp_model_.add_constraints();
@@ -622,6 +644,17 @@ Constraint CpModelBuilder::AddMaxEquality(IntVar target,
   return Constraint(proto);
 }
 
+Constraint CpModelBuilder::AddLinMaxEquality(
+    const LinearExpr& target, absl::Span<const LinearExpr> exprs) {
+  ConstraintProto* const proto = cp_model_.add_constraints();
+  LinearExprToProto(target, proto->mutable_lin_max()->mutable_target());
+  for (const LinearExpr& expr : exprs) {
+    LinearExpressionProto* expr_proto = proto->mutable_lin_max()->add_exprs();
+    LinearExprToProto(expr, expr_proto);
+  }
+  return Constraint(proto);
+}
+
 Constraint CpModelBuilder::AddDivisionEquality(IntVar target, IntVar numerator,
                                                IntVar denominator) {
   ConstraintProto* const proto = cp_model_.add_constraints();
diff --git a/ortools/sat/cp_model.h b/ortools/sat/cp_model.h
index 73e6e08f1b..2dafa94d46 100644
--- a/ortools/sat/cp_model.h
+++ b/ortools/sat/cp_model.h
@@ -737,9 +737,17 @@ class CpModelBuilder {
   /// Adds target == min(vars).
   Constraint AddMinEquality(IntVar target, absl::Span<const IntVar> vars);
 
+  /// Adds target == min(exprs).
+  Constraint AddLinMinEquality(const LinearExpr& target,
+                               absl::Span<const LinearExpr> exprs);
+
   /// Adds target == max(vars).
   Constraint AddMaxEquality(IntVar target, absl::Span<const IntVar> vars);
 
+  /// Adds target == max(exprs).
+  Constraint AddLinMaxEquality(const LinearExpr& target,
+                               absl::Span<const LinearExpr> exprs);
+
   /// Adds target = num / denom (integer division rounded towards 0).
   Constraint AddDivisionEquality(IntVar target, IntVar numerator,
                                  IntVar denominator);
@@ -812,6 +820,10 @@ class CpModelBuilder {
   friend class CumulativeConstraint;
   friend class ReservoirConstraint;
 
+  // Fills the 'expr_proto' with the linear expression represented by 'expr'.
+  void LinearExprToProto(const LinearExpr& expr,
+                         LinearExpressionProto* expr_proto);
+
   // Returns a (cached) integer variable index with a constant value.
   int IndexFromConstant(int64 value);
 
diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto
index 91d388065d..d525fca45c 100644
--- a/ortools/sat/cp_model.proto
+++ b/ortools/sat/cp_model.proto
@@ -76,6 +76,17 @@ message IntegerArgumentProto {
   repeated int32 vars = 2;
 }
 
+message LinearExpressionProto {
+  repeated int32 vars = 1;
+  repeated int32 coeffs = 2;
+  int64 offset = 3;
+}
+
+message LinearArgumentProto {
+  LinearExpressionProto target = 1;
+  repeated LinearExpressionProto exprs = 2;
+}
+
 // All variables must take different values.
 message AllDifferentConstraintProto {
   repeated int32 vars = 1;
@@ -246,7 +257,7 @@ message AutomatonConstraintProto {
   repeated int32 vars = 7;
 }
 
-// Next id: 27
+// Next id: 29
 message ConstraintProto {
   // For debug/logging only. Can be empty.
   string name = 1;
@@ -297,12 +308,22 @@ message ConstraintProto {
 
     // The int_max constraint forces the target to equal the maximum of all
     // variables.
+    // TODO(user): Remove int_max in favor of lin_max.
     IntegerArgumentProto int_max = 9;
 
+    // The lin_max constraint forces the target to equal the maximum of all
+    // linear expressions.
+    LinearArgumentProto lin_max = 27;
+
     // The int_min constraint forces the target to equal the minimum of all
     // variables.
+    // TODO(user): Remove int_min in favor of lin_min.
     IntegerArgumentProto int_min = 10;
 
+    // The lin_max constraint forces the target to equal the minimum of all
+    // linear expressions.
+    LinearArgumentProto lin_min = 28;
+
     // The int_min constraint forces the target to equal the product of all
     // variables.
     IntegerArgumentProto int_prod = 11;
diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc
index daeba5be2c..7ee056fbf7 100644
--- a/ortools/sat/cp_model_checker.cc
+++ b/ortools/sat/cp_model_checker.cc
@@ -194,6 +194,19 @@ std::string ValidateLinearConstraint(const CpModelProto& model,
   return "";
 }
 
+std::string ValidateLinearExpression(const CpModelProto& model,
+                                     const LinearExpressionProto& expr) {
+  if (expr.coeffs_size() != expr.vars_size()) {
+    return absl::StrCat("coeffs_size() != vars_size() in linear expression: ",
+                        ProtobufShortDebugString(expr));
+  }
+  if (PossibleIntegerOverflow(model, expr)) {
+    return absl::StrCat("Possible overflow in linear expression: ",
+                        ProtobufShortDebugString(expr));
+  }
+  return "";
+}
+
 std::string ValidateCircuitConstraint(const CpModelProto& model,
                                       const ConstraintProto& ct) {
   const int size = ct.circuit().tails().size();
@@ -387,6 +400,29 @@ std::string ValidateCpModel(const CpModelProto& model) {
         }
         RETURN_IF_NOT_EMPTY(ValidateLinearConstraint(model, ct));
         break;
+      case ConstraintProto::ConstraintCase::kLinMax: {
+        const std::string target_error =
+            ValidateLinearExpression(model, ct.lin_min().target());
+        if (!target_error.empty()) return target_error;
+        for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
+          const std::string expr_error =
+              ValidateLinearExpression(model, ct.lin_max().exprs(i));
+          if (!expr_error.empty()) return expr_error;
+        }
+        break;
+      }
+      case ConstraintProto::ConstraintCase::kLinMin: {
+        const std::string target_error =
+            ValidateLinearExpression(model, ct.lin_min().target());
+        if (!target_error.empty()) return target_error;
+        for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
+          const std::string expr_error =
+              ValidateLinearExpression(model, ct.lin_min().exprs(i));
+          if (!expr_error.empty()) return expr_error;
+        }
+        break;
+      }
+
       case ConstraintProto::ConstraintCase::kInterval:
         support_enforcement = true;
         RETURN_IF_NOT_EMPTY(ValidateIntervalConstraint(model, ct));
@@ -524,6 +560,25 @@ class ConstraintChecker {
     return max == actual_max;
   }
 
+  int64 LinearExpressionValue(const LinearExpressionProto& expr) {
+    int64 sum = expr.offset();
+    const int num_variables = expr.vars_size();
+    for (int i = 0; i < num_variables; ++i) {
+      sum += Value(expr.vars(i)) * expr.coeffs(i);
+    }
+    return sum;
+  }
+
+  bool LinMaxConstraintIsFeasible(const ConstraintProto& ct) {
+    const int64 max = LinearExpressionValue(ct.lin_max().target());
+    int64 actual_max = kint64min;
+    for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
+      const int64 expr_value = LinearExpressionValue(ct.lin_max().exprs(i));
+      actual_max = std::max(actual_max, expr_value);
+    }
+    return max == actual_max;
+  }
+
   bool IntProdConstraintIsFeasible(const ConstraintProto& ct) {
     const int64 prod = Value(ct.int_prod().target());
     int64 actual_prod = 1;
@@ -552,6 +607,16 @@ class ConstraintChecker {
     return min == actual_min;
   }
 
+  bool LinMinConstraintIsFeasible(const ConstraintProto& ct) {
+    const int64 min = LinearExpressionValue(ct.lin_min().target());
+    int64 actual_min = kint64max;
+    for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
+      const int64 expr_value = LinearExpressionValue(ct.lin_min().exprs(i));
+      actual_min = std::min(actual_min, expr_value);
+    }
+    return min == actual_min;
+  }
+
   bool AllDiffConstraintIsFeasible(const ConstraintProto& ct) {
     absl::flat_hash_set<int64> values;
     for (const int v : ct.all_diff().vars()) {
@@ -949,9 +1014,15 @@ bool SolutionIsFeasible(const CpModelProto& model,
       case ConstraintProto::ConstraintCase::kIntMin:
         is_feasible = checker.IntMinConstraintIsFeasible(ct);
         break;
+      case ConstraintProto::ConstraintCase::kLinMin:
+        is_feasible = checker.LinMinConstraintIsFeasible(ct);
+        break;
       case ConstraintProto::ConstraintCase::kIntMax:
         is_feasible = checker.IntMaxConstraintIsFeasible(ct);
         break;
+      case ConstraintProto::ConstraintCase::kLinMax:
+        is_feasible = checker.LinMaxConstraintIsFeasible(ct);
+        break;
       case ConstraintProto::ConstraintCase::kAllDiff:
         is_feasible = checker.AllDiffConstraintIsFeasible(ct);
         break;
diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc
index eec9b91ec2..40baa07d0a 100644
--- a/ortools/sat/cp_model_expand.cc
+++ b/ortools/sat/cp_model_expand.cc
@@ -817,6 +817,33 @@ void ExpandAutomaton(ConstraintProto* ct, PresolveContext* context) {
   ct->Clear();
 }
 
+// Fills the target as negated ref.
+void NegatedLinearExpression(const LinearExpressionProto& ref,
+                             LinearExpressionProto* target) {
+  for (int i = 0; i < ref.vars_size(); ++i) {
+    target->add_vars(NegatedRef(ref.vars(i)));
+    target->add_coeffs(ref.coeffs(i));
+  }
+  target->set_offset(-ref.offset());
+}
+
+void ExpandLinMax(ConstraintProto* ct, PresolveContext* context) {
+  ConstraintProto* const lin_min = context->working_model->add_constraints();
+  for (int i = 0; i < ct->enforcement_literal_size(); ++i) {
+    lin_min->add_enforcement_literal(ct->enforcement_literal(i));
+  }
+
+  // Target
+  NegatedLinearExpression(ct->lin_max().target(),
+                          lin_min->mutable_lin_min()->mutable_target());
+
+  for (int i = 0; i < ct->lin_max().exprs_size(); ++i) {
+    LinearExpressionProto* const expr = lin_min->mutable_lin_min()->add_exprs();
+    NegatedLinearExpression(ct->lin_max().exprs(i), expr);
+  }
+  ct->Clear();
+}
+
 }  // namespace
 
 void ExpandCpModel(PresolveOptions options, PresolveContext* context) {
@@ -838,6 +865,9 @@ void ExpandCpModel(PresolveOptions options, PresolveContext* context) {
       case ConstraintProto::ConstraintCase::kIntProd:
         ExpandIntProd(ct, context);
         break;
+      case ConstraintProto::ConstraintCase::kLinMax:
+        ExpandLinMax(ct, context);
+        break;
       case ConstraintProto::ConstraintCase::kElement:
         if (options.parameters.expand_element_constraints()) {
           ExpandElement(ct, context);
diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc
index f99a6ab943..2ee7319d7e 100644
--- a/ortools/sat/cp_model_lns.cc
+++ b/ortools/sat/cp_model_lns.cc
@@ -711,7 +711,9 @@ WeightedRandomRelaxationNeighborhoodGenerator::
       case ConstraintProto::kIntDiv:
       case ConstraintProto::kIntMod:
       case ConstraintProto::kIntMax:
+      case ConstraintProto::kLinMax:
       case ConstraintProto::kIntMin:
+      case ConstraintProto::kLinMin:
       case ConstraintProto::kNoOverlap:
       case ConstraintProto::kNoOverlap2D:
         constraint_weights_.push_back(2.0);
diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc
index 3e2e74dd72..057a52cbc9 100644
--- a/ortools/sat/cp_model_loader.cc
+++ b/ortools/sat/cp_model_loader.cc
@@ -19,6 +19,7 @@
 #include <set>
 #include <string>
 #include <utility>
+#include <vector>
 
 #include "absl/container/flat_hash_map.h"
 #include "absl/container/flat_hash_set.h"
@@ -30,6 +31,7 @@
 #include "ortools/sat/all_different.h"
 #include "ortools/sat/circuit.h"
 #include "ortools/sat/cp_constraints.h"
+#include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/cp_model_utils.h"
 #include "ortools/sat/cumulative.h"
 #include "ortools/sat/diffn.h"
@@ -1249,6 +1251,29 @@ void LoadIntMinConstraint(const ConstraintProto& ct, Model* m) {
   m->Add(IsEqualToMinOf(min, vars));
 }
 
+LinearExpression GetExprFromProto(const LinearExpressionProto& expr_proto,
+                                  const CpModelMapping& mapping) {
+  LinearExpression expr;
+  expr.vars = mapping.Integers(expr_proto.vars());
+  for (int j = 0; j < expr_proto.coeffs_size(); ++j) {
+    expr.coeffs.push_back(IntegerValue(expr_proto.coeffs(j)));
+  }
+  expr.offset = IntegerValue(expr_proto.offset());
+  return CanonicalizeExpr(expr);
+}
+
+void LoadLinMinConstraint(const ConstraintProto& ct, Model* m) {
+  auto* mapping = m->GetOrCreate<CpModelMapping>();
+  const LinearExpression min =
+      GetExprFromProto(ct.lin_min().target(), *mapping);
+  std::vector<LinearExpression> exprs;
+  exprs.reserve(ct.lin_min().exprs_size());
+  for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
+    exprs.push_back(GetExprFromProto(ct.lin_min().exprs(i), *mapping));
+  }
+  m->Add(IsEqualToMinOf(min, exprs));
+}
+
 void LoadIntMaxConstraint(const ConstraintProto& ct, Model* m) {
   auto* mapping = m->GetOrCreate<CpModelMapping>();
   const IntegerVariable max = mapping->Integer(ct.int_max().target());
@@ -1738,6 +1763,9 @@ bool LoadConstraint(const ConstraintProto& ct, Model* m) {
     case ConstraintProto::ConstraintProto::kIntMin:
       LoadIntMinConstraint(ct, m);
       return true;
+    case ConstraintProto::ConstraintProto::kLinMin:
+      LoadLinMinConstraint(ct, m);
+      return true;
     case ConstraintProto::ConstraintProto::kIntMax:
       LoadIntMaxConstraint(ct, m);
       return true;
diff --git a/ortools/sat/cp_model_loader.h b/ortools/sat/cp_model_loader.h
index d8f078b3fc..f932817838 100644
--- a/ortools/sat/cp_model_loader.h
+++ b/ortools/sat/cp_model_loader.h
@@ -227,6 +227,7 @@ void LoadAllDiffConstraint(const ConstraintProto& ct, Model* m);
 void LoadIntProdConstraint(const ConstraintProto& ct, Model* m);
 void LoadIntDivConstraint(const ConstraintProto& ct, Model* m);
 void LoadIntMinConstraint(const ConstraintProto& ct, Model* m);
+void LoadLinMinConstraint(const ConstraintProto& ct, Model* m);
 void LoadIntMaxConstraint(const ConstraintProto& ct, Model* m);
 void LoadNoOverlapConstraint(const ConstraintProto& ct, Model* m);
 void LoadNoOverlap2dConstraint(const ConstraintProto& ct, Model* m);
diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h
index a9172ad859..279ae7e4eb 100644
--- a/ortools/sat/cp_model_presolve.h
+++ b/ortools/sat/cp_model_presolve.h
@@ -104,6 +104,7 @@ class CpModelPresolver {
   bool PresolveInterval(int c, ConstraintProto* ct);
   bool PresolveIntDiv(ConstraintProto* ct);
   bool PresolveIntProd(ConstraintProto* ct);
+  // TODO(user) : Add presolve rules for LinMin.
   bool PresolveIntMin(ConstraintProto* ct);
   bool PresolveIntMax(ConstraintProto* ct);
   bool PresolveBoolXor(ConstraintProto* ct);
diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc
index bfd0968782..a5f7f14ee8 100644
--- a/ortools/sat/cp_model_solver.cc
+++ b/ortools/sat/cp_model_solver.cc
@@ -1968,6 +1968,7 @@ class LnsSolver : public SubSolver {
       SatParameters local_params(parameters_);
       local_params.set_max_deterministic_time(data.deterministic_limit);
       local_params.set_stop_after_first_solution(false);
+      local_params.set_log_search_progress(false);
 
       if (FLAGS_cp_model_dump_lns) {
         const std::string name =
diff --git a/ortools/sat/cp_model_utils.cc b/ortools/sat/cp_model_utils.cc
index f10dee9f14..e0ab6e187e 100644
--- a/ortools/sat/cp_model_utils.cc
+++ b/ortools/sat/cp_model_utils.cc
@@ -60,10 +60,24 @@ IndexReferences GetReferencesUsedByConstraint(const ConstraintProto& ct) {
       output.variables.push_back(ct.int_max().target());
       AddIndices(ct.int_max().vars(), &output.variables);
       break;
+    case ConstraintProto::ConstraintCase::kLinMax: {
+      AddIndices(ct.lin_max().target().vars(), &output.variables);
+      for (int i = 0; i < ct.lin_max().exprs_size(); ++i) {
+        AddIndices(ct.lin_max().exprs(i).vars(), &output.variables);
+      }
+      break;
+    }
     case ConstraintProto::ConstraintCase::kIntMin:
       output.variables.push_back(ct.int_min().target());
       AddIndices(ct.int_min().vars(), &output.variables);
       break;
+    case ConstraintProto::ConstraintCase::kLinMin: {
+      AddIndices(ct.lin_min().target().vars(), &output.variables);
+      for (int i = 0; i < ct.lin_min().exprs_size(); ++i) {
+        AddIndices(ct.lin_min().exprs(i).vars(), &output.variables);
+      }
+      break;
+    }
     case ConstraintProto::ConstraintCase::kIntProd:
       output.variables.push_back(ct.int_prod().target());
       AddIndices(ct.int_prod().vars(), &output.variables);
@@ -155,8 +169,12 @@ void ApplyToAllLiteralIndices(const std::function<void(int*)>& f,
       break;
     case ConstraintProto::ConstraintCase::kIntMax:
       break;
+    case ConstraintProto::ConstraintCase::kLinMax:
+      break;
     case ConstraintProto::ConstraintCase::kIntMin:
       break;
+    case ConstraintProto::ConstraintCase::kLinMin:
+      break;
     case ConstraintProto::ConstraintCase::kIntProd:
       break;
     case ConstraintProto::ConstraintCase::kLinear:
@@ -218,10 +236,22 @@ void ApplyToAllVariableIndices(const std::function<void(int*)>& f,
       APPLY_TO_SINGULAR_FIELD(int_max, target);
       APPLY_TO_REPEATED_FIELD(int_max, vars);
       break;
+    case ConstraintProto::ConstraintCase::kLinMax:
+      APPLY_TO_REPEATED_FIELD(lin_max, target()->mutable_vars);
+      for (int i = 0; i < ct->lin_max().exprs_size(); ++i) {
+        APPLY_TO_REPEATED_FIELD(lin_max, exprs(i)->mutable_vars);
+      }
+      break;
     case ConstraintProto::ConstraintCase::kIntMin:
       APPLY_TO_SINGULAR_FIELD(int_min, target);
       APPLY_TO_REPEATED_FIELD(int_min, vars);
       break;
+    case ConstraintProto::ConstraintCase::kLinMin:
+      APPLY_TO_REPEATED_FIELD(lin_min, target()->mutable_vars);
+      for (int i = 0; i < ct->lin_min().exprs_size(); ++i) {
+        APPLY_TO_REPEATED_FIELD(lin_min, exprs(i)->mutable_vars);
+      }
+      break;
     case ConstraintProto::ConstraintCase::kIntProd:
       APPLY_TO_SINGULAR_FIELD(int_prod, target);
       APPLY_TO_REPEATED_FIELD(int_prod, vars);
@@ -292,8 +322,12 @@ void ApplyToAllIntervalIndices(const std::function<void(int*)>& f,
       break;
     case ConstraintProto::ConstraintCase::kIntMax:
       break;
+    case ConstraintProto::ConstraintCase::kLinMax:
+      break;
     case ConstraintProto::ConstraintCase::kIntMin:
       break;
+    case ConstraintProto::ConstraintCase::kLinMin:
+      break;
     case ConstraintProto::ConstraintCase::kIntProd:
       break;
     case ConstraintProto::ConstraintCase::kLinear:
@@ -353,8 +387,12 @@ std::string ConstraintCaseName(
       return "kIntMod";
     case ConstraintProto::ConstraintCase::kIntMax:
       return "kIntMax";
+    case ConstraintProto::ConstraintCase::kLinMax:
+      return "kLinMax";
     case ConstraintProto::ConstraintCase::kIntMin:
       return "kIntMin";
+    case ConstraintProto::ConstraintCase::kLinMin:
+      return "kLinMin";
     case ConstraintProto::ConstraintCase::kIntProd:
       return "kIntProd";
     case ConstraintProto::ConstraintCase::kLinear:
@@ -422,8 +460,12 @@ std::vector<int> UsedIntervals(const ConstraintProto& ct) {
       break;
     case ConstraintProto::ConstraintCase::kIntMax:
       break;
+    case ConstraintProto::ConstraintCase::kLinMax:
+      break;
     case ConstraintProto::ConstraintCase::kIntMin:
       break;
+    case ConstraintProto::ConstraintCase::kLinMin:
+      break;
     case ConstraintProto::ConstraintCase::kIntProd:
       break;
     case ConstraintProto::ConstraintCase::kLinear:
diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc
index 924db292bc..7400987c42 100644
--- a/ortools/sat/integer.cc
+++ b/ortools/sat/integer.cc
@@ -791,7 +791,7 @@ void IntegerTrail::AppendRelaxedLinearReason(
   for (const IntegerVariable var : vars) {
     tmp_indices_.push_back(vars_[var].current_trail_index);
   }
-  RelaxLinearReason(slack, coeffs, &tmp_indices_);
+  if (slack > 0) RelaxLinearReason(slack, coeffs, &tmp_indices_);
   for (const int i : tmp_indices_) {
     reason->push_back(IntegerLiteral::GreaterOrEqual(integer_trail_[i].var,
                                                      integer_trail_[i].bound));
@@ -810,7 +810,7 @@ void IntegerTrail::RelaxLinearReason(IntegerValue slack,
   // We start by filtering *trail_indices:
   // - remove all level zero entries.
   // - keep the one that cannot be relaxed.
-  // - move the other one the the relax_heap_ (and creating the heap).
+  // - move the other one to the relax_heap_ (and creating the heap).
   int new_size = 0;
   const int size = coeffs.size();
   const int num_vars = vars_.size();
diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc
index 7af4a21e10..04a3341f47 100644
--- a/ortools/sat/integer_expr.cc
+++ b/ortools/sat/integer_expr.cc
@@ -15,9 +15,12 @@
 
 #include <algorithm>
 #include <memory>
+#include <vector>
 
 #include "absl/container/flat_hash_map.h"
+#include "absl/memory/memory.h"
 #include "ortools/base/stl_util.h"
+#include "ortools/sat/integer.h"
 #include "ortools/util/sorted_interval_list.h"
 
 namespace operations_research {
@@ -291,6 +294,146 @@ void MinPropagator::RegisterWith(GenericLiteralWatcher* watcher) {
   watcher->WatchUpperBound(min_var_, id);
 }
 
+LinearExpression CanonicalizeExpr(const LinearExpression& expr) {
+  LinearExpression canonical_expr;
+  canonical_expr.offset = expr.offset;
+  for (int i = 0; i < expr.vars.size(); ++i) {
+    if (expr.coeffs[i] < 0) {
+      canonical_expr.vars.push_back(NegationOf(expr.vars[i]));
+      canonical_expr.coeffs.push_back(-expr.coeffs[i]);
+    } else {
+      canonical_expr.vars.push_back(expr.vars[i]);
+      canonical_expr.coeffs.push_back(expr.coeffs[i]);
+    }
+  }
+  return canonical_expr;
+}
+
+IntegerValue LinExprLowerBound(const LinearExpression& expr,
+                               const IntegerTrail& integer_trail) {
+  IntegerValue lower_bound = expr.offset;
+  for (int i = 0; i < expr.vars.size(); ++i) {
+    DCHECK_GE(expr.coeffs[i], 0) << "The expression is not canonicalized";
+    lower_bound += expr.coeffs[i] * integer_trail.LowerBound(expr.vars[i]);
+  }
+  return lower_bound;
+}
+
+IntegerValue LinExprUpperBound(const LinearExpression& expr,
+                               const IntegerTrail& integer_trail) {
+  IntegerValue upper_bound = expr.offset;
+  for (int i = 0; i < expr.vars.size(); ++i) {
+    DCHECK_GE(expr.coeffs[i], 0) << "The expression is not canonicalized";
+    upper_bound += expr.coeffs[i] * integer_trail.UpperBound(expr.vars[i]);
+  }
+  return upper_bound;
+}
+
+LinMinPropagator::LinMinPropagator(const std::vector<LinearExpression>& exprs,
+                                   IntegerVariable min_var, Model* model)
+    : exprs_(exprs),
+      min_var_(min_var),
+      model_(model),
+      integer_trail_(model_->GetOrCreate<IntegerTrail>()) {}
+
+bool LinMinPropagator::Propagate() {
+  ub_constraints_.resize(rev_ub_constraints_size_);
+
+  if (exprs_.empty()) return true;
+
+  expr_lbs_.clear();
+
+  // Count the number of interval that are possible candidate for the min.
+  // Only the intervals for which lb > current_min_ub cannot.
+  const IntegerLiteral min_ub_literal =
+      integer_trail_->UpperBoundAsLiteral(min_var_);
+  const IntegerValue current_min_ub = integer_trail_->UpperBound(min_var_);
+  int num_intervals_that_can_be_min = 0;
+  int last_possible_min_interval = 0;
+
+  IntegerValue min_of_linear_expression_lb = kMaxIntegerValue;
+  for (int i = 0; i < exprs_.size(); ++i) {
+    const IntegerValue lb = LinExprLowerBound(exprs_[i], *integer_trail_);
+    expr_lbs_.push_back(lb);
+    min_of_linear_expression_lb = std::min(min_of_linear_expression_lb, lb);
+    if (lb <= current_min_ub) {
+      ++num_intervals_that_can_be_min;
+      last_possible_min_interval = i;
+    }
+  }
+
+  // Propagation a) lb(min) >= lb(MIN(exprs)) = MIN(lb(expr));
+
+  // Conflict will be detected by the fact that the [lb, ub] of the min is
+  // empty. In case of conflict, we just need the reason for pushing UB + 1.
+  if (min_of_linear_expression_lb > current_min_ub) {
+    min_of_linear_expression_lb = current_min_ub + 1;
+  }
+  if (min_of_linear_expression_lb > integer_trail_->LowerBound(min_var_)) {
+    integer_reason_.clear();
+    for (int i = 0; i < exprs_.size(); ++i) {
+      const IntegerValue slack = expr_lbs_[i] - min_of_linear_expression_lb;
+      integer_trail_->AppendRelaxedLinearReason(
+          slack, exprs_[i].coeffs, exprs_[i].vars, &integer_reason_);
+    }
+    if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(
+                                     min_var_, min_of_linear_expression_lb),
+                                 {}, integer_reason_)) {
+      return false;
+    }
+  }
+
+  // Propagation b) ub(min) >= ub(MIN(vars)) and we can't propagate anything
+  // here unless there is just one possible expression 'e' that can be the min:
+  //   for all u != e, lb(u) > ub(min);
+  // In this case, ub(min) >= ub(e).
+  if (num_intervals_that_can_be_min == 1) {
+    const IntegerValue ub_of_only_candidate =
+        LinExprUpperBound(exprs_[last_possible_min_interval], *integer_trail_);
+    if (current_min_ub < ub_of_only_candidate) {
+      integer_reason_.clear();
+
+      // The reason is that all the other interval start after current_min_ub.
+      // And that min_ub has its current value.
+      integer_reason_.push_back(min_ub_literal);
+      for (int i = 0; i < exprs_.size(); ++i) {
+        if (i == last_possible_min_interval) continue;
+        const IntegerValue slack = expr_lbs_[i] - (current_min_ub + 1);
+        integer_trail_->AppendRelaxedLinearReason(
+            slack, exprs_[i].coeffs, exprs_[i].vars, &integer_reason_);
+      }
+      std::vector<Literal> literal_reason;
+      integer_trail_->MergeReasonInto(integer_reason_, &literal_reason);
+      auto constraint = absl::make_unique<IntegerSumLE>(
+          literal_reason, exprs_[last_possible_min_interval].vars,
+          exprs_[last_possible_min_interval].coeffs,
+          current_min_ub - exprs_[last_possible_min_interval].offset, model_);
+      constraint->Propagate();
+      ub_constraints_.emplace_back(std::move(constraint));
+      rev_ub_constraints_size_ = ub_constraints_.size();
+    }
+  }
+
+  return true;
+}
+
+void LinMinPropagator::RegisterWith(GenericLiteralWatcher* watcher) {
+  const int id = watcher->Register(this);
+  for (const LinearExpression& expr : exprs_) {
+    for (int i = 0; i < expr.vars.size(); ++i) {
+      const IntegerVariable& var = expr.vars[i];
+      const IntegerValue coeff = expr.coeffs[i];
+      if (coeff > 0) {
+        watcher->WatchLowerBound(var, id);
+      } else {
+        watcher->WatchUpperBound(var, id);
+      }
+    }
+  }
+  watcher->RegisterReversibleInt(id, &rev_ub_constraints_size_);
+  watcher->WatchUpperBound(min_var_, id);
+}
+
 PositiveProductPropagator::PositiveProductPropagator(
     IntegerVariable a, IntegerVariable b, IntegerVariable p,
     IntegerTrail* integer_trail)
diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h
index dad78ce3f9..4c46107313 100644
--- a/ortools/sat/integer_expr.h
+++ b/ortools/sat/integer_expr.h
@@ -26,6 +26,7 @@
 #include "ortools/sat/precedences.h"
 #include "ortools/sat/sat_base.h"
 #include "ortools/sat/sat_solver.h"
+#include "ortools/util/rev.h"
 
 namespace operations_research {
 namespace sat {
@@ -142,6 +143,57 @@ class MinPropagator : public PropagatorInterface {
   DISALLOW_COPY_AND_ASSIGN(MinPropagator);
 };
 
+// Helper struct to model linear expression for lin_min/lin_max constraints. The
+// canonical expression should only contain positive coefficients.
+struct LinearExpression {
+  std::vector<IntegerVariable> vars;
+  std::vector<IntegerValue> coeffs;
+  IntegerValue offset = IntegerValue(0);
+};
+
+// Returns the same expression in the canonical form (all positive
+// coefficients).
+LinearExpression CanonicalizeExpr(const LinearExpression& expr);
+
+// Returns lower bound of linear expression using variable bounds of the
+// variables in expression. Assumes Canonical expression (all positive
+// coefficients).
+IntegerValue LinExprLowerBound(const LinearExpression& expr,
+                               const IntegerTrail& integer_trail);
+
+// Returns upper bound of linear expression using variable bounds of the
+// variables in expression. Assumes Canonical expression (all positive
+// coefficients).
+IntegerValue LinExprUpperBound(const LinearExpression& expr,
+                               const IntegerTrail& integer_trail);
+
+// Same as MinPropagator except this works on min = MIN(exprs) where exprs are
+// linear expressions. It uses IntegerSumLE to propagate bounds on the exprs.
+// Assumes Canonical expressions (all positive coefficients).
+// TODO(user): Remove the dependency on IntegerSumLE.
+class LinMinPropagator : public PropagatorInterface {
+ public:
+  LinMinPropagator(const std::vector<LinearExpression>& exprs,
+                   IntegerVariable min_var, Model* model);
+  LinMinPropagator(const LinMinPropagator&) = delete;
+  LinMinPropagator& operator=(const LinMinPropagator&) = delete;
+
+  bool Propagate() final;
+  void RegisterWith(GenericLiteralWatcher* watcher);
+
+ private:
+  const std::vector<LinearExpression> exprs_;
+  const IntegerVariable min_var_;
+  std::vector<IntegerValue> expr_lbs_;
+  Model* model_;
+  IntegerTrail* integer_trail_;
+
+  int rev_ub_constraints_size_ = 0;
+  std::vector<std::unique_ptr<IntegerSumLE>> ub_constraints_;
+
+  std::vector<IntegerLiteral> integer_reason_;
+};
+
 // Propagates a * b = c. Basic version, we don't extract any special cases, and
 // we only propagates the bounds.
 //
@@ -557,6 +609,59 @@ inline std::function<void(Model*)> IsEqualToMinOf(
   };
 }
 
+// Expresses the fact that an existing integer variable is equal to the minimum
+// of linear expressions. Assumes Canonical expressions (all positive
+// coefficients).
+inline std::function<void(Model*)> IsEqualToMinOf(
+    const LinearExpression& min_expr,
+    const std::vector<LinearExpression>& exprs) {
+  return [=](Model* model) {
+    IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
+
+    IntegerVariable min_var;
+    if (min_expr.vars.size() == 1 &&
+        std::abs(min_expr.coeffs[0].value()) == 1) {
+      if (min_expr.coeffs[0].value() == 1) {
+        min_var = min_expr.vars[0];
+      } else {
+        min_var = NegationOf(min_expr.vars[0]);
+      }
+    } else {
+      // Create a new variable if the expression is not just a single variable.
+      IntegerValue min_lb = LinExprLowerBound(min_expr, *integer_trail);
+      IntegerValue min_ub = LinExprUpperBound(min_expr, *integer_trail);
+      min_var = integer_trail->AddIntegerVariable(min_lb, min_ub);
+
+      // min_var = min_expr
+      std::vector<IntegerVariable> min_sum_vars = min_expr.vars;
+      std::vector<int64> min_sum_coeffs;
+      for (IntegerValue coeff : min_expr.coeffs) {
+        min_sum_coeffs.push_back(coeff.value());
+      }
+      min_sum_vars.push_back(min_var);
+      min_sum_coeffs.push_back(-1);
+
+      model->Add(FixedWeightedSum(min_sum_vars, min_sum_coeffs,
+                                  -min_expr.offset.value()));
+    }
+    for (const LinearExpression& expr : exprs) {
+      // min_var <= expr
+      std::vector<IntegerVariable> vars = expr.vars;
+      std::vector<int64> coeffs;
+      for (IntegerValue coeff : expr.coeffs) {
+        coeffs.push_back(coeff.value());
+      }
+      vars.push_back(min_var);
+      coeffs.push_back(-1);
+      model->Add(WeightedSumGreaterOrEqual(vars, coeffs, -expr.offset.value()));
+    }
+
+    LinMinPropagator* constraint = new LinMinPropagator(exprs, min_var, model);
+    constraint->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
+    model->TakeOwnership(constraint);
+  };
+}
+
 // Expresses the fact that an existing integer variable is equal to the maximum
 // of other integer variables.
 inline std::function<void(Model*)> IsEqualToMaxOf(
@@ -653,7 +758,7 @@ inline std::function<void(Model*)> ProductConstraint(IntegerVariable a,
   };
 }
 
-// Adds the constraint: a / b = d.
+// Adds the constraint: a / b = c.
 inline std::function<void(Model*)> DivisionConstraint(IntegerVariable a,
                                                       IntegerVariable b,
                                                       IntegerVariable c) {
@@ -666,7 +771,7 @@ inline std::function<void(Model*)> DivisionConstraint(IntegerVariable a,
   };
 }
 
-// Adds the constraint: a / b = d where b is a constant.
+// Adds the constraint: a / b = c where b is a constant.
 inline std::function<void(Model*)> FixedDivisionConstraint(IntegerVariable a,
                                                            IntegerValue b,
                                                            IntegerVariable c) {
diff --git a/ortools/sat/linear_constraint.cc b/ortools/sat/linear_constraint.cc
index b847978060..11cc6442c6 100644
--- a/ortools/sat/linear_constraint.cc
+++ b/ortools/sat/linear_constraint.cc
@@ -213,5 +213,32 @@ void MakeAllVariablesPositive(LinearConstraint* constraint) {
   }
 }
 
+// TODO(user): it would be better if LinearConstraint natively supported
+// term and not two separated vectors. Fix?
+//
+// TODO(user): This is really similar to CleanTermsAndFillConstraint(), maybe
+// we should just make the later switch negative variable to positive ones to
+// avoid an extra linear scan on each new cuts.
+void CanonicalizeConstraint(LinearConstraint* ct) {
+  std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
+
+  const int size = ct->vars.size();
+  for (int i = 0; i < size; ++i) {
+    if (VariableIsPositive(ct->vars[i])) {
+      terms.push_back({ct->vars[i], ct->coeffs[i]});
+    } else {
+      terms.push_back({NegationOf(ct->vars[i]), -ct->coeffs[i]});
+    }
+  }
+  std::sort(terms.begin(), terms.end());
+
+  ct->vars.clear();
+  ct->coeffs.clear();
+  for (const auto& term : terms) {
+    ct->vars.push_back(term.first);
+    ct->coeffs.push_back(term.second);
+  }
+}
+
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/linear_constraint.h b/ortools/sat/linear_constraint.h
index b673f30b2c..eeee724aa5 100644
--- a/ortools/sat/linear_constraint.h
+++ b/ortools/sat/linear_constraint.h
@@ -148,6 +148,12 @@ void CleanTermsAndFillConstraint(
     std::vector<std::pair<IntegerVariable, IntegerValue>>* terms,
     LinearConstraint* constraint);
 
+// Sorts the terms and makes all IntegerVariable positive. This assumes that a
+// variable or its negation only appear once.
+//
+// Note that currently this allocates some temporary memory.
+void CanonicalizeConstraint(LinearConstraint* ct);
+
 }  // namespace sat
 }  // namespace operations_research
 
diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc
index 4f4149ae98..e635652459 100644
--- a/ortools/sat/linear_constraint_manager.cc
+++ b/ortools/sat/linear_constraint_manager.cc
@@ -38,29 +38,6 @@ size_t ComputeHashOfTerms(const LinearConstraint& ct) {
   return hash;
 }
 
-// TODO(user): it would be better if LinearConstraint natively supported
-// term and not two separated vectors. Fix?
-void CanonicalizeConstraint(LinearConstraint* ct) {
-  std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
-
-  const int size = ct->vars.size();
-  for (int i = 0; i < size; ++i) {
-    if (VariableIsPositive(ct->vars[i])) {
-      terms.push_back({ct->vars[i], ct->coeffs[i]});
-    } else {
-      terms.push_back({NegationOf(ct->vars[i]), -ct->coeffs[i]});
-    }
-  }
-  std::sort(terms.begin(), terms.end());
-
-  ct->vars.clear();
-  ct->coeffs.clear();
-  for (const auto& term : terms) {
-    ct->vars.push_back(term.first);
-    ct->coeffs.push_back(term.second);
-  }
-}
-
 }  // namespace
 
 LinearConstraintManager::~LinearConstraintManager() {
@@ -132,7 +109,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
 // to detect duplicate constraints and merge bounds. This is also relevant if
 // we regenerate identical cuts for some reason.
 LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
-    LinearConstraint ct) {
+    LinearConstraint ct, bool* added) {
   CHECK(!ct.vars.empty());
   SimplifyConstraint(&ct);
   DivideByGCD(&ct);
@@ -145,25 +122,28 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
     const ConstraintIndex ct_index = equiv_constraints_[key];
     if (constraint_infos_[ct_index].constraint.vars == ct.vars &&
         constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) {
+      if (added != nullptr) *added = false;
       if (ct.lb > constraint_infos_[ct_index].constraint.lb) {
         if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
         constraint_infos_[ct_index].constraint.lb = ct.lb;
+        if (added != nullptr) *added = true;
       }
       if (ct.ub < constraint_infos_[ct_index].constraint.ub) {
         if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
         constraint_infos_[ct_index].constraint.ub = ct.ub;
+        if (added != nullptr) *added = true;
       }
       ++num_merged_constraints_;
       return ct_index;
     }
   }
 
+  if (added != nullptr) *added = true;
   const ConstraintIndex ct_index(constraint_infos_.size());
   ConstraintInfo ct_info;
   ct_info.constraint = std::move(ct);
   ct_info.l2_norm = ComputeL2Norm(ct_info.constraint);
   ct_info.is_in_lp = false;
-  ct_info.is_cut = false;
   ct_info.objective_parallelism_computed = false;
   ct_info.objective_parallelism = 0.0;
   ct_info.inactive_count = 0;
@@ -227,19 +207,20 @@ void LinearConstraintManager::AddCut(
   // Only add cut with sufficient efficacy.
   if (violation / l2_norm < 1e-5) return;
 
-  VLOG(1) << "Cut '" << type_name << "'"
-          << " size=" << ct.vars.size()
-          << " max_magnitude=" << ComputeInfinityNorm(ct) << " norm=" << l2_norm
-          << " violation=" << violation << " eff=" << violation / l2_norm;
-
   // Add the constraint. We only mark the constraint as a cut if it is not an
   // update of an already existing one.
-  const int64 prev_size = constraint_infos_.size();
-  const ConstraintIndex ct_index = Add(std::move(ct));
-  if (prev_size + 1 == constraint_infos_.size()) {
+  bool added = false;
+  const ConstraintIndex ct_index = Add(std::move(ct), &added);
+  if (added) {
+    VLOG(1) << "Cut '" << type_name << "'"
+            << " size=" << constraint_infos_[ct_index].constraint.vars.size()
+            << " max_magnitude="
+            << ComputeInfinityNorm(constraint_infos_[ct_index].constraint)
+            << " norm=" << l2_norm << " violation=" << violation
+            << " eff=" << violation / l2_norm;
+
     num_cuts_++;
     type_to_num_cuts_[type_name]++;
-    constraint_infos_[ct_index].is_cut = true;
   }
 }
 
diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h
index f88f56811b..59ca0a548f 100644
--- a/ortools/sat/linear_constraint_manager.h
+++ b/ortools/sat/linear_constraint_manager.h
@@ -45,7 +45,6 @@ class LinearConstraintManager {
   struct ConstraintInfo {
     LinearConstraint constraint;
     double l2_norm;
-    bool is_cut;
     int64 inactive_count;
     double objective_parallelism;
     bool objective_parallelism_computed;
@@ -66,9 +65,11 @@ class LinearConstraintManager {
 
   // Add a new constraint to the manager. Note that we canonicalize constraints
   // and merge the bounds of constraints with the same terms. We also perform
-  // basic preprocessing.
+  // basic preprocessing. If added is given, it will be set to true if this
+  // constraint was actually a new one and to false if it was dominated by an
+  // already existing one.
   DEFINE_INT_TYPE(ConstraintIndex, int32);
-  ConstraintIndex Add(LinearConstraint ct);
+  ConstraintIndex Add(LinearConstraint ct, bool* added = nullptr);
 
   // Same as Add(), but logs some information about the newly added constraint.
   // Cuts are also handled slightly differently than normal constraints.
diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc
index 5e4a94fce4..e33479cb2c 100644
--- a/ortools/sat/linear_programming_constraint.cc
+++ b/ortools/sat/linear_programming_constraint.cc
@@ -1574,8 +1574,13 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
     // If we add the row to the dense_terms, diff will indicate by how much
     // |upper_bound - ImpliedLB(dense_terms)| will change. That correspond to
     // increasing the multiplier by 1.
-    IntegerValue positive_diff = row_bound;
-    IntegerValue negative_diff = row_bound;
+    //
+    // At this stage, we are not sure computing sum coeff * bound will not
+    // overflow, so we use floating point numbers. It is fine to do so since
+    // this is not directly involved in the actual exact constraint generation:
+    // these variables are just used in an heuristic.
+    double positive_diff = ToDouble(row_bound);
+    double negative_diff = ToDouble(row_bound);
 
     // TODO(user): we could relax a bit some of the condition and allow a sign
     // change. It is just trickier to compute the diff when we allow such
@@ -1600,11 +1605,11 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
         positive_limit = std::min(positive_limit, overflow_limit);
         negative_limit = std::min(negative_limit, overflow_limit);
         if (coeff > 0) {
-          positive_diff -= coeff * lb;
-          negative_diff -= coeff * ub;
+          positive_diff -= ToDouble(coeff) * ToDouble(lb);
+          negative_diff -= ToDouble(coeff) * ToDouble(ub);
         } else {
-          positive_diff -= coeff * ub;
-          negative_diff -= coeff * lb;
+          positive_diff -= ToDouble(coeff) * ToDouble(ub);
+          negative_diff -= ToDouble(coeff) * ToDouble(lb);
         }
         continue;
       }
@@ -1629,20 +1634,24 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
       // This is how diff change.
       const IntegerValue implied = current > 0 ? lb : ub;
       if (implied != 0) {
-        positive_diff -= coeff * implied;
-        negative_diff -= coeff * implied;
+        positive_diff -= ToDouble(coeff) * ToDouble(implied);
+        negative_diff -= ToDouble(coeff) * ToDouble(implied);
       }
     }
 
     // Only add a multiple of this row if it tighten the final constraint.
+    // The positive_diff/negative_diff are supposed to be integer modulo the
+    // double precision, so we only add a multiple if they seems far away from
+    // zero.
     IntegerValue to_add(0);
-    if (positive_diff < 0 && positive_limit > 0) {
+    if (positive_diff <= -1.0 && positive_limit > 0) {
       to_add = positive_limit;
     }
-    if (negative_diff > 0 && negative_limit > 0) {
+    if (negative_diff >= 1.0 && negative_limit > 0) {
       // Pick this if it is better than the positive sign.
-      if (to_add == 0 || IntTypeAbs(negative_limit * negative_diff) >
-                             IntTypeAbs(positive_limit * positive_diff)) {
+      if (to_add == 0 ||
+          std::abs(ToDouble(negative_limit) * negative_diff) >
+              std::abs(ToDouble(positive_limit) * positive_diff)) {
         to_add = -negative_limit;
       }
     }
diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc
index 77967ee0c4..5ee4a38884 100644
--- a/ortools/sat/optimization.cc
+++ b/ortools/sat/optimization.cc
@@ -61,7 +61,7 @@ class Logger {
   explicit Logger(LogBehavior v) : use_stdout_(v == STDOUT_LOG) {}
   void Log(const std::string& message) {
     if (use_stdout_) {
-      printf("%s\n", message.c_str());
+      absl::PrintF("%s\n", message);
     } else {
       LOG(INFO) << message;
     }
diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc
index 1e77fd2752..e0e6b15f8f 100644
--- a/ortools/sat/presolve_context.cc
+++ b/ortools/sat/presolve_context.cc
@@ -50,7 +50,8 @@ void PresolveContext::AddImplication(int a, int b) {
 }
 
 // b => x in [lb, ub].
-void PresolveContext::AddImplyInDomain(int b, int x, const Domain& domain) {
+ConstraintProto* PresolveContext::AddImplyInDomain(int b, int x,
+                                                   const Domain& domain) {
   ConstraintProto* const imply = working_model->add_constraints();
 
   // Doing it like this seems to use slightly less memory.
@@ -60,6 +61,7 @@ void PresolveContext::AddImplyInDomain(int b, int x, const Domain& domain) {
   mutable_linear->mutable_vars()->Resize(1, x);
   mutable_linear->mutable_coeffs()->Resize(1, 1);
   FillDomainInProto(domain, mutable_linear);
+  return imply;
 }
 
 bool PresolveContext::DomainIsEmpty(int ref) const {
@@ -424,9 +426,16 @@ void PresolveContext::InsertVarValueEncoding(int literal, int ref,
         }
       }
     } else {
-      AddImplyInDomain(literal, var, Domain(var_value));
-      AddImplyInDomain(NegatedRef(literal), var,
-                       Domain(var_value).Complement());
+      VLOG(2) << "Insert lit(" << literal << ") <=> var(" << var
+              << ") == " << value;
+      const std::pair<int, int64> key{var, var_value};
+      eq_half_encoding[key][literal] =
+          AddImplyInDomain(literal, var, Domain(var_value));
+      VLOG(2) << "  imply_eq = " << eq_half_encoding[key][literal];
+      neq_half_encoding[key][NegatedRef(literal)] = AddImplyInDomain(
+          NegatedRef(literal), var, Domain(var_value).Complement());
+      VLOG(2) << "  imply_neq = "
+              << neq_half_encoding[key][NegatedRef(literal)];
     }
   } else {
     const int previous_literal = insert.first->second;
@@ -450,40 +459,38 @@ bool PresolveContext::InsertHalfVarValueEncoding(int literal, int var,
   const auto& new_info = direct_map.insert(std::make_pair(literal, ct));
   if (new_info.second) {
     VLOG(2) << "Collect lit(" << literal << ") implies var(" << var
-            << (imply_eq ? ") == " : ") != ") << value;
+            << (imply_eq ? ") == " : ") != ") << value << " [" << ct << "]";
     UpdateRuleStats("variables: detect half reified value encoding");
 
     if (gtl::ContainsKey(other_map, NegatedRef(literal))) {
-      const int ref_lit = imply_eq ? literal : NegatedRef(literal);
-      const auto insert_encoding_status =
-          encoding.insert(std::make_pair(std::make_pair(var, value), ref_lit));
-      if (insert_encoding_status.second) {
-        VLOG(2) << "Detect and store lit(" << ref_lit << ") <=> var(" << var
-                << ") == " << value;
-        UpdateRuleStats("variables: detect fully reified value encoding");
-      } else if (ref_lit != insert_encoding_status.first->second) {
-        const int previous_lit = insert_encoding_status.first->second;
-        VLOG(2) << "Detect duplicate encoding lit(" << ref_lit << ") == lit("
-                << previous_lit << ") <=> var(" << var << ") == " << value;
-        StoreBooleanEqualityRelation(previous_lit, ref_lit);
+      const int imply_eq_literal = imply_eq ? literal : NegatedRef(literal);
 
-        const AffineRelation::Relation r = GetAffineRelation(previous_lit);
-        if (r.representative != previous_lit) {
-          const int new_ref_lit =
-              r.coeff > 0 ? r.representative : NegatedRef(r.representative);
-          VLOG(2)
-              << "Different representative of previous literal, updating to "
-              << new_ref_lit;
+      const auto insert_encoding_status =
+          encoding.insert(std::make_pair(key, imply_eq_literal));
+      if (insert_encoding_status.second) {
+        VLOG(2) << "Detect and store lit(" << imply_eq_literal << ") <=> var("
+                << var << ") == " << value;
+        UpdateRuleStats("variables: detect fully reified value encoding");
+        encoding[key] = imply_eq_literal;
+      } else if (imply_eq_literal != insert_encoding_status.first->second) {
+        const int previous_imply_eq_literal =
+            insert_encoding_status.first->second;
+        VLOG(2) << "Detect duplicate encoding lit(" << imply_eq_literal
+                << ") == lit(" << previous_imply_eq_literal << ") <=> var("
+                << var << ") == " << value;
+        StoreBooleanEqualityRelation(imply_eq_literal,
+                                     previous_imply_eq_literal);
+
+        // Update reference literal.
+        const AffineRelation::Relation r =
+            GetAffineRelation(previous_imply_eq_literal);
+        const int new_ref_lit =
+            r.coeff > 0 ? r.representative : NegatedRef(r.representative);
+        if (new_ref_lit != previous_imply_eq_literal) {
+          VLOG(2) << "Updating reference encoding literal to " << new_ref_lit;
           encoding[key] = new_ref_lit;
         }
 
-        // Remove the two half encoding constraints from the model.
-        VLOG(2) << "Delete " << ct->DebugString();
-        ct->Clear();
-        direct_map.erase(literal);
-        VLOG(2) << "Delete " << other_map[NegatedRef(literal)]->DebugString();
-        other_map[NegatedRef(literal)]->Clear();
-        other_map.erase(NegatedRef(literal));
         UpdateRuleStats(
             "variables: merge equivalent var value encoding literals");
       }
@@ -534,6 +541,8 @@ int PresolveContext::GetOrCreateVarValueEncoding(int ref, int64 value) {
   const int64 var_min = MinOf(var);
   const int64 var_max = MaxOf(var);
   if (domains[var].Size() == 2) {
+    VLOG(2) << "GetOrCreate var(" << var << ") {" << var_min << ", " << var_max
+            << "}";
     // Checks if the other value is already encoded.
     const int64 other_value = var_value == var_min ? var_max : var_min;
     const std::pair<int, int64> other_key{var, other_value};
@@ -724,7 +733,8 @@ void PresolveContext::SubstituteVariableInObjective(
 
   // We can only "easily" substitute if the objective coefficient is a multiple
   // of the one in the constraint.
-  const int64 coeff_in_objective = gtl::FindOrDie(objective_map, var_in_equality);
+  const int64 coeff_in_objective =
+      gtl::FindOrDie(objective_map, var_in_equality);
   CHECK_NE(coeff_in_equality, 0);
   CHECK_EQ(coeff_in_objective % coeff_in_equality, 0);
   const int64 multiplier = coeff_in_objective / coeff_in_equality;
@@ -765,6 +775,11 @@ void PresolveContext::SubstituteVariableInObjective(
   // Tricky: The objective domain is without the offset, so we need to shift it
   objective_offset += static_cast<double>(offset.Min());
   objective_domain = objective_domain.AdditionWith(Domain(-offset.Min()));
+
+  // Because we can assume that the constraint we used was constraining
+  // (otherwise it would have been removed), the objective domain should be now
+  // constraining.
+  objective_domain_is_constraining = true;
 }
 
 void PresolveContext::WriteObjectiveToProto() {
diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h
index 847adc2065..c5666436ef 100644
--- a/ortools/sat/presolve_context.h
+++ b/ortools/sat/presolve_context.h
@@ -48,7 +48,7 @@ struct PresolveContext {
   void AddImplication(int a, int b);
 
   // b => x in [lb, ub].
-  void AddImplyInDomain(int b, int x, const Domain& domain);
+  ConstraintProto* AddImplyInDomain(int b, int x, const Domain& domain);
 
   // Helpers to query the current domain of a variable.
   bool DomainIsEmpty(int ref) const;