diff --git a/ortools/math_opt/core/BUILD.bazel b/ortools/math_opt/core/BUILD.bazel
index aa95a0c9e7..06da18fc67 100644
--- a/ortools/math_opt/core/BUILD.bazel
+++ b/ortools/math_opt/core/BUILD.bazel
@@ -105,6 +105,7 @@ cc_library(
     srcs = ["solver.cc"],
     hdrs = ["solver.h"],
     deps = [
+        ":base_solver",
         ":concurrent_calls_guard",
         ":math_opt_proto_utils",
         ":model_summary",
@@ -126,7 +127,6 @@ cc_library(
         "//ortools/math_opt/validators:result_validator",
         "//ortools/math_opt/validators:solve_parameters_validator",
         "//ortools/port:proto_utils",
-        "//ortools/util:solve_interrupter",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/memory",
         "@com_google_absl//absl/status",
@@ -230,3 +230,22 @@ cc_library(
         "@com_google_protobuf//:protobuf",
     ],
 )
+
+cc_library(
+    name = "base_solver",
+    srcs = ["base_solver.cc"],
+    hdrs = ["base_solver.h"],
+    deps = [
+        "//ortools/math_opt:callback_cc_proto",
+        "//ortools/math_opt:infeasible_subsystem_cc_proto",
+        "//ortools/math_opt:model_cc_proto",
+        "//ortools/math_opt:model_parameters_cc_proto",
+        "//ortools/math_opt:model_update_cc_proto",
+        "//ortools/math_opt:parameters_cc_proto",
+        "//ortools/math_opt:result_cc_proto",
+        "//ortools/util:solve_interrupter",
+        "@com_google_absl//absl/status:statusor",
+        "@com_google_absl//absl/strings",
+        "@com_google_absl//absl/strings:string_view",
+    ],
+)
diff --git a/ortools/math_opt/core/base_solver.cc b/ortools/math_opt/core/base_solver.cc
new file mode 100644
index 0000000000..a92d07f03c
--- /dev/null
+++ b/ortools/math_opt/core/base_solver.cc
@@ -0,0 +1,51 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/core/base_solver.h"
+
+#include <ostream>
+
+#include "absl/strings/str_cat.h"
+#include "absl/strings/string_view.h"
+
+namespace operations_research::math_opt {
+namespace {
+
+// Returns a string describing if the input pointer-like type is == nullptr.
+template <typename Ptr>
+absl::string_view IsNullString(const Ptr& ptr) {
+  return ptr == nullptr ? "<null>" : "<not null>";
+}
+
+}  // namespace
+
+std::ostream& operator<<(std::ostream& out, const BaseSolver::SolveArgs& args) {
+  out << "{ parameters: <" << absl::StrCat(args.parameters)
+      << ">, model_parameters: <" << absl::StrCat(args.model_parameters)
+      << ">, message_callback: " << IsNullString(args.message_callback)
+      << ", callback_registration: <"
+      << absl::StrCat(args.callback_registration)
+      << ">, user_cb: " << IsNullString(args.user_cb)
+      << ", interrupter: " << args.interrupter << " }";
+  return out;
+}
+
+std::ostream& operator<<(
+    std::ostream& out, const BaseSolver::ComputeInfeasibleSubsystemArgs& args) {
+  out << "{ parameters: <" << absl::StrCat(args.parameters)
+      << ">, message_callback: " << IsNullString(args.message_callback)
+      << ", interrupter: " << args.interrupter << " }";
+  return out;
+}
+
+}  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/core/base_solver.h b/ortools/math_opt/core/base_solver.h
new file mode 100644
index 0000000000..8768ac78d7
--- /dev/null
+++ b/ortools/math_opt/core/base_solver.h
@@ -0,0 +1,129 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_CORE_BASE_SOLVER_H_
+#define OR_TOOLS_MATH_OPT_CORE_BASE_SOLVER_H_
+
+#include <functional>
+#include <ostream>
+#include <string>
+#include <vector>
+
+#include "absl/status/statusor.h"
+#include "ortools/math_opt/callback.pb.h"
+#include "ortools/math_opt/infeasible_subsystem.pb.h"
+#include "ortools/math_opt/model.pb.h"
+#include "ortools/math_opt/model_parameters.pb.h"
+#include "ortools/math_opt/model_update.pb.h"
+#include "ortools/math_opt/parameters.pb.h"
+#include "ortools/math_opt/result.pb.h"
+#include "ortools/util/solve_interrupter.h"
+
+namespace operations_research::math_opt {
+
+// The API of solvers (in-process, sub-process and streaming RPC ones).
+//
+// Thread-safety: methods Solve() and Update() must not be called concurrently;
+// they should immediately return with an error status if this happens.
+//
+// TODO: b/350984134 - Rename `Solver` into `InProcessSolver` and then rename
+// `BaseSolver` into `Solver`.
+class BaseSolver {
+ public:
+  // Callback function for messages callback sent by the solver.
+  //
+  // Each message represents a single output line from the solver, and each
+  // message does not contain any '\n' character in it.
+  //
+  // Thread-safety: a callback may be called concurrently from multiple
+  // threads. The users is expected to use proper synchronization primitives to
+  // deal with that.
+  using MessageCallback = std::function<void(const std::vector<std::string>&)>;
+
+  // Callback function type for MIP/LP callbacks.
+  using Callback = std::function<CallbackResultProto(const CallbackDataProto&)>;
+
+  // Arguments used when calling Solve() to solve the problem.
+  struct SolveArgs {
+    SolveParametersProto parameters;
+    ModelSolveParametersProto model_parameters;
+
+    // An optional callback for messages emitted by the solver.
+    //
+    // When set it enables the solver messages and ignores the `enable_output`
+    // in solve parameters; messages are redirected to the callback and not
+    // printed on stdout/stderr/logs anymore.
+    MessageCallback message_callback = nullptr;
+
+    CallbackRegistrationProto callback_registration;
+    Callback user_cb = nullptr;
+
+    // An optional interrupter that the solver can use to interrupt the solve
+    // early.
+    const SolveInterrupter* interrupter = nullptr;
+
+    friend std::ostream& operator<<(std::ostream& out, const SolveArgs& args);
+  };
+
+  // Arguments used when calling ComputeInfeasibleSubsystem().
+  struct ComputeInfeasibleSubsystemArgs {
+    SolveParametersProto parameters;
+
+    // An optional callback for messages emitted by the solver.
+    //
+    // When set it enables the solver messages and ignores the `enable_output`
+    // in solve parameters; messages are redirected to the callback and not
+    // printed on stdout/stderr/logs anymore.
+    MessageCallback message_callback = nullptr;
+
+    // An optional interrupter that the solver can use to interrupt the solve
+    // early.
+    const SolveInterrupter* interrupter = nullptr;
+
+    friend std::ostream& operator<<(std::ostream& out,
+                                    const ComputeInfeasibleSubsystemArgs& args);
+  };
+
+  BaseSolver() = default;
+  BaseSolver(const BaseSolver&) = delete;
+  BaseSolver& operator=(const BaseSolver&) = delete;
+  virtual ~BaseSolver() = default;
+
+  // Solves the current model (including all updates).
+  virtual absl::StatusOr<SolveResultProto> Solve(
+      const SolveArgs& arguments) = 0;
+
+  // Computes an infeasible subsystem of `model` (including all updates).
+  virtual absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
+  ComputeInfeasibleSubsystem(
+      const ComputeInfeasibleSubsystemArgs& arguments) = 0;
+
+  // Updates the model to solve and returns true, or returns false if this
+  // update is not supported by the underlying solver.
+  //
+  // The model_update is passed by value. Non in-process implementations will
+  // move it in-place in the messages used to communicate with the other
+  // process. Thus if possible, the caller should std::move() this proto to this
+  // function.
+  //
+  // A status error will be returned if the model_update is invalid or the
+  // underlying solver has an internal error.
+  //
+  // When this function returns false, the BaseSolver object is in a failed
+  // state.
+  virtual absl::StatusOr<bool> Update(ModelUpdateProto model_update) = 0;
+};
+
+}  // namespace operations_research::math_opt
+
+#endif  // OR_TOOLS_MATH_OPT_CORE_BASE_SOLVER_H_
diff --git a/ortools/math_opt/core/solver.cc b/ortools/math_opt/core/solver.cc
index d447dadb18..b337ff2a4e 100644
--- a/ortools/math_opt/core/solver.cc
+++ b/ortools/math_opt/core/solver.cc
@@ -155,7 +155,7 @@ absl::StatusOr<SolveResultProto> Solver::Solve(const SolveArgs& arguments) {
   return result;
 }
 
-absl::StatusOr<bool> Solver::Update(const ModelUpdateProto& model_update) {
+absl::StatusOr<bool> Solver::Update(const ModelUpdateProto model_update) {
   ASSIGN_OR_RETURN(const auto guard,
                    ConcurrentCallsGuard::TryAcquire(concurrent_calls_tracker_));
 
diff --git a/ortools/math_opt/core/solver.h b/ortools/math_opt/core/solver.h
index 8c9181417d..92c5206264 100644
--- a/ortools/math_opt/core/solver.h
+++ b/ortools/math_opt/core/solver.h
@@ -14,12 +14,12 @@
 #ifndef OR_TOOLS_MATH_OPT_CORE_SOLVER_H_
 #define OR_TOOLS_MATH_OPT_CORE_SOLVER_H_
 
-#include <functional>
 #include <memory>
 
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
 #include "ortools/math_opt/callback.pb.h"
+#include "ortools/math_opt/core/base_solver.h"
 #include "ortools/math_opt/core/concurrent_calls_guard.h"
 #include "ortools/math_opt/core/model_summary.h"
 #include "ortools/math_opt/core/solver_interface.h"
@@ -29,10 +29,8 @@
 #include "ortools/math_opt/model_update.pb.h"
 #include "ortools/math_opt/parameters.pb.h"
 #include "ortools/math_opt/result.pb.h"
-#include "ortools/util/solve_interrupter.h"
 
-namespace operations_research {
-namespace math_opt {
+namespace operations_research::math_opt {
 
 // A solver for a given model and solver implementation.
 //
@@ -65,95 +63,15 @@ namespace math_opt {
 //   CHECK_OK(second_solution.status());
 //   // Use the second_solution of the updated problem here.
 //
-class Solver {
+class Solver : public BaseSolver {
  public:
   using InitArgs = SolverInterface::InitArgs;
 
-  // Callback function for messages callback sent by the solver.
-  //
-  // Each message represents a single output line from the solver, and each
-  // message does not contain any '\n' character in it.
-  //
-  // Thread-safety: a callback may be called concurrently from multiple
-  // threads. The users is expected to use proper synchronization primitives to
-  // deal with that.
-  using MessageCallback = SolverInterface::MessageCallback;
-
-  // Callback function type for MIP/LP callbacks.
-  using Callback = std::function<CallbackResultProto(const CallbackDataProto&)>;
-
-  // Arguments used when calling Solve() to solve the problem.
-  struct SolveArgs {
-    SolveParametersProto parameters;
-    ModelSolveParametersProto model_parameters;
-
-    // An optional callback for messages emitted by the solver.
-    //
-    // When set it enables the solver messages and ignores the `enable_output`
-    // in solve parameters; messages are redirected to the callback and not
-    // printed on stdout/stderr/logs anymore.
-    MessageCallback message_callback = nullptr;
-
-    CallbackRegistrationProto callback_registration;
-    Callback user_cb = nullptr;
-
-    // An optional interrupter that the solver can use to interrupt the solve
-    // early.
-    const SolveInterrupter* interrupter = nullptr;
-  };
-
-  // Arguments used when calling ComputeInfeasibleSubsystem().
-  struct ComputeInfeasibleSubsystemArgs {
-    SolveParametersProto parameters;
-
-    // An optional callback for messages emitted by the solver.
-    //
-    // When set it enables the solver messages and ignores the `enable_output`
-    // in solve parameters; messages are redirected to the callback and not
-    // printed on stdout/stderr/logs anymore.
-    MessageCallback message_callback = nullptr;
-
-    // An optional interrupter that the solver can use to interrupt the solve
-    // early.
-    const SolveInterrupter* interrupter = nullptr;
-  };
-
   // A shortcut for calling Solver::New() and then Solver::Solve().
   static absl::StatusOr<SolveResultProto> NonIncrementalSolve(
       const ModelProto& model, SolverTypeProto solver_type,
       const InitArgs& init_args, const SolveArgs& solve_args);
 
-  // Builds a solver of the given type with the provided model and
-  // initialization parameters.
-  static absl::StatusOr<std::unique_ptr<Solver>> New(
-      SolverTypeProto solver_type, const ModelProto& model,
-      const InitArgs& arguments);
-
-  Solver(const Solver&) = delete;
-  Solver& operator=(const Solver&) = delete;
-
-  ~Solver();
-
-  // Solves the current model (included all updates).
-  absl::StatusOr<SolveResultProto> Solve(const SolveArgs& arguments);
-
-  // Updates the model to solve and returns true, or returns false if this
-  // update is not supported by the underlying solver.
-  //
-  // A status error will be returned if the model_update is invalid or the
-  // underlying solver has an internal error.
-  //
-  // When this function returns false, the Solver object is in a failed
-  // state. In that case the underlying SolverInterface implementation has been
-  // destroyed (this enables the caller to instantiate a new Solver without
-  // destroying the previous one first even if they use Gurobi with a single-use
-  // license).
-  absl::StatusOr<bool> Update(const ModelUpdateProto& model_update);
-
-  // Computes an infeasible subsystem of `model`.
-  absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
-  ComputeInfeasibleSubsystem(const ComputeInfeasibleSubsystemArgs& arguments);
-
   // A shortcut for calling Solver::New() and then
   // Solver()::ComputeInfeasibleSubsystem()
   static absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
@@ -162,6 +80,29 @@ class Solver {
       const InitArgs& init_args,
       const ComputeInfeasibleSubsystemArgs& compute_infeasible_subsystem_args);
 
+  // Builds a solver of the given type with the provided model and
+  // initialization parameters.
+  static absl::StatusOr<std::unique_ptr<Solver>> New(
+      SolverTypeProto solver_type, const ModelProto& model,
+      const InitArgs& arguments);
+
+  ~Solver() override;
+
+  absl::StatusOr<SolveResultProto> Solve(const SolveArgs& arguments) override;
+
+  // See BaseSolver::Update.
+  //
+  // When this function returns false, the Solver object is in a failed
+  // state. In that case the underlying SolverInterface implementation has been
+  // destroyed (this enables the caller to instantiate a new Solver without
+  // destroying the previous one first even if they use Gurobi with a single-use
+  // license).
+  absl::StatusOr<bool> Update(ModelUpdateProto model_update) override;
+
+  absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
+  ComputeInfeasibleSubsystem(
+      const ComputeInfeasibleSubsystemArgs& arguments) override;
+
  private:
   Solver(std::unique_ptr<SolverInterface> underlying_solver,
          ModelSummary model_summary);
@@ -191,7 +132,6 @@ absl::Status ValidateInitArgs(const Solver::InitArgs& init_args,
                               SolverTypeProto solver_type);
 
 }  // namespace internal
-}  // namespace math_opt
-}  // namespace operations_research
+}  // namespace operations_research::math_opt
 
 #endif  // OR_TOOLS_MATH_OPT_CORE_SOLVER_H_
diff --git a/ortools/math_opt/cpp/BUILD.bazel b/ortools/math_opt/cpp/BUILD.bazel
index 39f5a6f2f0..a6063883b1 100644
--- a/ortools/math_opt/cpp/BUILD.bazel
+++ b/ortools/math_opt/cpp/BUILD.bazel
@@ -24,6 +24,7 @@ cc_library(
     deps = [
         ":model",
         ":solve",
+        ":solver_resources",
     ],
 )
 
@@ -52,17 +53,15 @@ cc_library(
         "//ortools/base:status_macros",
         "//ortools/math_opt:solution_cc_proto",
         "//ortools/math_opt:sparse_containers_cc_proto",
+        "//ortools/math_opt/constraints/quadratic:quadratic_constraint",
         "//ortools/math_opt/core:sparse_vector_view",
         "//ortools/math_opt/storage:model_storage",
         "//ortools/math_opt/storage:model_storage_types",
-        "//ortools/math_opt/validators:ids_validator",
         "//ortools/math_opt/validators:sparse_vector_validator",
-        "//ortools/util:status_macros",
         "@com_google_absl//absl/algorithm:container",
         "@com_google_absl//absl/container:flat_hash_map",
         "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
-        "@com_google_absl//absl/strings:string_view",
         "@com_google_absl//absl/types:span",
         "@com_google_protobuf//:protobuf",
     ],
@@ -168,13 +167,12 @@ cc_library(
         "//ortools/math_opt:result_cc_proto",
         "//ortools/math_opt:solution_cc_proto",
         "//ortools/math_opt:sparse_containers_cc_proto",
-        "//ortools/math_opt/core:sparse_vector_view",
+        "//ortools/math_opt/constraints/quadratic:quadratic_constraint",
         "//ortools/math_opt/storage:model_storage",
-        "//ortools/math_opt/validators:ids_validator",
-        "//ortools/math_opt/validators:sparse_vector_validator",
         "//ortools/util:status_macros",
         "@com_google_absl//absl/container:flat_hash_map",
         "@com_google_absl//absl/log",
+        "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
         "@com_google_absl//absl/strings",
@@ -224,6 +222,7 @@ cc_library(
         "//ortools/math_opt/storage:model_storage",
         "@com_google_absl//absl/algorithm:container",
         "@com_google_absl//absl/container:flat_hash_set",
+        "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
     ],
 )
@@ -352,30 +351,23 @@ cc_library(
     srcs = ["solve.cc"],
     hdrs = ["solve.h"],
     deps = [
-        ":callback",
         ":compute_infeasible_subsystem_arguments",
         ":compute_infeasible_subsystem_result",
-        ":enums",
+        ":incremental_solver",
         ":model",
-        ":model_solve_parameters",
         ":parameters",
         ":solve_arguments",
+        ":solve_impl",
         ":solve_result",
         ":solver_init_arguments",
         ":streamable_solver_init_arguments",
         ":update_result",
-        ":update_tracker",
-        "//ortools/base:status_macros",
         "//ortools/math_opt:callback_cc_proto",
         "//ortools/math_opt:infeasible_subsystem_cc_proto",
         "//ortools/math_opt:parameters_cc_proto",
+        "//ortools/math_opt/core:base_solver",
         "//ortools/math_opt/core:solver",
-        "//ortools/math_opt/storage:model_storage",
-        "//ortools/util:status_macros",
-        "@com_google_absl//absl/memory",
-        "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
-        "@com_google_absl//absl/synchronization",
     ],
 )
 
@@ -430,6 +422,7 @@ cc_library(
         ":variable_and_expressions",
         "//ortools/base:gmock",
         "//ortools/base:logging",
+        "//ortools/math_opt/constraints/quadratic:quadratic_constraint",
         "@com_google_absl//absl/container:flat_hash_map",
         "@com_google_absl//absl/log",
         "@com_google_absl//absl/log:check",
@@ -526,3 +519,45 @@ cc_library(
         "@com_google_absl//absl/strings:string_view",
     ],
 )
+
+cc_library(
+    name = "solve_impl",
+    srcs = ["solve_impl.cc"],
+    hdrs = ["solve_impl.h"],
+    deps = [
+        ":compute_infeasible_subsystem_arguments",
+        ":compute_infeasible_subsystem_result",
+        ":incremental_solver",
+        ":model",
+        ":parameters",
+        ":solve_arguments",
+        ":solve_result",
+        ":update_result",
+        ":update_tracker",
+        "//ortools/base:status_macros",
+        "//ortools/math_opt/core:base_solver",
+        "//ortools/math_opt/storage:model_storage",
+        "//ortools/util:solve_interrupter",
+        "//ortools/util:status_macros",
+        "@com_google_absl//absl/functional:any_invocable",
+        "@com_google_absl//absl/log:check",
+        "@com_google_absl//absl/memory",
+        "@com_google_absl//absl/status",
+        "@com_google_absl//absl/status:statusor",
+        "@com_google_absl//absl/synchronization",
+    ],
+)
+
+cc_library(
+    name = "incremental_solver",
+    hdrs = ["incremental_solver.h"],
+    deps = [
+        ":compute_infeasible_subsystem_arguments",
+        ":compute_infeasible_subsystem_result",
+        ":parameters",
+        ":solve_arguments",
+        ":solve_result",
+        ":update_result",
+        "@com_google_absl//absl/status:statusor",
+    ],
+)
diff --git a/ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h b/ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h
index 9d1b1d9410..6ccc4ca1c9 100644
--- a/ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h
+++ b/ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h
@@ -14,8 +14,8 @@
 #ifndef OR_TOOLS_MATH_OPT_CPP_COMPUTE_INFEASIBLE_SUBSYSTEM_ARGUMENTS_H_
 #define OR_TOOLS_MATH_OPT_CPP_COMPUTE_INFEASIBLE_SUBSYSTEM_ARGUMENTS_H_
 
-#include "ortools/math_opt/cpp/message_callback.h"    // IWYU pragma: export
-#include "ortools/math_opt/cpp/parameters.h"          // IWYU pragma: export
+#include "ortools/math_opt/cpp/message_callback.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/parameters.h"        // IWYU pragma: export
 #include "ortools/util/solve_interrupter.h"         // IWYU pragma: export
 
 namespace operations_research::math_opt {
diff --git a/ortools/math_opt/cpp/incremental_solver.h b/ortools/math_opt/cpp/incremental_solver.h
new file mode 100644
index 0000000000..bde0eb8872
--- /dev/null
+++ b/ortools/math_opt/cpp/incremental_solver.h
@@ -0,0 +1,176 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// IWYU pragma: private, include "ortools/math_opt/cpp/math_opt.h"
+// IWYU pragma: friend "ortools/math_opt/cpp/.*"
+
+#ifndef OR_TOOLS_MATH_OPT_CPP_INCREMENTAL_SOLVER_H_
+#define OR_TOOLS_MATH_OPT_CPP_INCREMENTAL_SOLVER_H_
+
+#include "absl/status/statusor.h"
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_result.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/parameters.h"       // IWYU pragma: export
+#include "ortools/math_opt/cpp/solve_arguments.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/solve_result.h"     // IWYU pragma: export
+#include "ortools/math_opt/cpp/update_result.h"    // IWYU pragma: export
+
+namespace operations_research::math_opt {
+
+// Interface for incrementally solving of a model.
+//
+// This is a feature for advance users. Most users should only use the
+// non-incremental Solve(), SubprocessSolve(),... functions.
+//
+// Here incremental means that the we try to reuse the existing underlying
+// solver internals between each solve. There is no guarantee though that the
+// solver supports all possible model changes. Hence there is not guarantee that
+// performances will be improved when using this class; this is solver
+// dependent. Typically LPs have more to gain from incremental solve than
+// MIPs. In both cases, even if the solver supports the model changes,
+// incremental solve may actually be slower.
+//
+// Implementation of this interface are returned by factories that can be found
+// next to non-incremental solve functions. See NewIncrementalSolver() in
+// solve.h and NewSubprocessIncrementalSolver() in subprocess_solve.h for
+// example.
+//
+// Those factories function instantiates the solver, setup it from the current
+// state of the Model and register on it to listen to changes. Calling Solve()
+// will update the underlying solver with latest model changes and solve this
+// model.
+//
+// Usage:
+//   Model model = ...;
+//   ASSIGN_OR_RETURN(
+//     const std::unique_ptr<IncrementalSolver> incremental_solve,
+//     NewIncrementalSolver(&model, SolverType::kXxx));
+//
+//   ASSIGN_OR_RETURN(const SolveResult result1, incremental_solve->Solve());
+//
+//   model.AddVariable(...);
+//   ...
+//
+//   ASSIGN_OR_RETURN(const SolveResult result2, incremental_solve->Solve());
+//
+//   ...
+//
+// Lifecycle: The IncrementalSolver is only keeping a std::weak_ref on Model's
+// internal data and thus it returns an error if Update() or Solve() are called
+// after the Model has been destroyed. It is fine to destroy the
+// IncrementalSolver after the associated Model though.
+//
+// Thread-safety: The New(), Solve() and Update() methods must not be called
+// while modifying the Model() (adding variables...). The user is expected to
+// use proper synchronization primitives to serialize changes to the model and
+// the use of this object. Note though that it is safe to call methods from
+// different IncrementalSolver instances on the same Model concurrently. Same
+// for calling factories of IncrementalSolver. The destructor is thread-safe and
+// can be called even during a modification of the Model.
+//
+// There is no problem calling SolveWithoutUpdate() concurrently on different
+// instances of IncrementalSolver or while the model is being modified (unless
+// of course the underlying solver itself is not thread-safe and can only be
+// called from a single-thread).
+//
+// Note that Solve(), Update() and SolveWithoutUpdate() are not reentrant so
+// they should not be called concurrently on the same instance of
+// IncrementalSolver.
+//
+// Some solvers may add more restrictions regarding threading. Please see
+// SolverType::kXxx documentation for details.
+class IncrementalSolver {
+ public:
+  IncrementalSolver() = default;
+  IncrementalSolver(const IncrementalSolver&) = delete;
+  IncrementalSolver& operator=(const IncrementalSolver&) = delete;
+  virtual ~IncrementalSolver() = default;
+
+  // Updates the underlying solver with latest model changes and runs the solve.
+  //
+  // A Status error will be returned if inputs are invalid or there is an
+  // unexpected failure in an underlying solver or for some internal math_opt
+  // errors. Otherwise, check SolveResult::termination.reason to see if an
+  // optimal solution was found.
+  //
+  // Memory model: the returned SolveResult owns its own memory (for solutions,
+  // solve stats, etc.), EXPECT for a pointer back to the model. As a result:
+  //  * Keep the model alive to access SolveResult,
+  //  * Avoid unnecessarily copying SolveResult,
+  //  * The result is generally accessible after mutating this, but some care
+  //    is needed if variables or linear constraints are added or deleted.
+  //
+  // See callback.h for documentation on arguments.callback and
+  // arguments.callback_registration.
+  virtual absl::StatusOr<SolveResult> Solve(
+      const SolveArguments& arguments) = 0;
+  absl::StatusOr<SolveResult> Solve() { return Solve({}); }
+
+  // Updates the underlying solver with latest model changes and runs the
+  // computation.
+  //
+  // Same as Solve() but compute the infeasible subsystem.
+  virtual absl::StatusOr<ComputeInfeasibleSubsystemResult>
+  ComputeInfeasibleSubsystem(
+      const ComputeInfeasibleSubsystemArguments& arguments) = 0;
+  absl::StatusOr<ComputeInfeasibleSubsystemResult>
+  ComputeInfeasibleSubsystem() {
+    return ComputeInfeasibleSubsystem({});
+  }
+
+  // Updates the model to solve.
+  //
+  // This is an advanced API, most users should use Solve() above that does the
+  // update and before calling the solver. Calling this function is only useful
+  // for users that want to access to update data or users that need to use
+  // SolveWithoutUpdate() (which should not be common).
+  //
+  // The returned value indicates if the update was possible or if the solver
+  // had to be recreated from scratch (which may happen when the solver does not
+  // support this specific update or any update at all). It also contains the
+  // attempted update data.
+  //
+  // A status error will be returned if the underlying solver has an internal
+  // error.
+  virtual absl::StatusOr<UpdateResult> Update() = 0;
+
+  // Same as Solve() but does not update the underlying solver with the latest
+  // changes to the model.
+  //
+  // This is an advanced API, most users should use Solve().
+  virtual absl::StatusOr<SolveResult> SolveWithoutUpdate(
+      const SolveArguments& arguments) const = 0;
+  absl::StatusOr<SolveResult> SolveWithoutUpdate() const {
+    return SolveWithoutUpdate({});
+  }
+
+  // Same as ComputeInfeasibleSubsystem() but does not update the underlying
+  // solver with the latest changes to the model.
+  //
+  // This is an advanced API, most users should use
+  // ComputeInfeasibleSubsystem().
+  virtual absl::StatusOr<ComputeInfeasibleSubsystemResult>
+  ComputeInfeasibleSubsystemWithoutUpdate(
+      const ComputeInfeasibleSubsystemArguments& arguments) const = 0;
+  absl::StatusOr<ComputeInfeasibleSubsystemResult>
+  ComputeInfeasibleSubsystemWithoutUpdate() const {
+    return ComputeInfeasibleSubsystemWithoutUpdate({});
+  }
+
+  // Returns the underlying solver used.
+  virtual SolverType solver_type() const = 0;
+};
+
+}  // namespace operations_research::math_opt
+
+#endif  // OR_TOOLS_MATH_OPT_CPP_INCREMENTAL_SOLVER_H_
diff --git a/ortools/math_opt/cpp/map_filter.cc b/ortools/math_opt/cpp/map_filter.cc
index 44f17106dc..a321973b45 100644
--- a/ortools/math_opt/cpp/map_filter.cc
+++ b/ortools/math_opt/cpp/map_filter.cc
@@ -63,4 +63,25 @@ absl::StatusOr<MapFilter<LinearConstraint>> LinearConstraintFilterFromProto(
   return result;
 }
 
+absl::StatusOr<MapFilter<QuadraticConstraint>>
+QuadraticConstraintFilterFromProto(const Model& model,
+                                   const SparseVectorFilterProto& proto) {
+  MapFilter<QuadraticConstraint> result = {.skip_zero_values =
+                                               proto.skip_zero_values()};
+  if (proto.filter_by_ids()) {
+    absl::flat_hash_set<QuadraticConstraint> filtered;
+    for (const int64_t id : proto.filtered_ids()) {
+      if (!model.has_quadratic_constraint(id)) {
+        return util::InvalidArgumentErrorBuilder()
+               << "cannot create MapFilter<QuadraticConstraint> from proto, "
+                  "quadratic constraint id: "
+               << id << " not in model";
+      }
+      filtered.insert(model.quadratic_constraint(id));
+    }
+    result.filtered_keys = std::move(filtered);
+  }
+  return result;
+}
+
 }  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/cpp/map_filter.h b/ortools/math_opt/cpp/map_filter.h
index ae7f372e30..13c348bcf9 100644
--- a/ortools/math_opt/cpp/map_filter.h
+++ b/ortools/math_opt/cpp/map_filter.h
@@ -17,11 +17,12 @@
 #ifndef OR_TOOLS_MATH_OPT_CPP_MAP_FILTER_H_
 #define OR_TOOLS_MATH_OPT_CPP_MAP_FILTER_H_
 
-#include <algorithm>
 #include <initializer_list>
 #include <optional>
 
 #include "absl/algorithm/container.h"
+#include "absl/container/flat_hash_set.h"
+#include "absl/status/status.h"
 #include "absl/status/statusor.h"
 #include "ortools/base/status_macros.h"
 #include "ortools/math_opt/cpp/key_types.h"
@@ -31,8 +32,7 @@
 #include "ortools/math_opt/sparse_containers.pb.h"
 #include "ortools/math_opt/storage/model_storage.h"
 
-namespace operations_research {
-namespace math_opt {
+namespace operations_research::math_opt {
 
 // A filter that only keeps some specific key-value pairs of a map.
 //
@@ -123,6 +123,14 @@ absl::StatusOr<MapFilter<Variable>> VariableFilterFromProto(
 absl::StatusOr<MapFilter<LinearConstraint>> LinearConstraintFilterFromProto(
     const Model& model, const SparseVectorFilterProto& proto);
 
+// Returns the MapFilter<QuadraticConstraint> equivalent to `proto`.
+//
+// Requires that (or returns a status error):
+//  * proto.filtered_ids has elements that are quadratic constraints in `model`.
+absl::StatusOr<MapFilter<QuadraticConstraint>>
+QuadraticConstraintFilterFromProto(const Model& model,
+                                   const SparseVectorFilterProto& proto);
+
 // Returns a filter that skips all key-value pairs.
 //
 // This is typically used to disable the dual data in SolveResult when these are
@@ -213,7 +221,6 @@ SparseVectorFilterProto MapFilter<KeyType>::Proto() const {
   return ret;
 }
 
-}  // namespace math_opt
-}  // namespace operations_research
+}  // namespace operations_research::math_opt
 
 #endif  // OR_TOOLS_MATH_OPT_CPP_MAP_FILTER_H_
diff --git a/ortools/math_opt/cpp/matchers.cc b/ortools/math_opt/cpp/matchers.cc
index 519e8bf78f..a6632c02da 100644
--- a/ortools/math_opt/cpp/matchers.cc
+++ b/ortools/math_opt/cpp/matchers.cc
@@ -92,6 +92,8 @@ void PrintTo(const PrimalSolution& primal_solution, std::ostream* const os) {
 
 void PrintTo(const DualSolution& dual_solution, std::ostream* const os) {
   *os << "{dual_values: " << Print(dual_solution.dual_values)
+      << ", quadratic_dual_values: "
+      << Print(dual_solution.quadratic_dual_values)
       << ", reduced_costs: " << Print(dual_solution.reduced_costs)
       << ", objective_value: " << Print(dual_solution.objective_value)
       << ", feasibility_status: " << Print(dual_solution.feasibility_status)
@@ -232,6 +234,23 @@ Matcher<LinearConstraintMap<double>> IsNear(
                                                /*all_keys=*/true, tolerance));
 }
 
+Matcher<absl::flat_hash_map<QuadraticConstraint, double>> IsNear(
+    absl::flat_hash_map<QuadraticConstraint, double> expected,
+    const double tolerance) {
+  return Matcher<absl::flat_hash_map<QuadraticConstraint, double>>(
+      new MapToDoubleMatcher<QuadraticConstraint>(
+          std::move(expected), /*all_keys=*/true, tolerance));
+}
+
+Matcher<absl::flat_hash_map<QuadraticConstraint, double>> IsNearlySubsetOf(
+    absl::flat_hash_map<QuadraticConstraint, double> expected,
+    double tolerance) {
+  return Matcher<absl::flat_hash_map<QuadraticConstraint, double>>(
+      new MapToDoubleMatcher<QuadraticConstraint>(std::move(expected),
+                                                  /*all_keys=*/false,
+                                                  tolerance));
+}
+
 template <typename K>
 Matcher<absl::flat_hash_map<K, double>> IsNear(
     absl::flat_hash_map<K, double> expected, const double tolerance) {
@@ -407,6 +426,8 @@ Matcher<DualSolution> IsNear(DualSolution expected, const double tolerance,
   return AllOf(
       Field("dual_values", &DualSolution::dual_values,
             IsNear(expected.dual_values, tolerance)),
+      Field("quadratic_dual_values", &DualSolution::quadratic_dual_values,
+            IsNear(expected.quadratic_dual_values, tolerance)),
       Field("reduced_costs", &DualSolution::reduced_costs,
             IsNear(expected.reduced_costs, tolerance)),
       Field("objective_value", &DualSolution::objective_value,
@@ -587,8 +608,16 @@ Matcher<SolveResult> TerminatesWith(const TerminationReason expected) {
 }
 
 namespace {
-testing::Matcher<SolveResult> LimitIs(const Limit expected,
-                                      const bool allow_limit_undetermined) {
+
+// Returns a matcher matching only Termination.limit.
+//
+// Note that this is different from LimitIs() which tests both the
+// Termination.limit and the Termination.reason.
+//
+// It matches if either the limit is the expected one, or if it is kUndetermined
+// and when allow_limit_undetermined is true.
+testing::Matcher<SolveResult> TerminationLimitIs(
+    const Limit expected, const bool allow_limit_undetermined) {
   if (allow_limit_undetermined) {
     return Field("termination", &SolveResult::termination,
                  Field("limit", &Termination::limit,
@@ -603,7 +632,7 @@ testing::Matcher<SolveResult> LimitIs(const Limit expected,
 testing::Matcher<SolveResult> TerminatesWithLimit(
     const Limit expected, const bool allow_limit_undetermined) {
   std::vector<Matcher<SolveResult>> matchers;
-  matchers.push_back(LimitIs(expected, allow_limit_undetermined));
+  matchers.push_back(TerminationLimitIs(expected, allow_limit_undetermined));
   matchers.push_back(TerminatesWithOneOf(
       {TerminationReason::kFeasible, TerminationReason::kNoSolutionFound}));
   return AllOfArray(matchers);
@@ -612,7 +641,7 @@ testing::Matcher<SolveResult> TerminatesWithLimit(
 testing::Matcher<SolveResult> TerminatesWithReasonFeasible(
     const Limit expected, const bool allow_limit_undetermined) {
   std::vector<Matcher<SolveResult>> matchers;
-  matchers.push_back(LimitIs(expected, allow_limit_undetermined));
+  matchers.push_back(TerminationLimitIs(expected, allow_limit_undetermined));
   matchers.push_back(TerminatesWith(TerminationReason::kFeasible));
   return AllOfArray(matchers);
 }
@@ -620,7 +649,7 @@ testing::Matcher<SolveResult> TerminatesWithReasonFeasible(
 testing::Matcher<SolveResult> TerminatesWithReasonNoSolutionFound(
     const Limit expected, const bool allow_limit_undetermined) {
   std::vector<Matcher<SolveResult>> matchers;
-  matchers.push_back(LimitIs(expected, allow_limit_undetermined));
+  matchers.push_back(TerminationLimitIs(expected, allow_limit_undetermined));
   matchers.push_back(TerminatesWith(TerminationReason::kNoSolutionFound));
   return AllOfArray(matchers);
 }
@@ -697,14 +726,15 @@ Matcher<Termination> TerminationIsOptimal() {
                              .primal_or_dual_infeasible = false})));
 }
 
-Matcher<Termination> TerminationIsOptimal(const double primal_objective_value,
-                                          const double dual_objective_value,
-                                          const double tolerance) {
+Matcher<Termination> TerminationIsOptimal(
+    const double primal_objective_value,
+    const std::optional<double> dual_objective_value, const double tolerance) {
   return AllOf(
       TerminationIsOptimal(),
       Field("objective_bounds", &Termination::objective_bounds,
             ObjectiveBoundsNear({.primal_bound = primal_objective_value,
-                                 .dual_bound = dual_objective_value},
+                                 .dual_bound = dual_objective_value.value_or(
+                                     primal_objective_value)},
                                 tolerance)));
 }
 
@@ -756,6 +786,24 @@ Matcher<SolveResult> IsOptimalWithDualSolution(
           tolerance));
 }
 
+Matcher<SolveResult> IsOptimalWithDualSolution(
+    const double expected_objective,
+    const LinearConstraintMap<double> expected_dual_values,
+    const absl::flat_hash_map<QuadraticConstraint, double>
+        expected_quadratic_dual_values,
+    const VariableMap<double> expected_reduced_costs, const double tolerance) {
+  return AllOf(
+      IsOptimal(std::make_optional(expected_objective), tolerance),
+      HasDualSolution(
+          DualSolution{
+              .dual_values = expected_dual_values,
+              .quadratic_dual_values = expected_quadratic_dual_values,
+              .reduced_costs = expected_reduced_costs,
+              .objective_value = std::make_optional(expected_objective),
+              .feasibility_status = SolutionStatus::kFeasible},
+          tolerance));
+}
+
 Matcher<SolveResult> HasSolution(PrimalSolution expected,
                                  const double tolerance) {
   return Field(
diff --git a/ortools/math_opt/cpp/matchers.h b/ortools/math_opt/cpp/matchers.h
index 95dd854e47..c99076feb3 100644
--- a/ortools/math_opt/cpp/matchers.h
+++ b/ortools/math_opt/cpp/matchers.h
@@ -101,6 +101,7 @@
 #include "absl/status/statusor.h"
 #include "gtest/gtest.h"
 #include "ortools/base/gmock.h"
+#include "ortools/math_opt/constraints/quadratic/quadratic_constraint.h"
 #include "ortools/math_opt/cpp/linear_constraint.h"
 #include "ortools/math_opt/cpp/math_opt.h"
 #include "ortools/math_opt/cpp/update_result.h"
@@ -142,25 +143,19 @@ testing::Matcher<LinearConstraintMap<double>> IsNearlySubsetOf(
     LinearConstraintMap<double> expected,
     double tolerance = kMatcherDefaultTolerance);
 
-// Checks that the maps have identical keys and values within tolerance. Works
-// for VariableMap, LinearConstraintMap, among other realizations of
-// absl::flat_hash_map. This factory will CHECK-fail if expected contains any
-// NaN values, and any NaN values in the expression compared against will result
-// in the matcher failing.
-template <typename K>
-testing::Matcher<absl::flat_hash_map<K, double>> IsNear(
-    absl::flat_hash_map<K, double> expected,
+// Checks that the maps have identical keys and values within tolerance. This
+// factory will CHECK-fail if expected contains any NaN values.
+testing::Matcher<absl::flat_hash_map<QuadraticConstraint, double>> IsNear(
+    absl::flat_hash_map<QuadraticConstraint, double> expected,
     double tolerance = kMatcherDefaultTolerance);
 
 // Checks that the keys of actual are a subset of the keys of expected, and that
-// for all shared keys, the values are within tolerance. Works for VariableMap,
-// LinearConstraintMap, among other realizations of absl::flat_hash_map. This
-// factory will CHECK-fail if expected contains any NaN values, and any NaN
-// values in the expression compared against will result in the matcher failing.
-template <typename K>
-testing::Matcher<absl::flat_hash_map<K, double>> IsNearlySubsetOf(
-    absl::flat_hash_map<K, double> expected,
-    double tolerance = kMatcherDefaultTolerance);
+// for all shared keys, the values are within tolerance. This factory will
+// CHECK-fail if expected contains any NaN values, and any NaN values in the
+// expression compared against will result in the matcher failing.
+testing::Matcher<absl::flat_hash_map<QuadraticConstraint, double>>
+IsNearlySubsetOf(absl::flat_hash_map<QuadraticConstraint, double> expected,
+                 double tolerance = kMatcherDefaultTolerance);
 
 ////////////////////////////////////////////////////////////////////////////////
 // Matchers for various Variable expressions (e.g. LinearExpression)
@@ -262,10 +257,21 @@ testing::Matcher<Termination> TerminationIsOptimal();
 // Checks that the termination reason is optimal, that the objective bounds
 // match the provided ones and that problem status is consistent with
 // optimality (i.e. status is primal and dual feasible).
+//
+// If dual_objective_value is not set, use the primal_objective_value.
 testing::Matcher<Termination> TerminationIsOptimal(
-    double primal_objective_value, double dual_objective_value,
+    double primal_objective_value,
+    std::optional<double> dual_objective_value = std::nullopt,
     double tolerance = kMatcherDefaultTolerance);
 
+// Checks the following:
+//  * The reason is kFeasible or kNoSolutionFound.
+//  * The limit is the expected one.
+//  * If the optional `detail_matcher` is provided, it matches the `detail`.
+//
+// See TerminatesWithLimit() for a similar matcher for SolveResult.
+//
+// TODO: b/343234961 - Maybe change the name and sync the signatures/features.
 testing::Matcher<Termination> LimitIs(
     Limit limit, testing::Matcher<std::string> detail_matcher = testing::_);
 
@@ -303,6 +309,13 @@ testing::Matcher<SolveResult> IsOptimalWithDualSolution(
     VariableMap<double> expected_reduced_costs,
     double tolerance = kMatcherDefaultTolerance);
 
+testing::Matcher<SolveResult> IsOptimalWithDualSolution(
+    double expected_objective, LinearConstraintMap<double> expected_dual_values,
+    absl::flat_hash_map<QuadraticConstraint, double>
+        expected_quadratic_dual_values,
+    VariableMap<double> expected_reduced_costs,
+    double tolerance = kMatcherDefaultTolerance);
+
 // Checks the following:
 //  * The result has the expected termination reason.
 testing::Matcher<SolveResult> TerminatesWith(TerminationReason expected);
@@ -314,6 +327,8 @@ testing::Matcher<SolveResult> TerminatesWithOneOf(
 // Checks the following:
 //  * The result has termination reason kFeasible or kNoSolutionFound.
 //  * The limit is expected, or is kUndetermined if allow_limit_undetermined.
+//
+// See LimitIs() for a matcher for Termination only.
 testing::Matcher<SolveResult> TerminatesWithLimit(
     Limit expected, bool allow_limit_undetermined = false);
 
diff --git a/ortools/math_opt/cpp/math_opt.h b/ortools/math_opt/cpp/math_opt.h
index 76f8e8af97..d470f2028c 100644
--- a/ortools/math_opt/cpp/math_opt.h
+++ b/ortools/math_opt/cpp/math_opt.h
@@ -16,7 +16,8 @@
 #ifndef OR_TOOLS_MATH_OPT_CPP_MATH_OPT_H_
 #define OR_TOOLS_MATH_OPT_CPP_MATH_OPT_H_
 
-#include "ortools/math_opt/cpp/model.h"  // IWYU pragma: export
-#include "ortools/math_opt/cpp/solve.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/model.h"             // IWYU pragma: export
+#include "ortools/math_opt/cpp/solve.h"             // IWYU pragma: export
+#include "ortools/math_opt/cpp/solver_resources.h"  // IWYU pragma: export
 
 #endif  // OR_TOOLS_MATH_OPT_CPP_MATH_OPT_H_
diff --git a/ortools/math_opt/cpp/model_solve_parameters.cc b/ortools/math_opt/cpp/model_solve_parameters.cc
index c7f8fd5128..4554dad921 100644
--- a/ortools/math_opt/cpp/model_solve_parameters.cc
+++ b/ortools/math_opt/cpp/model_solve_parameters.cc
@@ -25,6 +25,7 @@
 #include "google/protobuf/repeated_field.h"
 #include "ortools/base/status_macros.h"
 #include "ortools/math_opt/cpp/linear_constraint.h"
+#include "ortools/math_opt/cpp/model.h"
 #include "ortools/math_opt/cpp/solution.h"
 #include "ortools/math_opt/cpp/sparse_containers.h"
 #include "ortools/math_opt/cpp/variable_and_expressions.h"
@@ -42,6 +43,8 @@ using ::google::protobuf::RepeatedField;
 ModelSolveParameters ModelSolveParameters::OnlyPrimalVariables() {
   ModelSolveParameters parameters;
   parameters.dual_values_filter = MakeSkipAllFilter<LinearConstraint>();
+  parameters.quadratic_dual_values_filter =
+      MakeSkipAllFilter<QuadraticConstraint>();
   parameters.reduced_costs_filter = MakeSkipAllFilter<Variable>();
   return parameters;
 }
@@ -65,6 +68,9 @@ absl::Status ModelSolveParameters::CheckModelStorage(
       << "invalid variable_values_filter";
   RETURN_IF_ERROR(dual_values_filter.CheckModelStorage(expected_storage))
       << "invalid dual_values_filter";
+  RETURN_IF_ERROR(
+      quadratic_dual_values_filter.CheckModelStorage(expected_storage))
+      << "invalid quadratic_dual_values_filter";
   RETURN_IF_ERROR(reduced_costs_filter.CheckModelStorage(expected_storage))
       << "invalid reduced_costs_filter";
   for (const auto [var, unused] : branching_priorities) {
@@ -165,6 +171,8 @@ ModelSolveParametersProto ModelSolveParameters::Proto() const {
   ModelSolveParametersProto ret;
   *ret.mutable_variable_values_filter() = variable_values_filter.Proto();
   *ret.mutable_dual_values_filter() = dual_values_filter.Proto();
+  *ret.mutable_quadratic_dual_values_filter() =
+      quadratic_dual_values_filter.Proto();
   *ret.mutable_reduced_costs_filter() = reduced_costs_filter.Proto();
 
   if (initial_basis.has_value()) {
@@ -219,6 +227,10 @@ absl::StatusOr<ModelSolveParameters> ModelSolveParameters::FromProto(
       result.dual_values_filter,
       LinearConstraintFilterFromProto(model, proto.dual_values_filter()),
       _ << "invalid dual_values_filter");
+  OR_ASSIGN_OR_RETURN3(result.quadratic_dual_values_filter,
+                       QuadraticConstraintFilterFromProto(
+                           model, proto.quadratic_dual_values_filter()),
+                       _ << "invalid quadratic_dual_values_filter");
   OR_ASSIGN_OR_RETURN3(
       result.reduced_costs_filter,
       VariableFilterFromProto(model, proto.reduced_costs_filter()),
diff --git a/ortools/math_opt/cpp/model_solve_parameters.h b/ortools/math_opt/cpp/model_solve_parameters.h
index 8cedf8bba5..88a913adef 100644
--- a/ortools/math_opt/cpp/model_solve_parameters.h
+++ b/ortools/math_opt/cpp/model_solve_parameters.h
@@ -88,6 +88,10 @@ struct ModelSolveParameters {
   // The filter that is applied to dual_values of DualSolution and DualRay.
   MapFilter<LinearConstraint> dual_values_filter;
 
+  // The filter that is applied to quadratic_dual_values of DualSolution and
+  // DualRay.
+  MapFilter<QuadraticConstraint> quadratic_dual_values_filter;
+
   // The filter that is applied to reduced_costs of DualSolution and DualRay.
   MapFilter<Variable> reduced_costs_filter;
 
diff --git a/ortools/math_opt/cpp/objective.h b/ortools/math_opt/cpp/objective.h
index 3de02ca4b4..6ee9df002d 100644
--- a/ortools/math_opt/cpp/objective.h
+++ b/ortools/math_opt/cpp/objective.h
@@ -74,6 +74,12 @@ class Objective {
   // Returns the constant offset of the objective.
   inline double offset() const;
 
+  // Returns the number of linear terms in the objective.
+  inline int64_t num_linear_terms() const;
+
+  // Returns the number of quadratic terms in the objective.
+  inline int64_t num_quadratic_terms() const;
+
   // Returns the linear coefficient for the variable in the model.
   inline double coefficient(Variable variable) const;
   // Returns the quadratic coefficient for the pair of variables in the model.
@@ -150,6 +156,14 @@ absl::string_view Objective::name() const {
 
 double Objective::offset() const { return storage_->objective_offset(id_); }
 
+int64_t Objective::num_quadratic_terms() const {
+  return storage_->num_quadratic_objective_terms(id_);
+}
+
+int64_t Objective::num_linear_terms() const {
+  return storage_->num_linear_objective_terms(id_);
+}
+
 double Objective::coefficient(const Variable variable) const {
   CHECK_EQ(variable.storage(), storage_)
       << internal::kObjectsFromOtherModelStorage;
diff --git a/ortools/math_opt/cpp/solution.cc b/ortools/math_opt/cpp/solution.cc
index a065fd1c9b..0aa54bf95f 100644
--- a/ortools/math_opt/cpp/solution.cc
+++ b/ortools/math_opt/cpp/solution.cc
@@ -17,20 +17,18 @@
 #include <utility>
 
 #include "absl/container/flat_hash_map.h"
+#include "absl/log/check.h"
+#include "absl/status/status.h"
 #include "absl/strings/string_view.h"
 #include "absl/types/span.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/status_builder.h"
 #include "ortools/base/status_macros.h"
-#include "ortools/math_opt/core/sparse_vector_view.h"
-#include "ortools/math_opt/cpp/linear_constraint.h"
 #include "ortools/math_opt/cpp/sparse_containers.h"
 #include "ortools/math_opt/cpp/variable_and_expressions.h"
 #include "ortools/math_opt/solution.pb.h"
 #include "ortools/math_opt/sparse_containers.pb.h"
 #include "ortools/math_opt/storage/model_storage.h"
-#include "ortools/math_opt/validators/ids_validator.h"
-#include "ortools/math_opt/validators/sparse_vector_validator.h"
 #include "ortools/util/status_macros.h"
 
 namespace operations_research {
@@ -128,6 +126,10 @@ absl::StatusOr<DualSolution> DualSolution::FromProto(
       dual_solution.dual_values,
       LinearConstraintValuesFromProto(model, dual_solution_proto.dual_values()),
       _ << "invalid dual_values");
+  OR_ASSIGN_OR_RETURN3(dual_solution.quadratic_dual_values,
+                       QuadraticConstraintValuesFromProto(
+                           model, dual_solution_proto.quadratic_dual_values()),
+                       _ << "invalid quadratic_dual_values");
   OR_ASSIGN_OR_RETURN3(
       dual_solution.reduced_costs,
       VariableValuesFromProto(model, dual_solution_proto.reduced_costs()),
@@ -147,6 +149,8 @@ absl::StatusOr<DualSolution> DualSolution::FromProto(
 DualSolutionProto DualSolution::Proto() const {
   DualSolutionProto result;
   *result.mutable_dual_values() = LinearConstraintValuesToProto(dual_values);
+  *result.mutable_quadratic_dual_values() =
+      QuadraticConstraintValuesToProto(quadratic_dual_values);
   *result.mutable_reduced_costs() = VariableValuesToProto(reduced_costs);
   if (objective_value.has_value()) {
     result.set_objective_value(*objective_value);
diff --git a/ortools/math_opt/cpp/solution.h b/ortools/math_opt/cpp/solution.h
index 8a40823928..c188e5f170 100644
--- a/ortools/math_opt/cpp/solution.h
+++ b/ortools/math_opt/cpp/solution.h
@@ -19,8 +19,10 @@
 
 #include <optional>
 
+#include "absl/container/flat_hash_map.h"
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
+#include "ortools/math_opt/constraints/quadratic/quadratic_constraint.h"
 #include "ortools/math_opt/cpp/basis_status.h"
 #include "ortools/math_opt/cpp/enums.h"  // IWYU pragma: export
 #include "ortools/math_opt/cpp/linear_constraint.h"
@@ -143,6 +145,11 @@ struct DualSolution {
   DualSolutionProto Proto() const;
 
   LinearConstraintMap<double> dual_values;
+
+  // Note: Some solvers only return quadratic constraint duals when a
+  // solver-specific parameter is set
+  // (see go/mathopt-qcqp-dual#solver-specific).
+  absl::flat_hash_map<QuadraticConstraint, double> quadratic_dual_values;
   VariableMap<double> reduced_costs;
   std::optional<double> objective_value;
 
diff --git a/ortools/math_opt/cpp/solve.cc b/ortools/math_opt/cpp/solve.cc
index a357e05afc..aa5e2be3c5 100644
--- a/ortools/math_opt/cpp/solve.cc
+++ b/ortools/math_opt/cpp/solve.cc
@@ -13,33 +13,23 @@
 
 #include "ortools/math_opt/cpp/solve.h"
 
-#include <functional>
 #include <memory>
-#include <optional>
 #include <utility>
 
-#include "absl/memory/memory.h"
-#include "absl/status/status.h"
 #include "absl/status/statusor.h"
-#include "absl/synchronization/mutex.h"
-#include "ortools/base/status_macros.h"
 #include "ortools/math_opt/callback.pb.h"
+#include "ortools/math_opt/core/base_solver.h"
 #include "ortools/math_opt/core/solver.h"
-#include "ortools/math_opt/cpp/callback.h"
 #include "ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h"
 #include "ortools/math_opt/cpp/compute_infeasible_subsystem_result.h"
-#include "ortools/math_opt/cpp/enums.h"
 #include "ortools/math_opt/cpp/model.h"
-#include "ortools/math_opt/cpp/model_solve_parameters.h"
 #include "ortools/math_opt/cpp/parameters.h"
 #include "ortools/math_opt/cpp/solve_arguments.h"
+#include "ortools/math_opt/cpp/solve_impl.h"
 #include "ortools/math_opt/cpp/solve_result.h"
 #include "ortools/math_opt/cpp/solver_init_arguments.h"
 #include "ortools/math_opt/cpp/streamable_solver_init_arguments.h"
-#include "ortools/math_opt/cpp/update_tracker.h"
 #include "ortools/math_opt/infeasible_subsystem.pb.h"
-#include "ortools/math_opt/storage/model_storage.h"
-#include "ortools/util/status_macros.h"
 
 namespace operations_research {
 namespace math_opt {
@@ -53,55 +43,18 @@ Solver::InitArgs ToSolverInitArgs(const SolverInitArguments& arguments) {
   };
 }
 
-absl::StatusOr<SolveResult> CallSolve(
-    Solver& solver, const ModelStorage* const expected_storage,
-    const SolveArguments& arguments) {
-  RETURN_IF_ERROR(arguments.CheckModelStorageAndCallback(expected_storage));
-
-  Solver::Callback cb = nullptr;
-  absl::Mutex mutex;
-  absl::Status cb_status;  // Guarded by `mutex`.
-  if (arguments.callback != nullptr) {
-    cb = [&](const CallbackDataProto& callback_data_proto) {
-      const CallbackData data(expected_storage, callback_data_proto);
-      // TODO(b/249995436): offer a way for user-callback to return a Status as
-      // well and somehow label it as "user error" (annotation + payload?).
-      const CallbackResult result = arguments.callback(data);
-
-      if (const absl::Status status =
-              result.CheckModelStorage(expected_storage);
-          !status.ok()) {
-        // Note that we use util::StatusBuilder() here as util::Annotate() is
-        // not available in open-source code.
-        util::StatusBuilder builder(status);
-        builder << "invalid CallbackResult returned by user callback";
-
-        const absl::MutexLock lock(&mutex);
-        cb_status.Update(builder);
-
-        // Trigger early termination of the solve.
-        CallbackResultProto result_proto;
-        result_proto.set_terminate(true);
-        return result_proto;
-      }
-
-      return result.Proto();
-    };
-  }
-  ASSIGN_OR_RETURN(
-      const SolveResultProto solve_result,
-      solver.Solve(
-          {.parameters = arguments.parameters.Proto(),
-           .model_parameters = arguments.model_parameters.Proto(),
-           .message_callback = arguments.message_callback,
-           .callback_registration = arguments.callback_registration.Proto(),
-           .user_cb = std::move(cb),
-           .interrupter = arguments.interrupter}));
-
-  const absl::MutexLock lock(&mutex);
-  RETURN_IF_ERROR(cb_status);
-
-  return SolveResult::FromProto(expected_storage, solve_result);
+internal::BaseSolverFactory FactoryFromInitArguments(
+    const SolverInitArguments& arguments) {
+  return [arguments](const SolverTypeProto solver_type, ModelProto model,
+                     SolveInterrupter* const local_canceller)
+             -> absl::StatusOr<std::unique_ptr<BaseSolver>> {
+    // We don't use the local_canceller as in-process solve can't be
+    // cancelled. If an error happens in the callback, the solve_impl code will
+    // use CallbackResultProto::set_terminate() to trigger a cooperative
+    // interruption.
+    return Solver::New(solver_type, std::move(model),
+                       ToSolverInitArgs(arguments));
+  };
 }
 
 }  // namespace
@@ -110,91 +63,28 @@ absl::StatusOr<SolveResult> Solve(const Model& model,
                                   const SolverType solver_type,
                                   const SolveArguments& solve_args,
                                   const SolverInitArguments& init_args) {
-  ASSIGN_OR_RETURN(const std::unique_ptr<Solver> solver,
-                   Solver::New(EnumToProto(solver_type),
-                               model.ExportModel(init_args.remove_names),
-                               ToSolverInitArgs(init_args)));
-  return CallSolve(*solver, model.storage(), solve_args);
+  return internal::SolveImpl(FactoryFromInitArguments(init_args), model,
+                             solver_type, solve_args,
+                             /*user_canceller=*/nullptr,
+                             /*remove_names=*/init_args.remove_names);
 }
 
 absl::StatusOr<ComputeInfeasibleSubsystemResult> ComputeInfeasibleSubsystem(
     const Model& model, const SolverType solver_type,
-    const ComputeInfeasibleSubsystemArguments& infeasible_subsystem_args,
+    const ComputeInfeasibleSubsystemArguments& compute_args,
     const SolverInitArguments& init_args) {
-  ASSIGN_OR_RETURN(
-      const ComputeInfeasibleSubsystemResultProto result_proto,
-      Solver::NonIncrementalComputeInfeasibleSubsystem(
-          model.ExportModel(init_args.remove_names), EnumToProto(solver_type),
-          ToSolverInitArgs(init_args),
-          {.parameters = infeasible_subsystem_args.parameters.Proto(),
-           .message_callback = infeasible_subsystem_args.message_callback,
-           .interrupter = infeasible_subsystem_args.interrupter}));
-  return ComputeInfeasibleSubsystemResult::FromProto(model.storage(),
-                                                     result_proto);
+  return internal::ComputeInfeasibleSubsystemImpl(
+      FactoryFromInitArguments(init_args), model, solver_type, compute_args,
+      /*user_canceller=*/nullptr,
+      /*remove_names=*/init_args.remove_names);
 }
 
-absl::StatusOr<std::unique_ptr<IncrementalSolver>> IncrementalSolver::New(
-    Model* const model, const SolverType solver_type,
-    SolverInitArguments arguments) {
-  if (model == nullptr) {
-    return absl::InvalidArgumentError("input model can't be null");
-  }
-  std::unique_ptr<UpdateTracker> update_tracker = model->NewUpdateTracker();
-  ASSIGN_OR_RETURN(const ModelProto model_proto,
-                   update_tracker->ExportModel(arguments.remove_names));
-  ASSIGN_OR_RETURN(std::unique_ptr<Solver> solver,
-                   Solver::New(EnumToProto(solver_type), model_proto,
-                               ToSolverInitArgs(arguments)));
-  return absl::WrapUnique<IncrementalSolver>(
-      new IncrementalSolver(solver_type, std::move(arguments), model->storage(),
-                            std::move(update_tracker), std::move(solver)));
-}
-
-IncrementalSolver::IncrementalSolver(
-    SolverType solver_type, SolverInitArguments init_args,
-    const ModelStorage* const expected_storage,
-    std::unique_ptr<UpdateTracker> update_tracker,
-    std::unique_ptr<Solver> solver)
-    : solver_type_(solver_type),
-      init_args_(std::move(init_args)),
-      expected_storage_(expected_storage),
-      update_tracker_(std::move(update_tracker)),
-      solver_(std::move(solver)) {}
-
-absl::StatusOr<SolveResult> IncrementalSolver::Solve(
-    const SolveArguments& arguments) {
-  RETURN_IF_ERROR(Update().status());
-  return SolveWithoutUpdate(arguments);
-}
-
-absl::StatusOr<UpdateResult> IncrementalSolver::Update() {
-  ASSIGN_OR_RETURN(std::optional<ModelUpdateProto> model_update,
-                   update_tracker_->ExportModelUpdate(init_args_.remove_names));
-  if (!model_update) {
-    return UpdateResult(true, std::move(model_update));
-  }
-
-  OR_ASSIGN_OR_RETURN3(const bool did_update, solver_->Update(*model_update),
-                       _ << "update failed");
-  RETURN_IF_ERROR(update_tracker_->AdvanceCheckpoint());
-
-  if (did_update) {
-    return UpdateResult(true, std::move(model_update));
-  }
-
-  ASSIGN_OR_RETURN(const ModelProto model_proto,
-                   update_tracker_->ExportModel(init_args_.remove_names));
-  OR_ASSIGN_OR_RETURN3(solver_,
-                       Solver::New(EnumToProto(solver_type_), model_proto,
-                                   ToSolverInitArgs(init_args_)),
-                       _ << "solver re-creation failed");
-
-  return UpdateResult(false, std::move(model_update));
-}
-
-absl::StatusOr<SolveResult> IncrementalSolver::SolveWithoutUpdate(
-    const SolveArguments& arguments) const {
-  return CallSolve(*solver_, expected_storage_, arguments);
+absl::StatusOr<std::unique_ptr<IncrementalSolver>> NewIncrementalSolver(
+    Model* model, SolverType solver_type, SolverInitArguments arguments) {
+  return internal::IncrementalSolverImpl::New(
+      FactoryFromInitArguments(arguments), model, solver_type,
+      /*user_canceller=*/nullptr,
+      /*remove_names=*/arguments.remove_names);
 }
 
 }  // namespace math_opt
diff --git a/ortools/math_opt/cpp/solve.h b/ortools/math_opt/cpp/solve.h
index 648787c989..35b18ab4da 100644
--- a/ortools/math_opt/cpp/solve.h
+++ b/ortools/math_opt/cpp/solve.h
@@ -28,18 +28,16 @@
 #include <memory>
 
 #include "absl/status/statusor.h"
-#include "ortools/math_opt/core/solver.h"
 #include "ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h"  // IWYU pragma: export
 #include "ortools/math_opt/cpp/compute_infeasible_subsystem_result.h"  // IWYU pragma: export
+#include "ortools/math_opt/cpp/incremental_solver.h"  // IWYU pragma: export
 #include "ortools/math_opt/cpp/model.h"
 #include "ortools/math_opt/cpp/parameters.h"             // IWYU pragma: export
 #include "ortools/math_opt/cpp/solve_arguments.h"        // IWYU pragma: export
 #include "ortools/math_opt/cpp/solve_result.h"           // IWYU pragma: export
 #include "ortools/math_opt/cpp/solver_init_arguments.h"  // IWYU pragma: export
 #include "ortools/math_opt/cpp/update_result.h"          // IWYU pragma: export
-#include "ortools/math_opt/cpp/update_tracker.h"         // IWYU pragma: export
 #include "ortools/math_opt/parameters.pb.h"              // IWYU pragma: export
-#include "ortools/math_opt/storage/model_storage.h"
 
 namespace operations_research {
 namespace math_opt {
@@ -106,136 +104,20 @@ using SolveFunction =
 // Thread-safety: this method is safe to call concurrently on the same Model.
 absl::StatusOr<ComputeInfeasibleSubsystemResult> ComputeInfeasibleSubsystem(
     const Model& model, SolverType solver_type,
-    const ComputeInfeasibleSubsystemArguments& infeasible_subsystem_args = {},
+    const ComputeInfeasibleSubsystemArguments& compute_args = {},
     const SolverInitArguments& init_args = {});
 
-// Incremental solve of a model.
+// Creates a new incremental solve for the given model. It may returns an
+// error if the parameters are invalid (for example if the selected solver is
+// not linked in the binary).
 //
-// This is a feature for advance users. Most users should only use the Solve()
-// function above.
-//
-// Here incremental means that the we try to reuse the existing underlying
-// solver internals between each solve. There is no guarantee though that the
-// solver supports all possible model changes. Hence there is not guarantee that
-// performances will be improved when using this class; this is solver
-// dependent. Typically LPs have more to gain from incremental solve than
-// MIPs. In both cases, even if the solver supports the model changes,
-// incremental solve may actually be slower.
-//
-// The New() function instantiates the solver, setup it from the current state
-// of the Model and register on it to listen to changes. Calling Solve() will
-// update the underlying solver with latest model changes and solve this model.
-//
-// Usage:
-//   Model model = ...;
-//   ASSIGN_OR_RETURN(
-//     const std::unique_ptr<IncrementalSolver> incremental_solve,
-//     IncrementalSolver::New(&model, SolverType::kXxx));
-//
-//   ASSIGN_OR_RETURN(const SolveResult result1, incremental_solve->Solve());
-//
-//   model.AddVariable(...);
-//   ...
-//
-//   ASSIGN_OR_RETURN(const SolveResult result2, incremental_solve->Solve());
-//
-//   ...
-//
-// Lifecycle: The IncrementalSolver is only keeping a std::weak_ref on Model's
-// internal data and thus it returns an error if Update() or Solve() are called
-// after the Model has been destroyed. It is fine to destroy the
-// IncrementalSolver after the associated Model though.
-//
-// Thread-safety: The New(), Solve() and Update() methods must not be called
-// while modifying the Model() (adding variables...). The user is expected to
-// use proper synchronization primitives to serialize changes to the model and
-// the use of this object. Note though that it is safe to call methods from
-// different IncrementalSolver instances on the same Model concurrently. Same
-// for calling IncrementalSolver::New(). The destructor is thread-safe and can
-// be called even during a modification of the Model.
-//
-// There is no problem calling SolveWithoutUpdate() concurrently on different
-// instances of IncrementalSolver or while the model is being modified (unless
-// of course the underlying solver itself is not thread-safe and can only be
-// called from a single-thread).
-//
-// Note that Solve(), Update() and SolveWithoutUpdate() are not reentrant so
-// they should not be called concurrently on the same instance of
-// IncrementalSolver.
-//
-// Some solvers may add more restrictions regarding threading. Please see
-// SolverType::kXxx documentation for details.
-class IncrementalSolver {
- public:
-  // Creates a new incremental solve for the given model. It may returns an
-  // error if the parameters are invalid (for example if the selected solver is
-  // not linked in the binary).
-  //
-  // The returned IncrementalSolver keeps a copy of `arguments`. Thus the
-  // content of arguments.non_streamable (for example pointers to solver
-  // specific struct) must be valid until the destruction of the
-  // IncrementalSolver. It also registers on the Model to keep track of updates
-  // (see class documentation for details).
-  static absl::StatusOr<std::unique_ptr<IncrementalSolver>> New(
-      Model* model, SolverType solver_type, SolverInitArguments arguments = {});
-
-  // Updates the underlying solver with latest model changes and runs the solve.
-  //
-  // A Status error will be returned if inputs are invalid or there is an
-  // unexpected failure in an underlying solver or for some internal math_opt
-  // errors. Otherwise, check SolveResult::termination.reason to see if an
-  // optimal solution was found.
-  //
-  // Memory model: the returned SolveResult owns its own memory (for solutions,
-  // solve stats, etc.), EXPECT for a pointer back to the model. As a result:
-  //  * Keep the model alive to access SolveResult,
-  //  * Avoid unnecessarily copying SolveResult,
-  //  * The result is generally accessible after mutating this, but some care
-  //    is needed if variables or linear constraints are added or deleted.
-  //
-  // See callback.h for documentation on arguments.callback and
-  // arguments.callback_registration.
-  absl::StatusOr<SolveResult> Solve(const SolveArguments& arguments = {});
-
-  // Updates the model to solve.
-  //
-  // This is an advanced API, most users should use Solve() above that does the
-  // update and before calling the solver. Calling this function is only useful
-  // for users that want to access to update data or users that need to use
-  // SolveWithoutUpdate() (which should not be common).
-  //
-  // The returned value indicates if the update was possible or if the solver
-  // had to be recreated from scratch (which may happen when the solver does not
-  // support this specific update or any update at all). It also contains the
-  // attempted update data.
-  //
-  // A status error will be returned if the underlying solver has an internal
-  // error.
-  absl::StatusOr<UpdateResult> Update();
-
-  // Same as Solve() but does not update the underlying solver with the latest
-  // changes to the model.
-  //
-  // This is an advanced API, most users should use Solve().
-  absl::StatusOr<SolveResult> SolveWithoutUpdate(
-      const SolveArguments& arguments = {}) const;
-
-  SolverType solver_type() const { return solver_type_; }
-
-  // TODO(b/273961536): Add ComputeInfeasibleSubsystem() member function.
-
- private:
-  IncrementalSolver(SolverType solver_type, SolverInitArguments init_args,
-                    const ModelStorage* expected_storage,
-                    std::unique_ptr<UpdateTracker> update_tracker,
-                    std::unique_ptr<Solver> solver);
-
-  const SolverType solver_type_;
-  const SolverInitArguments init_args_;
-  const ModelStorage* const expected_storage_;
-  const std::unique_ptr<UpdateTracker> update_tracker_;
-  std::unique_ptr<Solver> solver_;
-};
+// The returned IncrementalSolver keeps a copy of `arguments`. Thus the
+// content of arguments.non_streamable (for example pointers to solver
+// specific struct) must be valid until the destruction of the
+// IncrementalSolver. It also registers on the Model to keep track of updates
+// (see class documentation for details).
+absl::StatusOr<std::unique_ptr<IncrementalSolver>> NewIncrementalSolver(
+    Model* model, SolverType solver_type, SolverInitArguments arguments = {});
 
 }  // namespace math_opt
 }  // namespace operations_research
diff --git a/ortools/math_opt/cpp/solve_impl.cc b/ortools/math_opt/cpp/solve_impl.cc
new file mode 100644
index 0000000000..322bb10497
--- /dev/null
+++ b/ortools/math_opt/cpp/solve_impl.cc
@@ -0,0 +1,250 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/cpp/solve_impl.h"
+
+#include <memory>
+#include <optional>
+#include <utility>
+
+#include "absl/functional/any_invocable.h"
+#include "absl/log/check.h"
+#include "absl/memory/memory.h"
+#include "absl/status/status.h"
+#include "absl/status/statusor.h"
+#include "absl/synchronization/mutex.h"
+#include "ortools/base/status_macros.h"
+#include "ortools/math_opt/core/base_solver.h"
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h"
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_result.h"
+#include "ortools/math_opt/cpp/model.h"
+#include "ortools/math_opt/cpp/parameters.h"
+#include "ortools/math_opt/cpp/solve_arguments.h"
+#include "ortools/math_opt/cpp/solve_result.h"
+#include "ortools/math_opt/cpp/update_result.h"
+#include "ortools/math_opt/cpp/update_tracker.h"
+#include "ortools/math_opt/storage/model_storage.h"
+#include "ortools/util/solve_interrupter.h"
+#include "ortools/util/status_macros.h"
+
+namespace operations_research::math_opt::internal {
+namespace {
+
+absl::StatusOr<SolveResult> CallSolve(
+    BaseSolver& solver, const ModelStorage* const expected_storage,
+    const SolveArguments& arguments, SolveInterrupter& local_canceller) {
+  RETURN_IF_ERROR(arguments.CheckModelStorageAndCallback(expected_storage));
+
+  BaseSolver::Callback cb = nullptr;
+  absl::Mutex mutex;
+  absl::Status cb_status;  // Guarded by `mutex`.
+  if (arguments.callback != nullptr) {
+    cb = [&](const CallbackDataProto& callback_data_proto) {
+      const CallbackData data(expected_storage, callback_data_proto);
+      const CallbackResult result = arguments.callback(data);
+      if (const absl::Status status =
+              result.CheckModelStorage(expected_storage);
+          !status.ok()) {
+        // Note that we use util::StatusBuilder() here as util::Annotate() is
+        // not available in open-source code.
+        util::StatusBuilder builder(status);
+        builder << "invalid CallbackResult returned by user callback";
+
+        {  // Limit `lock` scope.
+          const absl::MutexLock lock(&mutex);
+          cb_status.Update(builder);
+        }
+
+        // Trigger subprocess cancellation.
+        local_canceller.Interrupt();
+
+        // Trigger early termination of the solve if cancellation is not
+        // supported (i.e. in in-process solve).
+        CallbackResultProto result_proto;
+        result_proto.set_terminate(true);
+        return result_proto;
+      }
+      return result.Proto();
+    };
+  }
+
+  const absl::StatusOr<SolveResultProto> solve_result_proto = solver.Solve(
+      {.parameters = arguments.parameters.Proto(),
+       .model_parameters = arguments.model_parameters.Proto(),
+       .message_callback = arguments.message_callback,
+       .callback_registration = arguments.callback_registration.Proto(),
+       .user_cb = std::move(cb),
+       .interrupter = arguments.interrupter});
+
+  // solver.Solve() returns an error on cancellation by local_canceller but in
+  // that case we want to ignore this error and return status generated in the
+  // callback instead.
+  {  // Limit `lock` scope.
+    const absl::MutexLock lock(&mutex);
+    RETURN_IF_ERROR(cb_status);
+  }
+
+  if (!solve_result_proto.ok()) {
+    return solve_result_proto.status();
+  }
+
+  return SolveResult::FromProto(expected_storage, solve_result_proto.value());
+}
+
+absl::StatusOr<ComputeInfeasibleSubsystemResult> CallComputeInfeasibleSubsystem(
+    BaseSolver& solver, const ModelStorage* const expected_storage,
+    const ComputeInfeasibleSubsystemArguments& arguments,
+    SolveInterrupter& local_canceller) {
+  ASSIGN_OR_RETURN(
+      const ComputeInfeasibleSubsystemResultProto compute_result_proto,
+      solver.ComputeInfeasibleSubsystem(
+          {.parameters = arguments.parameters.Proto(),
+           .message_callback = arguments.message_callback,
+           .interrupter = arguments.interrupter}));
+
+  return ComputeInfeasibleSubsystemResult::FromProto(expected_storage,
+                                                     compute_result_proto);
+}
+
+}  // namespace
+
+absl::StatusOr<SolveResult> SolveImpl(
+    const BaseSolverFactory solver_factory, const Model& model,
+    const SolverType solver_type, const SolveArguments& solve_args,
+    const SolveInterrupter* const user_canceller, const bool remove_names) {
+  SolveInterrupter local_canceller;
+  const ScopedSolveInterrupterCallback user_canceller_cb(
+      user_canceller, [&]() { local_canceller.Interrupt(); });
+  ASSIGN_OR_RETURN(
+      const std::unique_ptr<BaseSolver> solver,
+      solver_factory(EnumToProto(solver_type), model.ExportModel(remove_names),
+                     &local_canceller));
+  return CallSolve(*solver, model.storage(), solve_args, local_canceller);
+}
+
+absl::StatusOr<ComputeInfeasibleSubsystemResult> ComputeInfeasibleSubsystemImpl(
+    const BaseSolverFactory solver_factory, const Model& model,
+    const SolverType solver_type,
+    const ComputeInfeasibleSubsystemArguments& compute_args,
+    const SolveInterrupter* const user_canceller, const bool remove_names) {
+  SolveInterrupter local_canceller;
+  const ScopedSolveInterrupterCallback user_canceller_cb(
+      user_canceller, [&]() { local_canceller.Interrupt(); });
+  ASSIGN_OR_RETURN(
+      const std::unique_ptr<BaseSolver> subprocess_solver,
+      solver_factory(EnumToProto(solver_type), model.ExportModel(remove_names),
+                     &local_canceller));
+  return CallComputeInfeasibleSubsystem(*subprocess_solver, model.storage(),
+                                        compute_args, local_canceller);
+}
+
+absl::StatusOr<std::unique_ptr<IncrementalSolverImpl>>
+IncrementalSolverImpl::New(BaseSolverFactory solver_factory, Model* const model,
+                           const SolverType solver_type,
+                           const SolveInterrupter* const user_canceller,
+                           const bool remove_names) {
+  if (model == nullptr) {
+    return absl::InvalidArgumentError("input model can't be null");
+  }
+  auto local_canceller = std::make_shared<SolveInterrupter>();
+  auto user_canceller_cb =
+      std::make_unique<const ScopedSolveInterrupterCallback>(
+          user_canceller,
+          [local_canceller]() { local_canceller->Interrupt(); });
+  std::unique_ptr<UpdateTracker> update_tracker = model->NewUpdateTracker();
+  ASSIGN_OR_RETURN(ModelProto model_proto,
+                   update_tracker->ExportModel(remove_names));
+  ASSIGN_OR_RETURN(
+      std::unique_ptr<BaseSolver> solver,
+      solver_factory(EnumToProto(solver_type), std::move(model_proto),
+                     local_canceller.get()));
+  return absl::WrapUnique<IncrementalSolverImpl>(new IncrementalSolverImpl(
+      std::move(solver_factory), solver_type, remove_names,
+      std::move(local_canceller), std::move(user_canceller_cb),
+      model->storage(), std::move(update_tracker), std::move(solver)));
+}
+
+IncrementalSolverImpl::IncrementalSolverImpl(
+    BaseSolverFactory solver_factory, SolverType solver_type,
+    const bool remove_names, std::shared_ptr<SolveInterrupter> local_canceller,
+    std::unique_ptr<const ScopedSolveInterrupterCallback> user_canceller_cb,
+    const ModelStorage* const expected_storage,
+    std::unique_ptr<UpdateTracker> update_tracker,
+    std::unique_ptr<BaseSolver> solver)
+    : solver_factory_(std::move(solver_factory)),
+      solver_type_(solver_type),
+      remove_names_(remove_names),
+      local_canceller_(std::move(local_canceller)),
+      user_canceller_cb_(std::move(user_canceller_cb)),
+      expected_storage_(expected_storage),
+      update_tracker_(std::move(update_tracker)),
+      solver_(std::move(solver)) {}
+
+absl::StatusOr<SolveResult> IncrementalSolverImpl::Solve(
+    const SolveArguments& arguments) {
+  // TODO: b/260337466 - Add permanent errors and concurrency protection.
+  RETURN_IF_ERROR(Update().status());
+  return SolveWithoutUpdate(arguments);
+}
+
+absl::StatusOr<ComputeInfeasibleSubsystemResult>
+IncrementalSolverImpl::ComputeInfeasibleSubsystem(
+    const ComputeInfeasibleSubsystemArguments& arguments) {
+  // TODO: b/260337466 - Add permanent errors and concurrency protection.
+  RETURN_IF_ERROR(Update().status());
+  return ComputeInfeasibleSubsystemWithoutUpdate(arguments);
+}
+
+absl::StatusOr<UpdateResult> IncrementalSolverImpl::Update() {
+  // TODO: b/260337466 - Add permanent errors and concurrency protection.
+  ASSIGN_OR_RETURN(std::optional<ModelUpdateProto> model_update,
+                   update_tracker_->ExportModelUpdate(remove_names_));
+  if (!model_update.has_value()) {
+    return UpdateResult(true);
+  }
+
+  OR_ASSIGN_OR_RETURN3(const bool did_update,
+                       solver_->Update(*std::move(model_update)),
+                       _ << "update failed");
+  RETURN_IF_ERROR(update_tracker_->AdvanceCheckpoint());
+
+  if (did_update) {
+    return UpdateResult(true);
+  }
+
+  ASSIGN_OR_RETURN(ModelProto model_proto,
+                   update_tracker_->ExportModel(remove_names_));
+  OR_ASSIGN_OR_RETURN3(
+      solver_,
+      solver_factory_(EnumToProto(solver_type_), std::move(model_proto),
+                      local_canceller_.get()),
+      _ << "solver re-creation failed");
+
+  return UpdateResult(false);
+}
+
+absl::StatusOr<SolveResult> IncrementalSolverImpl::SolveWithoutUpdate(
+    const SolveArguments& arguments) const {
+  // TODO: b/260337466 - Add permanent errors and concurrency protection.
+  return CallSolve(*solver_, expected_storage_, arguments, *local_canceller_);
+}
+
+absl::StatusOr<ComputeInfeasibleSubsystemResult>
+IncrementalSolverImpl::ComputeInfeasibleSubsystemWithoutUpdate(
+    const ComputeInfeasibleSubsystemArguments& arguments) const {
+  // TODO: b/260337466 - Add permanent errors and concurrency protection.
+  return CallComputeInfeasibleSubsystem(*solver_, expected_storage_, arguments,
+                                        *local_canceller_);
+}
+
+}  // namespace operations_research::math_opt::internal
diff --git a/ortools/math_opt/cpp/solve_impl.h b/ortools/math_opt/cpp/solve_impl.h
new file mode 100644
index 0000000000..8782b0523b
--- /dev/null
+++ b/ortools/math_opt/cpp/solve_impl.h
@@ -0,0 +1,124 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_CPP_SOLVE_IMPL_H_
+#define OR_TOOLS_MATH_OPT_CPP_SOLVE_IMPL_H_
+
+#include <memory>
+
+#include "absl/functional/any_invocable.h"
+#include "absl/status/statusor.h"
+#include "ortools/math_opt/core/base_solver.h"
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_arguments.h"
+#include "ortools/math_opt/cpp/compute_infeasible_subsystem_result.h"
+#include "ortools/math_opt/cpp/incremental_solver.h"
+#include "ortools/math_opt/cpp/model.h"
+#include "ortools/math_opt/cpp/parameters.h"
+#include "ortools/math_opt/cpp/solve_arguments.h"
+#include "ortools/math_opt/cpp/solve_result.h"
+#include "ortools/math_opt/cpp/update_result.h"
+#include "ortools/math_opt/cpp/update_tracker.h"
+#include "ortools/math_opt/storage/model_storage.h"
+#include "ortools/util/solve_interrupter.h"
+
+namespace operations_research::math_opt::internal {
+
+// A factory of solver.
+//
+// The local_canceller is a local interrupter that exists in the scope of
+// SolveImpl(), ComputeInfeasibleSubsystemImpl() or IncrementalSolverImpl. It is
+// triggered:
+// * either when the user_canceller is triggered
+// * or when the BaseSolver::Callback returns an invalid CallbackResultProto; in
+//   that case a new CallbackResultProto with its `terminate` set to true is
+//   also returned instead.
+//
+// Solvers that don't support cancellation (i.e. in-process solving) should
+// ignore the local_canceller: this use case won't have a user_canceller and the
+// CallbackResultProto.terminate will terminate the solve as soon as possible if
+// the CallbackResultProto is invalid.
+using BaseSolverFactory =
+    absl::AnyInvocable<absl::StatusOr<std::unique_ptr<BaseSolver>>(
+        SolverTypeProto solver_type, ModelProto model,
+        SolveInterrupter* local_canceller) const>;
+
+// Solves the input model.
+//
+// The `user_canceller` parameter is optional.
+absl::StatusOr<SolveResult> SolveImpl(BaseSolverFactory solver_factory,
+                                      const Model& model,
+                                      SolverType solver_type,
+                                      const SolveArguments& solve_args,
+                                      const SolveInterrupter* user_canceller,
+                                      bool remove_names);
+
+// ComputeInfeasibleSubsystems the input model in a subprocess.
+//
+// The `user_canceller` parameter is optional.
+absl::StatusOr<ComputeInfeasibleSubsystemResult> ComputeInfeasibleSubsystemImpl(
+    BaseSolverFactory solver_factory, const Model& model,
+    SolverType solver_type,
+    const ComputeInfeasibleSubsystemArguments& compute_args,
+    const SolveInterrupter* user_canceller, bool remove_names);
+
+// Incremental solve of a model.
+class IncrementalSolverImpl : public IncrementalSolver {
+ public:
+  // Creates a new incremental solve.
+  //
+  // The `user_canceller` parameter is optional.
+  static absl::StatusOr<std::unique_ptr<IncrementalSolverImpl>> New(
+      BaseSolverFactory solver_factory, Model* model, SolverType solver_type,
+      const SolveInterrupter* user_canceller, bool remove_names);
+
+  absl::StatusOr<SolveResult> Solve(const SolveArguments& arguments) override;
+
+  absl::StatusOr<ComputeInfeasibleSubsystemResult> ComputeInfeasibleSubsystem(
+      const ComputeInfeasibleSubsystemArguments& arguments) override;
+
+  absl::StatusOr<UpdateResult> Update() override;
+
+  absl::StatusOr<SolveResult> SolveWithoutUpdate(
+      const SolveArguments& arguments) const override;
+
+  absl::StatusOr<ComputeInfeasibleSubsystemResult>
+  ComputeInfeasibleSubsystemWithoutUpdate(
+      const ComputeInfeasibleSubsystemArguments& arguments) const override;
+
+  SolverType solver_type() const override { return solver_type_; }
+
+ private:
+  IncrementalSolverImpl(
+      BaseSolverFactory solver_factory, SolverType solver_type,
+      bool remove_names, std::shared_ptr<SolveInterrupter> local_canceller,
+      std::unique_ptr<const ScopedSolveInterrupterCallback> user_canceller_cb,
+      const ModelStorage* expected_storage,
+      std::unique_ptr<UpdateTracker> update_tracker,
+      std::unique_ptr<BaseSolver> solver);
+
+  const BaseSolverFactory solver_factory_;
+  const SolverType solver_type_;
+  const bool remove_names_;
+  // Here we use a shared_ptr so that we don't have to make sure that
+  // user_canceller_cb_, which points to local_canceller_ via a lambda-capture,
+  // can be destroyed after local_canceller_ without risk.
+  std::shared_ptr<SolveInterrupter> local_canceller_;
+  std::unique_ptr<const ScopedSolveInterrupterCallback> user_canceller_cb_;
+  const ModelStorage* const expected_storage_;
+  const std::unique_ptr<UpdateTracker> update_tracker_;
+  std::unique_ptr<BaseSolver> solver_;
+};
+
+}  // namespace operations_research::math_opt::internal
+
+#endif  // OR_TOOLS_MATH_OPT_CPP_SOLVE_IMPL_H_
diff --git a/ortools/math_opt/cpp/solver_init_arguments.h b/ortools/math_opt/cpp/solver_init_arguments.h
index 975cd28a62..503b1b8be5 100644
--- a/ortools/math_opt/cpp/solver_init_arguments.h
+++ b/ortools/math_opt/cpp/solver_init_arguments.h
@@ -17,14 +17,12 @@
 #ifndef OR_TOOLS_MATH_OPT_CPP_SOLVER_INIT_ARGUMENTS_H_
 #define OR_TOOLS_MATH_OPT_CPP_SOLVER_INIT_ARGUMENTS_H_
 
-#include <memory>
-
 #include "ortools/math_opt/core/non_streamable_solver_init_arguments.h"  // IWYU pragma: export
 #include "ortools/math_opt/cpp/streamable_solver_init_arguments.h"  // IWYU pragma: export
 
 namespace operations_research::math_opt {
 
-// Arguments passed to Solve() and IncrementalSolver::New() to control the
+// Arguments passed to Solve() and NewIncrementalSolver() to control the
 // instantiation of the solver.
 //
 // Usage with streamable arguments:
diff --git a/ortools/math_opt/cpp/solver_resources.h b/ortools/math_opt/cpp/solver_resources.h
index 7779180249..045dce3d07 100644
--- a/ortools/math_opt/cpp/solver_resources.h
+++ b/ortools/math_opt/cpp/solver_resources.h
@@ -11,6 +11,9 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+// IWYU pragma: private, include "ortools/math_opt/cpp/math_opt.h"
+// IWYU pragma: friend "ortools/math_opt/cpp/.*"
+
 #ifndef OR_TOOLS_MATH_OPT_CPP_SOLVER_RESOURCES_H_
 #define OR_TOOLS_MATH_OPT_CPP_SOLVER_RESOURCES_H_
 
diff --git a/ortools/math_opt/cpp/sparse_containers.cc b/ortools/math_opt/cpp/sparse_containers.cc
index 70820c19eb..e3382b9710 100644
--- a/ortools/math_opt/cpp/sparse_containers.cc
+++ b/ortools/math_opt/cpp/sparse_containers.cc
@@ -22,9 +22,11 @@
 #include "absl/container/flat_hash_map.h"
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
+#include "absl/types/span.h"
 #include "google/protobuf/map.h"
 #include "ortools/base/status_builder.h"
 #include "ortools/base/status_macros.h"
+#include "ortools/math_opt/constraints/quadratic/quadratic_constraint.h"
 #include "ortools/math_opt/core/sparse_vector_view.h"
 #include "ortools/math_opt/cpp/basis_status.h"
 #include "ortools/math_opt/cpp/linear_constraint.h"
@@ -32,6 +34,7 @@
 #include "ortools/math_opt/cpp/variable_and_expressions.h"
 #include "ortools/math_opt/storage/model_storage.h"
 #include "ortools/math_opt/storage/model_storage_types.h"
+#include "ortools/math_opt/validators/sparse_vector_validator.h"
 
 namespace operations_research::math_opt {
 namespace {
@@ -123,6 +126,17 @@ absl::Status LinearConstraintIdsExist(const ModelStorage* const model,
   return absl::OkStatus();
 }
 
+absl::Status QuadraticConstraintIdsExist(const ModelStorage* const model,
+                                         const absl::Span<const int64_t> ids) {
+  for (const int64_t id : ids) {
+    if (!model->has_constraint(QuadraticConstraintId(id))) {
+      return util::InvalidArgumentErrorBuilder()
+             << "no quadratic constraint with id " << id << " exists";
+    }
+  }
+  return absl::OkStatus();
+}
+
 }  // namespace
 
 absl::StatusOr<VariableMap<double>> VariableValuesFromProto(
@@ -185,6 +199,21 @@ SparseDoubleVectorProto LinearConstraintValuesToProto(
   return MapToProto(linear_constraint_values);
 }
 
+absl::StatusOr<absl::flat_hash_map<QuadraticConstraint, double>>
+QuadraticConstraintValuesFromProto(
+    const ModelStorage* const model,
+    const SparseDoubleVectorProto& quad_cons_proto) {
+  RETURN_IF_ERROR(CheckSparseVectorProto(quad_cons_proto));
+  RETURN_IF_ERROR(QuadraticConstraintIdsExist(model, quad_cons_proto.ids()));
+  return MakeView(quad_cons_proto).as_map<QuadraticConstraint>(model);
+}
+
+SparseDoubleVectorProto QuadraticConstraintValuesToProto(
+    const absl::flat_hash_map<QuadraticConstraint, double>&
+        quadratic_constraint_values) {
+  return MapToProto(quadratic_constraint_values);
+}
+
 absl::StatusOr<VariableMap<BasisStatus>> VariableBasisFromProto(
     const ModelStorage* const model,
     const SparseBasisStatusVector& basis_proto) {
diff --git a/ortools/math_opt/cpp/sparse_containers.h b/ortools/math_opt/cpp/sparse_containers.h
index 76c2b7a416..6fa634b8e9 100644
--- a/ortools/math_opt/cpp/sparse_containers.h
+++ b/ortools/math_opt/cpp/sparse_containers.h
@@ -20,12 +20,8 @@
 #include <cstdint>
 
 #include "absl/container/flat_hash_map.h"
-#include "absl/strings/string_view.h"
-#include "absl/types/span.h"
 #include "ortools/base/logging.h"
-#include "ortools/base/status_builder.h"
-#include "ortools/base/status_macros.h"
-#include "ortools/math_opt/core/sparse_vector_view.h"
+#include "ortools/math_opt/constraints/quadratic/quadratic_constraint.h"
 #include "ortools/math_opt/cpp/basis_status.h"
 #include "ortools/math_opt/cpp/linear_constraint.h"
 #include "ortools/math_opt/cpp/objective.h"
@@ -33,9 +29,6 @@
 #include "ortools/math_opt/solution.pb.h"
 #include "ortools/math_opt/sparse_containers.pb.h"
 #include "ortools/math_opt/storage/model_storage.h"
-#include "ortools/math_opt/validators/ids_validator.h"
-#include "ortools/math_opt/validators/sparse_vector_validator.h"
-#include "ortools/util/status_macros.h"
 
 namespace operations_research::math_opt {
 
@@ -101,6 +94,26 @@ absl::StatusOr<LinearConstraintMap<double>> LinearConstraintValuesFromProto(
 SparseDoubleVectorProto LinearConstraintValuesToProto(
     const LinearConstraintMap<double>& linear_constraint_values);
 
+// Returns the absl::flat_hash_map<QuadraticConstraint, double> equivalent to
+// `quad_cons_proto`.
+//
+// Requires that (or returns a status error):
+//  * quad_cons_proto.ids and quad_cons_proto.values have equal size.
+//  * quad_cons_proto.ids is sorted.
+//  * quad_cons_proto.ids has elements that are quadratic constraints in `model`
+//    (this implies that each id is in [0, max(int64_t))).
+//
+// Note that the values of quad_cons_proto.values are not checked (it may have
+// NaNs).
+absl::StatusOr<absl::flat_hash_map<QuadraticConstraint, double>>
+QuadraticConstraintValuesFromProto(
+    const ModelStorage* model, const SparseDoubleVectorProto& quad_cons_proto);
+
+// Returns the proto equivalent of quadratic_constraint_values.
+SparseDoubleVectorProto QuadraticConstraintValuesToProto(
+    const absl::flat_hash_map<QuadraticConstraint, double>&
+        quadratic_constraint_values);
+
 // Returns the VariableMap<BasisStatus> equivalent to `basis_proto`.
 //
 // Requires that (or returns a status error):
diff --git a/ortools/math_opt/cpp/update_result.h b/ortools/math_opt/cpp/update_result.h
index 98e6a4bb79..ea3f4bf643 100644
--- a/ortools/math_opt/cpp/update_result.h
+++ b/ortools/math_opt/cpp/update_result.h
@@ -14,26 +14,18 @@
 #ifndef OR_TOOLS_MATH_OPT_CPP_UPDATE_RESULT_H_
 #define OR_TOOLS_MATH_OPT_CPP_UPDATE_RESULT_H_
 
-#include <optional>
-#include <utility>
-
 #include "ortools/math_opt/model_update.pb.h"
 
 namespace operations_research::math_opt {
 
 // Result of the Update() on an incremental solver.
 struct UpdateResult {
-  UpdateResult(const bool did_update, std::optional<ModelUpdateProto> update)
-      : did_update(did_update), update(std::move(update)) {}
+  explicit UpdateResult(const bool did_update) : did_update(did_update) {}
 
   // True if the solver has been successfully updated or if no update was
   // necessary (in which case `update` will be nullopt). False if the solver
   // had to be recreated.
   bool did_update;
-
-  // The update that was attempted on the solver. Can be nullopt when no
-  // update was needed (the model was not changed).
-  std::optional<ModelUpdateProto> update;
 };
 
 }  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/io/BUILD.bazel b/ortools/math_opt/io/BUILD.bazel
index 4600ff6dfe..428beafedd 100644
--- a/ortools/math_opt/io/BUILD.bazel
+++ b/ortools/math_opt/io/BUILD.bazel
@@ -76,3 +76,26 @@ cc_library(
         "@com_google_absl//absl/status:statusor",
     ],
 )
+
+cc_library(
+    name = "lp_parser",
+    srcs = ["lp_parser.cc"],
+    hdrs = ["lp_parser.h"],
+    deps = [
+        ":mps_converter",
+        "//ortools/base",
+        "//ortools/base:file",
+        "//ortools/base:path",
+        "//ortools/base:status_macros",
+        "//ortools/base:temp_path",
+        "//ortools/gscip",
+        "//ortools/linear_solver:scip_helper_macros",
+        "//ortools/math_opt:model_cc_proto",
+        "//ortools/util:status_macros",
+        "@com_google_absl//absl/memory",
+        "@com_google_absl//absl/status",
+        "@com_google_absl//absl/status:statusor",
+        "@com_google_absl//absl/strings",
+        "@scip//:libscip",
+    ],
+)
diff --git a/ortools/math_opt/io/lp/BUILD.bazel b/ortools/math_opt/io/lp/BUILD.bazel
new file mode 100644
index 0000000000..6f912cb4f2
--- /dev/null
+++ b/ortools/math_opt/io/lp/BUILD.bazel
@@ -0,0 +1,52 @@
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+cc_library(
+    name = "lp_model",
+    srcs = ["lp_model.cc"],
+    hdrs = ["lp_model.h"],
+    deps = [
+        ":lp_name",
+        "//ortools/base:intops",
+        "//ortools/base:status_macros",
+        "//ortools/base:strong_vector",
+        "//ortools/util:fp_roundtrip_conv",
+        "@com_google_absl//absl/container:flat_hash_map",
+        "@com_google_absl//absl/status",
+        "@com_google_absl//absl/status:statusor",
+        "@com_google_absl//absl/strings",
+    ],
+)
+
+cc_library(
+    name = "matchers",
+    testonly = 1,
+    srcs = ["matchers.cc"],
+    hdrs = ["matchers.h"],
+    deps = [
+        ":lp_model",
+        "//ortools/base:gmock",
+    ],
+)
+
+cc_library(
+    name = "lp_name",
+    srcs = ["lp_name.cc"],
+    hdrs = ["lp_name.h"],
+    deps = [
+        "//ortools/base:status_macros",
+        "@com_google_absl//absl/status",
+        "@com_google_absl//absl/strings",
+        "@com_google_absl//absl/strings:string_view",
+    ],
+)
diff --git a/ortools/math_opt/io/lp/lp_model.cc b/ortools/math_opt/io/lp/lp_model.cc
new file mode 100644
index 0000000000..c9276d5994
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_model.cc
@@ -0,0 +1,149 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/lp_model.h"
+
+#include <cmath>
+#include <cstdlib>
+#include <ostream>
+#include <string>
+#include <utility>
+
+#include "absl/status/status.h"
+#include "absl/status/statusor.h"
+#include "absl/strings/escaping.h"
+#include "absl/strings/str_cat.h"
+#include "absl/strings/str_join.h"
+#include "ortools/base/status_builder.h"
+#include "ortools/base/status_macros.h"
+#include "ortools/math_opt/io/lp/lp_name.h"
+#include "ortools/util/fp_roundtrip_conv.h"
+
+namespace operations_research::lp_format {
+namespace {
+
+absl::Status ValidateRelation(Relation relation) {
+  switch (relation) {
+    case Relation::kLessOrEqual:
+    case Relation::kEqual:
+    case Relation::kGreaterOrEqual:
+      return absl::OkStatus();
+  }
+  return util::InvalidArgumentErrorBuilder()
+         << "Invalid Relation: " << static_cast<int>(relation);
+}
+
+}  // namespace
+
+std::ostream& operator<<(std::ostream& ostr, const Relation relation) {
+  switch (relation) {
+    case Relation::kGreaterOrEqual:
+      ostr << ">=";
+      return ostr;
+    case Relation::kEqual:
+      ostr << "=";
+      return ostr;
+    case Relation::kLessOrEqual:
+      ostr << "<=";
+      return ostr;
+  }
+  ostr << "__invalid_Relation_" << static_cast<int>(relation) << "__";
+  return ostr;
+}
+
+std::ostream& operator<<(std::ostream& ostr, const Constraint& constraint) {
+  auto term_format = [](std::string* out,
+                        const std::pair<double, VariableIndex>& term) {
+    absl::StrAppend(out, "{", RoundTripDoubleFormat(term.first), ", ",
+                    term.second, "}");
+  };
+  ostr << "terms: {" << absl::StrJoin(constraint.terms, ", ", term_format)
+       << "} relation: " << constraint.relation
+       << " rhs: " << RoundTripDoubleFormat(constraint.rhs) << " name: \""
+       << absl::CEscape(constraint.name) << "\"";
+  return ostr;
+}
+
+absl::StatusOr<ConstraintIndex> LpModel::AddConstraint(Constraint constraint) {
+  if (!constraint.name.empty()) {
+    RETURN_IF_ERROR(ValidateName(constraint.name)) << "invalid constraint name";
+  }
+  if (constraint.terms.empty()) {
+    return absl::InvalidArgumentError("constraint must have at least one term");
+  }
+  for (const auto [coef, var] : constraint.terms) {
+    if (var < VariableIndex{0} || var >= variables_.end_index()) {
+      return util::InvalidArgumentErrorBuilder()
+             << "variable ids should be in [0," << variables_.end_index()
+             << ") but found: " << var;
+    }
+    if (!std::isfinite(coef)) {
+      return util::InvalidArgumentErrorBuilder()
+             << "All coefficients in constraints must be finite and not NaN "
+                "but found: "
+             << coef;
+    }
+  }
+  RETURN_IF_ERROR(ValidateRelation(constraint.relation));
+  if (std::isnan(constraint.rhs)) {
+    return absl::InvalidArgumentError("rhs of constraint was NaN");
+  }
+  ConstraintIndex result = constraints_.end_index();
+  constraints_.push_back(std::move(constraint));
+  return result;
+}
+
+absl::StatusOr<VariableIndex> LpModel::AddVariable(absl::string_view name) {
+  if (variable_names_.contains(name)) {
+    return util::InvalidArgumentErrorBuilder()
+           << "duplicate variable name: " << name;
+  }
+  RETURN_IF_ERROR(ValidateName(name)) << "invalid variable name";
+  const VariableIndex index = variables_.end_index();
+  variables_.push_back(std::string(name));
+  variable_names_[name] = index;
+  return index;
+}
+
+std::ostream& operator<<(std::ostream& ostr, const LpModel& model) {
+  ostr << "SUBJECT TO" << std::endl;
+  for (const Constraint& constraint : model.constraints()) {
+    ostr << "  ";
+    if (!constraint.name.empty()) {
+      ostr << constraint.name << ": ";
+    }
+    bool first = true;
+    for (const auto [coef, var] : constraint.terms) {
+      if (first) {
+        if (coef != 1.0) {
+          ostr << RoundTripDoubleFormat(coef) << " ";
+        }
+      } else {
+        if (coef > 0) {
+          ostr << " + ";
+        } else {
+          ostr << " - ";
+        }
+        ostr << RoundTripDoubleFormat(std::abs(coef)) << " ";
+      }
+      ostr << model.variables()[var];
+      first = false;
+    }
+    ostr << " " << constraint.relation << " "
+         << RoundTripDoubleFormat(constraint.rhs) << std::endl;
+  }
+  ostr << "END" << std::endl;
+  return ostr;
+}
+
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp/lp_model.h b/ortools/math_opt/io/lp/lp_model.h
new file mode 100644
index 0000000000..df7c869c38
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_model.h
@@ -0,0 +1,97 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_IO_LP_LP_MODEL_H_
+#define OR_TOOLS_MATH_OPT_IO_LP_LP_MODEL_H_
+
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "absl/container/flat_hash_map.h"
+#include "ortools/base/strong_int.h"
+#include "ortools/base/strong_vector.h"
+
+namespace operations_research::lp_format {
+
+DEFINE_STRONG_INT_TYPE(VariableIndex, int64_t);
+DEFINE_STRONG_INT_TYPE(ConstraintIndex, int64_t);
+using Term = std::pair<double, VariableIndex>;
+enum class Relation { kLessOrEqual, kGreaterOrEqual, kEqual };
+
+std::ostream& operator<<(std::ostream& ostr, Relation relation);
+
+struct Constraint {
+  std::vector<Term> terms;
+  Relation relation = Relation::kLessOrEqual;
+  double rhs = 0.0;
+  std::string name;
+};
+
+// Note: this prints an exact representation of the data in Constraint, not
+// the string form of the constraint in LP format.
+std::ostream& operator<<(std::ostream& ostr, const Constraint& constraint);
+
+// The contents of an optimization model in LP file format.
+//
+// You can convert this to a string in the LP file format using <<, and read
+// from a string in the LP file format using ParseLp from parse_lp.h.
+class LpModel {
+ public:
+  LpModel() = default;
+
+  // Adds a new variable to the model and returns it. Errors if name:
+  //  * is empty
+  //  * is the same as any existing variable name
+  //  * has invalid characters for the LP file format
+  //
+  // Variable names are case sensitive.
+  absl::StatusOr<VariableIndex> AddVariable(absl::string_view name);
+
+  // Adds a new constraint to the model and returns its index.
+  //
+  // Errors if:
+  //  * a variable id from constraint.bounds is out of bounds
+  //  * constraint.relation is an invalid enum
+  //  * a coefficient on constraint.terms is Inf or NaN
+  //  * the name has invalid characters
+  //  * there are no terms in the constraint
+  //
+  // Constraint names can be repeated but this is not recommended.
+  absl::StatusOr<ConstraintIndex> AddConstraint(Constraint constraint);
+
+  const absl::flat_hash_map<std::string, VariableIndex>& variable_names()
+      const {
+    return variable_names_;
+  }
+  const util_intops::StrongVector<VariableIndex, std::string>& variables()
+      const {
+    return variables_;
+  }
+  const util_intops::StrongVector<ConstraintIndex, Constraint>& constraints()
+      const {
+    return constraints_;
+  }
+
+ private:
+  absl::flat_hash_map<std::string, VariableIndex> variable_names_;
+  util_intops::StrongVector<VariableIndex, std::string> variables_;
+  util_intops::StrongVector<ConstraintIndex, Constraint> constraints_;
+};
+
+// Prints the model in LP format to ostr.
+std::ostream& operator<<(std::ostream& ostr, const LpModel& model);
+
+}  // namespace operations_research::lp_format
+
+#endif  // OR_TOOLS_MATH_OPT_IO_LP_LP_MODEL_H_
diff --git a/ortools/math_opt/io/lp/lp_model_test.cc b/ortools/math_opt/io/lp/lp_model_test.cc
new file mode 100644
index 0000000000..2387ddc4a9
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_model_test.cc
@@ -0,0 +1,191 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/lp_model.h"
+
+#include <limits>
+#include <sstream>
+#include <string>
+
+#include "absl/status/status.h"
+#include "gtest/gtest.h"
+#include "ortools/base/gmock.h"
+#include "ortools/math_opt/io/lp/matchers.h"
+
+namespace operations_research::lp_format {
+namespace {
+
+using ::testing::AllOf;
+using ::testing::ElementsAre;
+using ::testing::HasSubstr;
+using ::testing::IsEmpty;
+using ::testing::Pair;
+using ::testing::StrEq;
+using ::testing::UnorderedElementsAre;
+using ::testing::status::IsOkAndHolds;
+using ::testing::status::StatusIs;
+
+constexpr double kInf = std::numeric_limits<double>::infinity();
+
+TEST(LpModelTest, EmptyModel) {
+  LpModel model;
+  EXPECT_THAT(model.variables(), IsEmpty());
+  EXPECT_THAT(model.variable_names(), IsEmpty());
+  EXPECT_THAT(model.constraints(), IsEmpty());
+}
+
+TEST(LpModelTest, AddVariableSuccess) {
+  LpModel model;
+  EXPECT_THAT(model.AddVariable("x"), IsOkAndHolds(VariableIndex{0}));
+  EXPECT_THAT(model.AddVariable("y78"), IsOkAndHolds(VariableIndex{1}));
+  EXPECT_THAT(model.variables(), ElementsAre("x", "y78"));
+  EXPECT_THAT(model.variable_names(),
+              UnorderedElementsAre(Pair("x", VariableIndex{0}),
+                                   Pair("y78", VariableIndex{1})));
+}
+
+TEST(LpModelTest, AddVariableRepeatNameError) {
+  LpModel model;
+  EXPECT_OK(model.AddVariable("xyz"));
+  EXPECT_THAT(model.AddVariable("xyz"),
+              StatusIs(absl::StatusCode::kInvalidArgument,
+                       AllOf(HasSubstr("duplicate"), HasSubstr("xyz"))));
+}
+
+TEST(LpModelTest, AddVariableInvalidNameError) {
+  LpModel model;
+  EXPECT_THAT(model.AddVariable("4x"),
+              StatusIs(absl::StatusCode::kInvalidArgument,
+                       AllOf(HasSubstr("invalid variable"), HasSubstr("4x"))));
+}
+
+TEST(LpModelTest, AddConstraintSuccess) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {.terms = {{4.0, x}},
+                        .relation = Relation::kEqual,
+                        .rhs = 2.5,
+                        .name = "c"};
+  EXPECT_THAT(model.AddConstraint(c), IsOkAndHolds(ConstraintIndex{0}));
+  EXPECT_THAT(model.constraints(), ElementsAre(ConstraintEquals(c)));
+}
+
+TEST(LpModelTest, AddConstraintNoTermsError) {
+  LpModel model;
+  const Constraint c = {.relation = Relation::kEqual, .rhs = 2.5, .name = "c"};
+  EXPECT_THAT(model.AddConstraint(c),
+              StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("term")));
+}
+
+TEST(LpModelTest, AddConstraintBadRelationError) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {.terms = {{4.0, x}},
+                        .relation = static_cast<Relation>(-12),
+                        .rhs = 2.5,
+                        .name = "c"};
+  EXPECT_THAT(
+      model.AddConstraint(c),
+      StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("Relation")));
+}
+
+TEST(LpModelTest, AddConstraintNanInTermsError) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {
+      .terms = {{std::numeric_limits<double>::quiet_NaN(), x}},
+      .relation = Relation::kEqual,
+      .rhs = 2.5,
+      .name = "c"};
+  EXPECT_THAT(model.AddConstraint(c),
+              StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("NaN")));
+}
+
+TEST(LpModelTest, AddConstraintInfInTermsError) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {.terms = {{kInf, x}},
+                        .relation = Relation::kEqual,
+                        .rhs = 2.5,
+                        .name = "c"};
+  EXPECT_THAT(
+      model.AddConstraint(c),
+      StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("finite")));
+}
+
+TEST(LpModelTest, AddConstraintNanInRhsError) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {.terms = {{2.0, x}},
+                        .relation = Relation::kEqual,
+                        .rhs = std::numeric_limits<double>::quiet_NaN(),
+                        .name = "c"};
+  EXPECT_THAT(model.AddConstraint(c),
+              StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("NaN")));
+}
+
+TEST(LpModelTest, AddConstraintVariableNotInModelError) {
+  LpModel model;
+  const VariableIndex x(0);
+  const Constraint c = {.terms = {{2.0, x}},
+                        .relation = Relation::kEqual,
+                        .rhs = 4.0,
+                        .name = "c"};
+  EXPECT_THAT(
+      model.AddConstraint(c),
+      StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("variable ids")));
+}
+
+template <typename T>
+std::string ToString(const T& t) {
+  std::stringstream ss;
+  ss << t;
+  return ss.str();
+}
+
+TEST(RelationStream, PrintsString) {
+  EXPECT_EQ(ToString(Relation::kEqual), "=");
+  EXPECT_EQ(ToString(Relation::kGreaterOrEqual), ">=");
+  EXPECT_EQ(ToString(Relation::kLessOrEqual), "<=");
+  EXPECT_THAT(ToString(static_cast<Relation>(-1)), HasSubstr("invalid"));
+}
+
+TEST(ConstraintStream, PrintsString) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.5, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  // StrEq gives better output than EXPECT_EQ on failure
+  EXPECT_THAT(
+      ToString(c),
+      StrEq("terms: {{1, 0}, {4.5, 3}} relation: = rhs: 5 name: \"cat\""));
+}
+
+// These tests are intentionally not exhaustive. Better to test the stream
+// operator by round tripping with parsing. We do not have a strong contract
+// on the string produced, and testing against a large string is brittle.
+TEST(LpModelTest, BasicStreamTest) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  const Constraint c = {.terms = {{4.0, x}},
+                        .relation = Relation::kEqual,
+                        .rhs = 2.5,
+                        .name = "c"};
+  ASSERT_OK(model.AddConstraint(c));
+  // StrEq gives better output than EXPECT_EQ on failure
+  EXPECT_THAT(ToString(model), StrEq("SUBJECT TO\n  c: 4 x = 2.5\nEND\n"));
+}
+
+}  // namespace
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp/lp_name.cc b/ortools/math_opt/io/lp/lp_name.cc
new file mode 100644
index 0000000000..ad6fbc9458
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_name.cc
@@ -0,0 +1,71 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/lp_name.h"
+
+#include "absl/status/status.h"
+#include "absl/strings/ascii.h"
+#include "absl/strings/string_view.h"
+#include "ortools/base/status_builder.h"
+
+namespace operations_research::lp_format {
+
+bool ValidateCharInName(unsigned char c, bool leading) {
+  if (absl::ascii_isalpha(c)) {
+    return true;
+  }
+  if (!leading) {
+    if (c == '.' || absl::ascii_isdigit(c)) {
+      return true;
+    }
+  }
+  switch (c) {
+    case '!':
+    case '"':
+    case '#':
+    case '$':
+    case '%':
+    case '&':
+    case '(':
+    case ')':
+    case ',':
+    case ';':
+    case '?':
+    case '@':
+    case '_':
+    case '`':
+    case '\'':
+    case '{':
+    case '}':
+    case '~':
+      return true;
+    default:
+      return false;
+  }
+}
+
+absl::Status ValidateName(absl::string_view name) {
+  if (name.empty()) {
+    return absl::InvalidArgumentError("empty name invalid");
+  }
+  for (int i = 0; i < name.size(); ++i) {
+    if (!ValidateCharInName(name[i], i == 0)) {
+      return util::InvalidArgumentErrorBuilder()
+             << "invalid character: " << name[i] << " at index: " << i
+             << " in: " << name;
+    }
+  }
+  return absl::OkStatus();
+}
+
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp/lp_name.h b/ortools/math_opt/io/lp/lp_name.h
new file mode 100644
index 0000000000..14620800f4
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_name.h
@@ -0,0 +1,31 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_IO_LP_LP_NAME_H_
+#define OR_TOOLS_MATH_OPT_IO_LP_LP_NAME_H_
+
+#include "absl/status/status.h"
+
+namespace operations_research::lp_format {
+
+// Returns true if c is valid character to be included the name of a variable
+// or constraint in an LP file, where is_leading indicates if c is the first
+// character of the name.
+bool ValidateCharInName(unsigned char c, bool is_leading);
+
+// Checks if `name` is a valid name for a variable or constraint in an LP file.
+absl::Status ValidateName(absl::string_view name);
+
+}  // namespace operations_research::lp_format
+
+#endif  // OR_TOOLS_MATH_OPT_IO_LP_LP_NAME_H_
diff --git a/ortools/math_opt/io/lp/lp_name_test.cc b/ortools/math_opt/io/lp/lp_name_test.cc
new file mode 100644
index 0000000000..d28c08577f
--- /dev/null
+++ b/ortools/math_opt/io/lp/lp_name_test.cc
@@ -0,0 +1,67 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/lp_name.h"
+
+#include "absl/status/status.h"
+#include "absl/strings/str_cat.h"
+#include "gtest/gtest.h"
+#include "ortools/base/gmock.h"
+
+namespace operations_research::lp_format {
+namespace {
+
+using ::testing::AllOf;
+using ::testing::HasSubstr;
+using ::testing::status::StatusIs;
+
+TEST(ValidateCharInNameTest, BasicUse) {
+  for (bool is_leading : {false, true}) {
+    SCOPED_TRACE(absl::StrCat("is leading: ", is_leading));
+    for (unsigned char c : {'a', 'A', 'b', 'B', 'z', 'Z', '_', '{', '}'}) {
+      SCOPED_TRACE(absl::StrCat("testing character: ", c));
+      EXPECT_TRUE(ValidateCharInName(c, is_leading));
+    }
+    for (unsigned char c : {'+', '-', '*', '/', ':', '\0'}) {
+      SCOPED_TRACE(absl::StrCat("testing character: ", c));
+      EXPECT_FALSE(ValidateCharInName(c, is_leading));
+    }
+  }
+}
+
+TEST(ValidateCharInNameTest, LeadingChars) {
+  for (bool is_leading : {false, true}) {
+    SCOPED_TRACE(absl::StrCat("is leading: ", is_leading));
+    for (unsigned char c : {'.', '0', '1', '9'}) {
+      SCOPED_TRACE(absl::StrCat("testing character: ", c));
+      const bool should_be_allowed = !is_leading;
+      EXPECT_EQ(ValidateCharInName(c, is_leading), should_be_allowed);
+    }
+  }
+}
+
+TEST(ValidateNameTest, BasicUse) {
+  EXPECT_OK(ValidateName("x8"));
+  EXPECT_OK(ValidateName("A_b_C"));
+  EXPECT_THAT(ValidateName(""),
+              StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("empty")));
+  EXPECT_THAT(ValidateName("8x"), StatusIs(absl::StatusCode::kInvalidArgument,
+                                           AllOf(HasSubstr("index: 0"),
+                                                 HasSubstr("character: 8"))));
+  EXPECT_THAT(ValidateName("x-8"), StatusIs(absl::StatusCode::kInvalidArgument,
+                                            AllOf(HasSubstr("index: 1"),
+                                                  HasSubstr("character: -"))));
+}
+
+}  // namespace
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp/matchers.cc b/ortools/math_opt/io/lp/matchers.cc
new file mode 100644
index 0000000000..6ae55fb822
--- /dev/null
+++ b/ortools/math_opt/io/lp/matchers.cc
@@ -0,0 +1,50 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/matchers.h"
+
+#include <vector>
+
+#include "gtest/gtest.h"
+#include "ortools/base/gmock.h"
+#include "ortools/math_opt/io/lp/lp_model.h"
+
+namespace operations_research::lp_format {
+namespace {
+using ::testing::AllOf;
+using ::testing::ElementsAreArray;
+using ::testing::Field;
+using ::testing::Matcher;
+using ::testing::Property;
+}  // namespace
+
+Matcher<Constraint> ConstraintEquals(const Constraint& expected) {
+  return AllOf(
+      Field("terms", &Constraint::terms, ElementsAreArray(expected.terms)),
+      Field("relation", &Constraint::relation, expected.relation),
+      Field("rhs", &Constraint::rhs, expected.rhs),
+      Field("name", &Constraint::name, expected.name));
+}
+
+Matcher<LpModel> ModelEquals(const LpModel& expected) {
+  std::vector<Matcher<Constraint>> expected_constraints;
+  for (const Constraint& constraint : expected.constraints()) {
+    expected_constraints.push_back(ConstraintEquals(constraint));
+  }
+  return AllOf(Property("variables", &LpModel::variables,
+                        ElementsAreArray(expected.variables())),
+               Property("constraints", &LpModel::constraints,
+                        ElementsAreArray(expected_constraints)));
+}
+
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp/matchers.h b/ortools/math_opt/io/lp/matchers.h
new file mode 100644
index 0000000000..b018122ccb
--- /dev/null
+++ b/ortools/math_opt/io/lp/matchers.h
@@ -0,0 +1,28 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_IO_LP_MATCHERS_H_
+#define OR_TOOLS_MATH_OPT_IO_LP_MATCHERS_H_
+
+#include "ortools/base/gmock.h"
+#include "ortools/math_opt/io/lp/lp_model.h"
+
+namespace operations_research::lp_format {
+
+testing::Matcher<Constraint> ConstraintEquals(const Constraint& expected);
+
+testing::Matcher<LpModel> ModelEquals(const LpModel& expected);
+
+}  // namespace operations_research::lp_format
+
+#endif  // OR_TOOLS_MATH_OPT_IO_LP_MATCHERS_H_
diff --git a/ortools/math_opt/io/lp/matchers_test.cc b/ortools/math_opt/io/lp/matchers_test.cc
new file mode 100644
index 0000000000..d26769fad9
--- /dev/null
+++ b/ortools/math_opt/io/lp/matchers_test.cc
@@ -0,0 +1,123 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp/matchers.h"
+
+#include "gtest/gtest.h"
+#include "ortools/base/gmock.h"
+#include "ortools/math_opt/io/lp/lp_model.h"
+
+namespace operations_research::lp_format {
+namespace {
+
+TEST(ConstraintEqualsTest, Equal) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.0, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  EXPECT_THAT(c, ConstraintEquals(c));
+}
+
+TEST(ConstraintEqualsTest, WrongNameNoMatch) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.0, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  Constraint d = c;
+  d.name = "dog";
+  EXPECT_THAT(c, Not(ConstraintEquals(d)));
+}
+
+TEST(ConstraintEqualsTest, WrongRhsNoMatch) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.0, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  Constraint d = c;
+  d.rhs = 4.0;
+  EXPECT_THAT(c, Not(ConstraintEquals(d)));
+}
+
+TEST(ConstraintEqualsTest, WrongRelationNoMatch) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.0, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  Constraint d = c;
+  d.relation = Relation::kGreaterOrEqual;
+  EXPECT_THAT(c, Not(ConstraintEquals(d)));
+}
+
+TEST(ConstraintEqualsTest, WrongTermsNoMatch) {
+  const Constraint c = {
+      .terms = {{1.0, VariableIndex(0)}, {4.0, VariableIndex(3)}},
+      .relation = Relation::kEqual,
+      .rhs = 5.0,
+      .name = "cat"};
+  Constraint d = c;
+  d.terms.clear();
+  EXPECT_THAT(c, Not(ConstraintEquals(d)));
+}
+
+TEST(ModelEqualsTest, ModelEqualsSelf) {
+  LpModel model;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x, model.AddVariable("x"));
+  ASSERT_OK(model.AddConstraint({.terms = {{2.0, x}},
+                                 .relation = Relation::kLessOrEqual,
+                                 .rhs = 4.0,
+                                 .name = "c"}));
+  EXPECT_THAT(model, ModelEquals(model));
+}
+
+TEST(ModelEqualsTest, EmptyModelsEqual) {
+  LpModel actual;
+  LpModel expected;
+  EXPECT_THAT(actual, ModelEquals(expected));
+}
+
+TEST(ModelEqualsTest, DifferentVariablesNotEqual) {
+  LpModel actual;
+  ASSERT_OK(actual.AddVariable("x"));
+
+  LpModel expected;
+  ASSERT_OK(expected.AddVariable("y"));
+
+  EXPECT_THAT(actual, Not(ModelEquals(expected)));
+}
+
+TEST(ModelEqualsTest, DifferentConstraintsNotEqual) {
+  LpModel actual;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x_actual, actual.AddVariable("x"));
+  ASSERT_OK(actual.AddConstraint({.terms = {{2.0, x_actual}},
+                                  .relation = Relation::kLessOrEqual,
+                                  .rhs = 4.0,
+                                  .name = "c"}));
+
+  LpModel expected;
+  ASSERT_OK_AND_ASSIGN(const VariableIndex x_expected,
+                       expected.AddVariable("x"));
+  // RHS is different
+  ASSERT_OK(actual.AddConstraint({.terms = {{2.0, x_expected}},
+                                  .relation = Relation::kLessOrEqual,
+                                  .rhs = 5.0,
+                                  .name = "c"}));
+
+  EXPECT_THAT(actual, Not(ModelEquals(expected)));
+}
+
+}  // namespace
+}  // namespace operations_research::lp_format
diff --git a/ortools/math_opt/io/lp_converter.h b/ortools/math_opt/io/lp_converter.h
index 830eae736f..99068fbe42 100644
--- a/ortools/math_opt/io/lp_converter.h
+++ b/ortools/math_opt/io/lp_converter.h
@@ -26,6 +26,10 @@ namespace operations_research::math_opt {
 // The RemoveNames() function can be used on the model to remove names if they
 // should not be exported.
 //
+// Warning: LP format does not support "range constraints" (linear constraints
+// with both upper and lower bounds). An equivalent model without range
+// constraints is exported instead.
+//
 // For more information about the different LP file formats:
 // http://lpsolve.sourceforge.net/5.5/lp-format.htm
 // http://lpsolve.sourceforge.net/5.5/CPLEX-format.htm
diff --git a/ortools/math_opt/io/lp_parser.cc b/ortools/math_opt/io/lp_parser.cc
new file mode 100644
index 0000000000..dcc7700fcb
--- /dev/null
+++ b/ortools/math_opt/io/lp_parser.cc
@@ -0,0 +1,75 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/lp_parser.h"
+
+#include <memory>
+#include <string>
+
+#include "absl/memory/memory.h"
+#include "absl/status/status.h"
+#include "absl/status/statusor.h"
+#include "absl/strings/string_view.h"
+#include "ortools/base/helpers.h"
+#include "ortools/base/options.h"
+#include "ortools/base/path.h"
+#include "ortools/base/status_macros.h"
+#include "ortools/base/temp_path.h"
+#include "ortools/gscip/gscip.h"
+#include "ortools/linear_solver/scip_helper_macros.h"
+#include "ortools/math_opt/io/mps_converter.h"
+#include "ortools/math_opt/model.pb.h"
+#include "ortools/util/status_macros.h"
+#include "scip/def.h"
+#include "scip/scip_prob.h"
+
+namespace operations_research::math_opt {
+namespace {
+
+absl::Status ScipConvertLpToMps(const std::string& lp_filename_in,
+                                const std::string& mps_filename_out) {
+  ASSIGN_OR_RETURN(const std::unique_ptr<GScip> gscip, GScip::Create(""));
+  // Warning: puts gscip into an invalid state, but the underlying SCIP is fine.
+  RETURN_IF_SCIP_ERROR(
+      SCIPreadProb(gscip->scip(), lp_filename_in.c_str(), "lp"));
+  RETURN_IF_SCIP_ERROR(SCIPwriteOrigProblem(
+      gscip->scip(), mps_filename_out.c_str(), "mps", /*genericnames=*/FALSE));
+  return absl::OkStatus();
+}
+
+}  // namespace
+
+absl::StatusOr<ModelProto> ModelProtoFromLp(const absl::string_view lp_data) {
+  // Set up temporary files
+  const auto dir = absl::WrapUnique(TempPath::Create(TempPath::Local));
+  if (dir == nullptr) {
+    return absl::InternalError(
+        "creating temporary directory when parsing LP file failed");
+  }
+  const std::string lp_file = file::JoinPath(dir->path(), "model.lp");
+  RETURN_IF_ERROR(file::SetContents(lp_file, lp_data, file::Defaults()));
+  const std::string mps_file = file::JoinPath(dir->path(), "model.mps");
+
+  // Do the conversion
+  RETURN_IF_ERROR(ScipConvertLpToMps(lp_file, mps_file))
+      << "failed to convert LP file with SCIP";
+  ASSIGN_OR_RETURN(const std::string mps_data,
+                   file::GetContents(mps_file, file::Defaults()));
+  OR_ASSIGN_OR_RETURN3(
+      ModelProto result, MpsToModelProto(mps_data),
+      _ << "failed to parse MPS (produced by SCIP from LP file)");
+  result.clear_name();
+  return result;
+}
+
+}  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/io/lp_parser.h b/ortools/math_opt/io/lp_parser.h
new file mode 100644
index 0000000000..4c90e9bf83
--- /dev/null
+++ b/ortools/math_opt/io/lp_parser.h
@@ -0,0 +1,63 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Parses a the model in "CPLEX LP" format using SCIP.
+//
+// Note that ../../lp_data/lp_parser.h parses ".lp" files in the LPSolve
+// version of the LP format, which is different from the (now) more standard
+// CPLEX version of the LP format. These formats are not compatible. See
+//  https://lpsolve.sourceforge.net/5.5/lp-format.htm
+//  https://lpsolve.sourceforge.net/5.5/CPLEX-format.htm
+// for a comparison.
+#ifndef OR_TOOLS_MATH_OPT_IO_LP_PARSER_H_
+#define OR_TOOLS_MATH_OPT_IO_LP_PARSER_H_
+
+#include "absl/status/statusor.h"
+#include "absl/strings/string_view.h"
+#include "ortools/math_opt/model.pb.h"
+
+namespace operations_research::math_opt {
+
+// Parses a the model in "CPLEX LP" format.
+//
+// This function creates and destroys local temporary files and thus is not
+// portable.
+//
+// For large models, this will not work on diskless jobs in prod.
+//
+// Warnings:
+//  * Only a linear objective and linear constraints are supported. When SCIP is
+//    used, indicator constraints are also supported.
+//  * The names of indicator constraints are not preserved when using SCIP
+//  * The variables may be permuted.
+//  * Two sided constraints are not in the LP format. If you round trip a Model
+//    Proto with lp_converter.h, the two sided constraints are and rewritten as
+//    two one sided constraints with new names.
+//
+// OR-Tools does not have an LP file parser, so we go from LP file to SCIP, then
+// export to MPS, parse the MPS to ModelProto. This is not efficient, but
+// usually still much faster than solving an LP or MIP. Note the SCIP LP parser
+// actually supports SOS and quadratics, but the OR-tools MPS reader does not.
+//
+// It would be preferable to write an LP Parser from scratch and delete this.
+//
+// For more information about the different LP file formats:
+// http://lpsolve.sourceforge.net/5.5/lp-format.htm
+// http://lpsolve.sourceforge.net/5.5/CPLEX-format.htm
+// https://www.ibm.com/docs/en/icos/12.8.0.0?topic=cplex-lp-file-format-algebraic-representation
+// http://www.gurobi.com/documentation/5.1/reference-manual/node871
+absl::StatusOr<ModelProto> ModelProtoFromLp(absl::string_view lp_data);
+
+}  // namespace operations_research::math_opt
+
+#endif  // OR_TOOLS_MATH_OPT_IO_LP_PARSER_H_
diff --git a/ortools/math_opt/io/mps_converter.cc b/ortools/math_opt/io/mps_converter.cc
index af89c82401..d31141c834 100644
--- a/ortools/math_opt/io/mps_converter.cc
+++ b/ortools/math_opt/io/mps_converter.cc
@@ -40,4 +40,10 @@ absl::StatusOr<ModelProto> ReadMpsFile(const absl::string_view filename) {
   return MPModelProtoToMathOptModel(mp_model);
 }
 
+absl::StatusOr<ModelProto> MpsToModelProto(absl::string_view mps_data) {
+  ASSIGN_OR_RETURN(const MPModelProto mp_model,
+                   glop::MpsDataToMPModelProto(mps_data));
+  return MPModelProtoToMathOptModel(mp_model);
+}
+
 }  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/io/mps_converter.h b/ortools/math_opt/io/mps_converter.h
index 228e8e7d8e..7fdbdc395b 100644
--- a/ortools/math_opt/io/mps_converter.h
+++ b/ortools/math_opt/io/mps_converter.h
@@ -28,6 +28,9 @@ namespace operations_research::math_opt {
 // should not be exported.
 absl::StatusOr<std::string> ModelProtoToMps(const ModelProto& model);
 
+// Converts a string with the contents of an MPS file into a ModelProto.
+absl::StatusOr<ModelProto> MpsToModelProto(absl::string_view mps_data);
+
 // Reads an MPS file and converts it to a ModelProto (like MpsToModelProto
 // above, but takes a file name instead of the file contents and reads the file.
 //
diff --git a/ortools/math_opt/io/python/BUILD.bazel b/ortools/math_opt/io/python/BUILD.bazel
new file mode 100644
index 0000000000..b0ebb9d16b
--- /dev/null
+++ b/ortools/math_opt/io/python/BUILD.bazel
@@ -0,0 +1,32 @@
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+load("@pybind11_bazel//:build_defs.bzl", "pybind_extension")
+
+package(default_visibility = ["//visibility:public"])
+
+pybind_extension(
+    name = "mps_converter",
+    srcs = ["mps_converter.cc"],
+    deps =
+        [
+            "//ortools/math_opt:model_cc_proto",
+            "//ortools/math_opt:model_parameters_cc_proto",
+            "//ortools/math_opt:model_update_cc_proto",
+            "//ortools/math_opt:result_cc_proto",
+            "//ortools/math_opt/io:mps_converter",
+            "@org_pybind11_abseil//pybind11_abseil:import_status_module",
+            "@org_pybind11_abseil//pybind11_abseil:status_casters",
+            "@pybind11_protobuf//pybind11_protobuf:native_proto_caster",
+        ],
+)
diff --git a/ortools/math_opt/io/python/mps_converter.cc b/ortools/math_opt/io/python/mps_converter.cc
new file mode 100644
index 0000000000..01311499c8
--- /dev/null
+++ b/ortools/math_opt/io/python/mps_converter.cc
@@ -0,0 +1,40 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/math_opt/io/mps_converter.h"
+
+#include <pybind11/functional.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "ortools/math_opt/model.pb.h"
+#include "ortools/math_opt/model_parameters.pb.h"
+#include "ortools/math_opt/model_update.pb.h"
+#include "ortools/math_opt/result.pb.h"
+#include "pybind11/attr.h"
+#include "pybind11/cast.h"
+#include "pybind11/gil.h"
+#include "pybind11_abseil/import_status_module.h"
+#include "pybind11_abseil/status_casters.h"  // IWYU pragma: keep
+#include "pybind11_protobuf/native_proto_caster.h"
+
+namespace operations_research {
+namespace math_opt {
+
+PYBIND11_MODULE(mps_converter, m) {
+  m.def("model_proto_to_mps", &ModelProtoToMps);
+  m.def("mps_to_model_proto", &MpsToModelProto);
+}
+
+}  // namespace math_opt
+}  // namespace operations_research
diff --git a/ortools/math_opt/io/python/mps_converter_test.py b/ortools/math_opt/io/python/mps_converter_test.py
new file mode 100644
index 0000000000..a9231057ac
--- /dev/null
+++ b/ortools/math_opt/io/python/mps_converter_test.py
@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from absl.testing import absltest
+from ortools.math_opt import model_pb2
+from ortools.math_opt.io.python import mps_converter
+from ortools.math_opt.python import mathopt
+from ortools.math_opt.python.testing import compare_proto
+
+_MODEL_MPS = """* Generated by MPModelProtoExporter
+*   Name             : unbounded_integers
+*   Format           : Free
+*   Constraints      : 1
+*   Variables        : 2
+*     Binary         : 0
+*     Integer        : 2
+*     Continuous     : 0
+NAME          unbounded_integers
+ROWS
+ N  COST
+ G  c
+COLUMNS
+  INTSTART  'MARKER'                            'INTORG'
+    x       COST         1  c            1
+    y       COST         1  c            1
+  INTEND    'MARKER'                            'INTEND'
+RHS
+    RHS     c            2
+BOUNDS
+ LI BOUND   x            1
+ LI BOUND   y            4
+ENDATA
+"""
+
+
+def simple_model_proto() -> model_pb2.ModelProto:
+    model = mathopt.Model(name="unbounded_integers")
+    x = model.add_variable(name="x", lb=1, ub=float("inf"), is_integer=True)
+    y = model.add_variable(name="y", lb=4, ub=float("inf"), is_integer=True)
+    model.add_linear_constraint(x + y >= 2, name="c")
+    model.minimize(x + y)
+    return model.export_model()
+
+
+class MPSConverterTest(absltest.TestCase, compare_proto.MathOptProtoAssertions):
+
+    def test_convert_empty_mps_to_model_proto(self) -> None:
+        simple_mps = "NAME MIN_SIZE_MAX_FEATURES"
+        model_proto = mps_converter.mps_to_model_proto(simple_mps)
+        self.assertEqual(model_proto.name, "MIN_SIZE_MAX_FEATURES")
+
+    def test_convert_simple_mps_to_model(self) -> None:
+        model_proto = mps_converter.mps_to_model_proto(_MODEL_MPS)
+        expected_model_proto = simple_model_proto()
+
+        self.assert_protos_equiv(model_proto, expected_model_proto)
+
+    def test_convert_model_proto_to_mps(self) -> None:
+        model_proto = simple_model_proto()
+        mps = mps_converter.model_proto_to_mps(model_proto)
+        self.assertEqual(mps, _MODEL_MPS)
+
+
+if __name__ == "__main__":
+    absltest.main()
diff --git a/ortools/math_opt/labs/BUILD.bazel b/ortools/math_opt/labs/BUILD.bazel
index 7b436d27a7..d048d84df1 100644
--- a/ortools/math_opt/labs/BUILD.bazel
+++ b/ortools/math_opt/labs/BUILD.bazel
@@ -70,3 +70,19 @@ cc_library(
         "@com_google_absl//absl/types:span",
     ],
 )
+
+cc_library(
+    name = "dualizer",
+    srcs = [
+        "dualizer.cc",
+    ],
+    hdrs = ["dualizer.h"],
+    visibility = ["//visibility:public"],
+    deps = [
+        "//ortools/base:map_util",
+        "//ortools/math_opt/cpp:math_opt",
+        "@com_google_absl//absl/container:flat_hash_map",
+        "@com_google_absl//absl/types:span",
+    ],
+    alwayslink = 1,
+)
diff --git a/ortools/math_opt/labs/dualizer.cc b/ortools/math_opt/labs/dualizer.cc
new file mode 100644
index 0000000000..7e0f418ea0
--- /dev/null
+++ b/ortools/math_opt/labs/dualizer.cc
@@ -0,0 +1,211 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Let L be a matrix and b a vector so that a(w) = L * w + b. Then
+//
+// max_w{ a(w) * x : w in W} = max_w{ w' * L' * x : w in W} + b * x
+//
+// where ' is the transpose operation. Because of this we can focus on
+// max_w{ l(w) * x : w in W}.
+//
+// We need the dual to be an LP even when uncertainty_model contains ranged
+// constraints, so we use the LP reformulation of go/mathopt-dual from
+// go/mathopt-traditional-dual#lp-reformulation-split. Using
+// that reformulation, for any fixed x the dual of max_w{ w' * L' * x : w in W}
+// is
+//
+// min_{y, yp, yn, r, rp, rn}    obj(yp, yn, rp, rn)
+//
+//                       A' y + r == L' * x
+//                       sign constraints on y and r
+//                        yp + yn == y
+//                        rp + rn == r
+//                         yp, rp >= 0
+//                         yn, rn <= 0
+//
+// where
+//
+//   obj(yp, yn, rp, rn) = uc * yp + lc * yn + uv * rp + lv * rn
+//
+// with the convention 0 * infinity = 0 * -infinity = 0.
+//
+// In this dual form x is not multiplied with w so we can consider x a variable
+// instead of a fixed value.
+//
+// Then max_w{ a(w) * x : w in W} <= rhs is equivalent to
+//
+//          obj(yp, yn, rp, rn) + b * x <= rhs
+//                             A' y + r == L' * x
+//                             sign constraints on y and r
+//                              yp + yn == y
+//                              rp + rn == r
+//                               yp, rp >= 0
+//                               yn, rn <= 0
+//
+// Note that we can use the equalities yp + yn == y and rp + rn == r to
+// eliminate variables y and r to reduce the number of constraints and variables
+// in the reformulation.
+
+#include "ortools/math_opt/labs/dualizer.h"
+
+#include <limits>
+#include <utility>
+#include <vector>
+
+#include "absl/container/flat_hash_map.h"
+#include "absl/types/span.h"
+#include "ortools/base/map_util.h"
+#include "ortools/math_opt/cpp/math_opt.h"
+
+namespace operations_research {
+namespace math_opt {
+
+namespace {
+
+constexpr double kInf = std::numeric_limits<double>::infinity();
+
+class RobustConstraintDualizer {
+ public:
+  RobustConstraintDualizer(
+      const Model& uncertainty_model, Variable rhs,
+      absl::Span<const std::pair<LinearExpression, Variable>>
+          uncertain_coefficients,
+      Model& main_model);
+
+ private:
+  LinearExpression AddDualizedVariable(double lower_bound, double upper_bound);
+  Model& main_model_;
+
+  // Over the variables of main_model
+  LinearExpression objective_expression_;
+  // The keys are from the uncertain model, the values are from main_model
+  absl::flat_hash_map<LinearConstraint, LinearExpression> y_;
+  // The keys are from the uncertain model, the values are from main_model
+  absl::flat_hash_map<Variable, LinearExpression> r_;
+};
+
+// Let
+// (var, varp, varn, lower_bound, upper_bound) = (y_i, yp_i, yn_i, lc_i, uc_i)
+// or
+// (var, varp, varn, lower_bound, upper_bound) = (r_j, rp_j, rn_j, lv_j, uv_j).
+//
+// The constraints from go/mathopt-traditional-dual#lp-reformulation-split that
+// only involve var, varp and varn are (note that our dual has a max objective):
+//
+//           var >= 0    if    lower_bound = -infinity
+//           var <= 0    if    upper_bound = +infinity
+//   varp + varn == var
+//          varp >= 0
+//          varn <= 0
+//
+// and the corresponding term in obj(yp, yn, rp, rn) is
+//
+//   upper_bound * varp + lower_bound * varn
+//
+// The following function adds varp and varn, updates the expression for
+// obj(yp, yn, rp, rn) with the associated term and returns the expression for
+// var. The function uses the sign constraints on var, varp and varn and the
+// values of lower_bound and upper_bound to minimize the number of created
+// variables.
+LinearExpression RobustConstraintDualizer::AddDualizedVariable(
+    const double lower_bound, const double upper_bound) {
+  if ((lower_bound <= -kInf) && (upper_bound >= kInf)) {
+    return 0.0;
+  } else if (lower_bound <= -kInf) {
+    const Variable varp = main_model_.AddContinuousVariable(0.0, kInf);
+    objective_expression_ += upper_bound * varp;
+    return varp;
+  } else if (upper_bound >= kInf) {
+    const Variable varn = main_model_.AddContinuousVariable(-kInf, 0.0);
+    objective_expression_ += lower_bound * varn;
+    return varn;
+  } else if (lower_bound == upper_bound) {
+    const Variable var = main_model_.AddContinuousVariable(-kInf, kInf);
+    objective_expression_ += lower_bound * var;
+    return var;
+  } else {
+    const Variable varp = main_model_.AddContinuousVariable(0.0, kInf);
+    const Variable varn = main_model_.AddContinuousVariable(-kInf, 0.0);
+    objective_expression_ += upper_bound * varp + lower_bound * varn;
+    return varp + varn;
+  }
+}
+
+// L' * x
+absl::flat_hash_map<Variable, LinearExpression> TransposeUncertainCoefficients(
+    absl::Span<const std::pair<LinearExpression, Variable>>
+        uncertain_coefficients) {
+  absl::flat_hash_map<Variable, LinearExpression> result;
+  for (const auto& [expression, main_model_variable] : uncertain_coefficients) {
+    for (const auto [v, coefficient] : expression.terms()) {
+      result[v] += coefficient * main_model_variable;
+    }
+  }
+  return result;
+}
+
+RobustConstraintDualizer::RobustConstraintDualizer(
+    const Model& uncertainty_model, const Variable rhs,
+    absl::Span<const std::pair<LinearExpression, Variable>>
+        uncertain_coefficients,
+    Model& main_model)
+    : main_model_(main_model) {
+  const std::vector<Variable> uncertainty_variables =
+      uncertainty_model.SortedVariables();
+  const std::vector<LinearConstraint> uncertainty_constraints =
+      uncertainty_model.SortedLinearConstraints();
+  for (const LinearConstraint c : uncertainty_constraints) {
+    y_.insert({c, AddDualizedVariable(c.lower_bound(), c.upper_bound())});
+  }
+  for (const Variable v : uncertainty_variables) {
+    r_.insert({v, AddDualizedVariable(v.lower_bound(), v.upper_bound())});
+  }
+
+  // Add obj(yp, yn, rp, rn) + b * x <= rhs
+  {
+    LinearExpression offset_expression;
+    for (const auto& [expression, variable] : uncertain_coefficients) {
+      offset_expression += expression.offset() * variable;
+    }
+    main_model_.AddLinearConstraint(objective_expression_ + offset_expression <=
+                                    rhs);
+  }
+  {  // Add A' y + r = L' * x
+    const absl::flat_hash_map<Variable, LinearExpression>
+        equality_rhs_expressions =
+            TransposeUncertainCoefficients(uncertain_coefficients);
+    for (const Variable v : uncertainty_variables) {
+      LinearExpression equality_lhs_expression = r_.at(v);
+      for (const LinearConstraint c : uncertainty_model.ColumnNonzeros(v)) {
+        equality_lhs_expression += c.coefficient(v) * y_.at(c);
+      }
+      main_model_.AddLinearConstraint(
+          equality_lhs_expression ==
+          gtl::FindWithDefault(equality_rhs_expressions, v));
+    }
+  }
+}
+
+}  // namespace
+
+void AddRobustConstraint(
+    const Model& uncertainty_model, const Variable rhs,
+    const std::vector<std::pair<LinearExpression, Variable>>&
+        uncertain_coefficients,
+    Model& main_model) {
+  RobustConstraintDualizer dualizer(uncertainty_model, rhs,
+                                    uncertain_coefficients, main_model);
+}
+
+}  // namespace math_opt
+}  // namespace operations_research
diff --git a/ortools/math_opt/labs/dualizer.h b/ortools/math_opt/labs/dualizer.h
new file mode 100644
index 0000000000..9d2ed10562
--- /dev/null
+++ b/ortools/math_opt/labs/dualizer.h
@@ -0,0 +1,53 @@
+// Copyright 2010-2024 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_MATH_OPT_LABS_DUALIZER_H_
+#define OR_TOOLS_MATH_OPT_LABS_DUALIZER_H_
+
+#include <utility>
+#include <vector>
+
+#include "ortools/math_opt/cpp/math_opt.h"
+
+namespace operations_research {
+namespace math_opt {
+
+// Uses LP duality to construct an extended formulation of
+//
+//   max_w{ a(w) * x : w in W} <= rhs
+//
+// where W is described by uncertainty_model (the variables of uncertainty_model
+// are w). All the variables and constraints of the extended formulation are
+// added to main_model.
+//
+// Requirements:
+//   * x must be variables of main_model
+//   * rhs must be a variable of main_model
+//   * uncertainty_model must be an LP
+//   * uncertain coefficient a(w)_i for x_i should be a LinearExpression of w.
+//
+// Input-only arguments:
+//   * uncertainty_model
+//   * rhs
+//   * uncertain_coefficients: pairs [a(w)_i, x_i] for all i
+// Input-output argument:
+//   * main_model
+void AddRobustConstraint(
+    const Model& uncertainty_model, Variable rhs,
+    const std::vector<std::pair<LinearExpression, Variable>>&
+        uncertain_coefficients,
+    Model& main_model);
+
+}  // namespace math_opt
+}  // namespace operations_research
+#endif  // OR_TOOLS_MATH_OPT_LABS_DUALIZER_H_
diff --git a/ortools/math_opt/model.proto b/ortools/math_opt/model.proto
index 1b71df8658..b53b9467e2 100644
--- a/ortools/math_opt/model.proto
+++ b/ortools/math_opt/model.proto
@@ -269,7 +269,7 @@ message ModelProto {
   //  * linear_constraint_matrix.row_ids are elements of linear_constraints.ids.
   //  * linear_constraint_matrix.column_ids are elements of variables.ids.
   //  * Matrix entries not specified are zero.
-  //  * linear_constraint_matrix.values must all be finite.
+  //  * linear_constraint_matrix.coefficients must all be finite.
   SparseDoubleMatrixProto linear_constraint_matrix = 5;
 
   // Mapped constraints (i.e., stored in "constraint ID"-to-"constraint data"
diff --git a/ortools/math_opt/model_parameters.proto b/ortools/math_opt/model_parameters.proto
index 7d9fb34700..5628674746 100644
--- a/ortools/math_opt/model_parameters.proto
+++ b/ortools/math_opt/model_parameters.proto
@@ -111,6 +111,14 @@ message ModelSolveParametersProto {
   //  * filtered_ids are elements of LinearConstraints.ids.
   SparseVectorFilterProto dual_values_filter = 2;
 
+  // Filter that is applied to all returned sparse containers keyed by quadratic
+  // constraints in DualSolutionProto and DualRay
+  // (DualSolutionProto.quadratic_dual_values, DualRay.quadratic_dual_values).
+  //
+  // Requirements:
+  //  * filtered_ids are keys of ModelProto.quadratic_constraints.
+  SparseVectorFilterProto quadratic_dual_values_filter = 10;
+
   // Filter that is applied to all returned sparse containers keyed by variables
   // in DualSolutionProto and DualRay (DualSolutionProto.reduced_costs,
   // DualRay.reduced_costs).
diff --git a/ortools/math_opt/python/BUILD.bazel b/ortools/math_opt/python/BUILD.bazel
index 2219ec5456..170edd664b 100644
--- a/ortools/math_opt/python/BUILD.bazel
+++ b/ortools/math_opt/python/BUILD.bazel
@@ -25,6 +25,7 @@ py_library(
     deps = [
         ":callback",
         ":compute_infeasible_subsystem_result",
+        ":errors",
         ":expressions",
         ":hash_model_storage",
         ":message_callback",
@@ -158,12 +159,14 @@ py_library(
     deps = [
         ":callback",
         ":compute_infeasible_subsystem_result",
+        ":errors",
         ":message_callback",
         ":model",
         ":model_parameters",
         ":parameters",
         ":result",
         "//ortools/math_opt:parameters_py_pb2",
+        "//ortools/math_opt:rpc_py_pb2",
         "//ortools/math_opt/core/python:solver",
     ],
 )
@@ -201,3 +204,9 @@ py_library(
     srcs = ["solver_resources.py"],
     deps = ["//ortools/math_opt:rpc_py_pb2"],
 )
+
+py_library(
+    name = "errors",
+    srcs = ["errors.py"],
+    deps = ["//ortools/math_opt:rpc_py_pb2"],
+)
diff --git a/ortools/math_opt/python/errors.py b/ortools/math_opt/python/errors.py
new file mode 100644
index 0000000000..ae80c14f52
--- /dev/null
+++ b/ortools/math_opt/python/errors.py
@@ -0,0 +1,104 @@
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Translate C++'s absl::Status errors to Python standard errors.
+
+Here we try to use the standard Python errors we would use if the C++ code was
+instead implemented in Python. This will give Python users a more familiar API.
+"""
+
+import enum
+from typing import Optional, Type
+from ortools.math_opt import rpc_pb2
+
+
+class _StatusCode(enum.Enum):
+    """The C++ absl::Status::code() values."""
+
+    OK = 0
+    CANCELLED = 1
+    UNKNOWN = 2
+    INVALID_ARGUMENT = 3
+    DEADLINE_EXCEEDED = 4
+    NOT_FOUND = 5
+    ALREADY_EXISTS = 6
+    PERMISSION_DENIED = 7
+    UNAUTHENTICATED = 16
+    RESOURCE_EXHAUSTED = 8
+    FAILED_PRECONDITION = 9
+    ABORTED = 10
+    OUT_OF_RANGE = 11
+    UNIMPLEMENTED = 12
+    INTERNAL = 13
+    UNAVAILABLE = 14
+    DATA_LOSS = 15
+
+
+class InternalMathOptError(RuntimeError):
+    """Some MathOpt internal error.
+
+    This error is usually raised because of a bug in MathOpt or one of the solver
+    library it wraps.
+    """
+
+
+def status_proto_to_exception(
+    status_proto: rpc_pb2.StatusProto,
+) -> Optional[Exception]:
+    """Returns the Python exception that best match the input absl::Status.
+
+    There are some Status that we expect the MathOpt code to return, for those the
+    matching exceptions are:
+    - InvalidArgument: ValueError
+    - FailedPrecondition: AssertionError
+    - Unimplemented: NotImplementedError
+    - Internal: InternalMathOptError
+
+    Other Status's are not used by MathOpt, if they are seen a
+    InternalMathOptError is raised (as if the Status was Internal) and the error
+    message contains the unexpected code.
+
+    Args:
+      status_proto: The input proto to convert to an exception.
+
+    Returns:
+      The corresponding exception. None if the input status is OK.
+    """
+    try:
+        code = _StatusCode(status_proto.code)
+    except ValueError:
+        return InternalMathOptError(
+            f"unknown C++ error (code = {status_proto.code}):"
+            f" {status_proto.message}"
+        )
+
+    if code == _StatusCode.OK:
+        return None
+
+    # For expected errors we compute the corresponding class.
+    error_type: Optional[Type[Exception]] = None
+    if code == _StatusCode.INVALID_ARGUMENT:
+        error_type = ValueError
+    if code == _StatusCode.FAILED_PRECONDITION:
+        error_type = AssertionError
+    if code == _StatusCode.UNIMPLEMENTED:
+        error_type = NotImplementedError
+    if code == _StatusCode.INTERNAL:
+        error_type = InternalMathOptError
+
+    if error_type is not None:
+        return error_type(f"{status_proto.message} (was C++ {code.name})")
+
+    return InternalMathOptError(
+        f"unexpected C++ error {code.name}: {status_proto.message}"
+    )
diff --git a/ortools/math_opt/python/errors_test.py b/ortools/math_opt/python/errors_test.py
new file mode 100644
index 0000000000..e75af6a420
--- /dev/null
+++ b/ortools/math_opt/python/errors_test.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Tests of the `errors` package."""
+
+from absl.testing import absltest
+from ortools.math_opt import rpc_pb2
+from ortools.math_opt.python import errors
+
+
+class StatusProtoToExceptionTest(absltest.TestCase):
+
+    def test_ok(self) -> None:
+        self.assertIsNone(
+            errors.status_proto_to_exception(
+                rpc_pb2.StatusProto(code=errors._StatusCode.OK.value)
+            )
+        )
+
+    def test_invalid_argument(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(
+                code=errors._StatusCode.INVALID_ARGUMENT.value, message="something"
+            )
+        )
+        self.assertIsInstance(error, ValueError)
+        self.assertEqual(str(error), "something (was C++ INVALID_ARGUMENT)")
+
+    def test_failed_precondition(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(
+                code=errors._StatusCode.FAILED_PRECONDITION.value,
+                message="something",
+            )
+        )
+        self.assertIsInstance(error, AssertionError)
+        self.assertEqual(str(error), "something (was C++ FAILED_PRECONDITION)")
+
+    def test_unimplemented(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(
+                code=errors._StatusCode.UNIMPLEMENTED.value, message="something"
+            )
+        )
+        self.assertIsInstance(error, NotImplementedError)
+        self.assertEqual(str(error), "something (was C++ UNIMPLEMENTED)")
+
+    def test_internal(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(
+                code=errors._StatusCode.INTERNAL.value, message="something"
+            )
+        )
+        self.assertIsInstance(error, errors.InternalMathOptError)
+        self.assertEqual(str(error), "something (was C++ INTERNAL)")
+
+    def test_unexpected_code(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(
+                code=errors._StatusCode.DEADLINE_EXCEEDED.value, message="something"
+            )
+        )
+        self.assertIsInstance(error, errors.InternalMathOptError)
+        self.assertEqual(
+            str(error), "unexpected C++ error DEADLINE_EXCEEDED: something"
+        )
+
+    def test_unknown_code(self) -> None:
+        error = errors.status_proto_to_exception(
+            rpc_pb2.StatusProto(code=-5, message="something")
+        )
+        self.assertIsInstance(error, errors.InternalMathOptError)
+        self.assertEqual(str(error), "unknown C++ error (code = -5): something")
+
+
+if __name__ == "__main__":
+    absltest.main()
diff --git a/ortools/math_opt/python/ipc/remote_http_solve.py b/ortools/math_opt/python/ipc/remote_http_solve.py
index 8e210f5ad7..f45abd16de 100644
--- a/ortools/math_opt/python/ipc/remote_http_solve.py
+++ b/ortools/math_opt/python/ipc/remote_http_solve.py
@@ -96,7 +96,7 @@ def create_optimization_service_session(
 
     Returns:
       requests.Session a session with the necessary headers to call the
-      optimization serive.
+      optimization service.
     """
     session = requests.Session()
     server_timeout = deadline_sec * (1 - _RELATIVE_TIME_BUFFER)
diff --git a/ortools/math_opt/python/mathopt.py b/ortools/math_opt/python/mathopt.py
index 4a09d9f1a4..8191122344 100644
--- a/ortools/math_opt/python/mathopt.py
+++ b/ortools/math_opt/python/mathopt.py
@@ -60,6 +60,8 @@ from ortools.math_opt.python.compute_infeasible_subsystem_result import (
 from ortools.math_opt.python.compute_infeasible_subsystem_result import (
     parse_model_subset_bounds,
 )
+from ortools.math_opt.python.errors import InternalMathOptError
+from ortools.math_opt.python.errors import status_proto_to_exception
 from ortools.math_opt.python.expressions import evaluate_expression
 from ortools.math_opt.python.expressions import fast_sum
 from ortools.math_opt.python.hash_model_storage import HashModelStorage
diff --git a/ortools/math_opt/python/solve.py b/ortools/math_opt/python/solve.py
index 49da39f0c8..500453fe12 100644
--- a/ortools/math_opt/python/solve.py
+++ b/ortools/math_opt/python/solve.py
@@ -16,9 +16,11 @@ import types
 from typing import Callable, Optional
 
 from ortools.math_opt import parameters_pb2
+from ortools.math_opt import rpc_pb2
 from ortools.math_opt.core.python import solver
 from ortools.math_opt.python import callback
 from ortools.math_opt.python import compute_infeasible_subsystem_result
+from ortools.math_opt.python import errors
 from ortools.math_opt.python import message_callback
 from ortools.math_opt.python import model
 from ortools.math_opt.python import model_parameters
@@ -86,7 +88,7 @@ def solve(
             None,
         )
     except StatusNotOk as e:
-        raise RuntimeError(str(e)) from None
+        raise _status_not_ok_to_exception(e) from None
     return result.parse_solve_result(proto_result, opt_model)
 
 
@@ -126,7 +128,7 @@ def compute_infeasible_subsystem(
             None,
         )
     except StatusNotOk as e:
-        raise RuntimeError(str(e)) from None
+        raise _status_not_ok_to_exception(e) from None
     return (
         compute_infeasible_subsystem_result.parse_compute_infeasible_subsystem_result(
             proto_result, opt_model
@@ -172,7 +174,7 @@ class IncrementalSolver:
                 parameters_pb2.SolverInitializerProto(),
             )
         except StatusNotOk as e:
-            raise RuntimeError(str(e)) from None
+            raise _status_not_ok_to_exception(e) from None
         self._closed = False
 
     def solve(
@@ -212,7 +214,7 @@ class IncrementalSolver:
                         parameters_pb2.SolverInitializerProto(),
                     )
             except StatusNotOk as e:
-                raise RuntimeError(str(e)) from None
+                raise _status_not_ok_to_exception(e) from None
             self._update_tracker.advance_checkpoint()
         params = params or parameters.SolveParameters()
         model_params = model_params or model_parameters.ModelSolveParameters()
@@ -232,7 +234,7 @@ class IncrementalSolver:
                 None,
             )
         except StatusNotOk as e:
-            raise RuntimeError(str(e)) from None
+            raise _status_not_ok_to_exception(e) from None
         return result.parse_solve_result(result_proto, self._model)
 
     def close(self) -> None:
@@ -268,3 +270,20 @@ class IncrementalSolver:
     ) -> None:
         """Closes the solver."""
         self.close()
+
+
+def _status_not_ok_to_exception(err: StatusNotOk) -> Exception:
+    """Converts a StatusNotOk to the best matching Python exception.
+
+    Args:
+      err: The input errors.
+
+    Returns:
+      The corresponding exception.
+    """
+    ret = errors.status_proto_to_exception(
+        rpc_pb2.StatusProto(code=err.canonical_code, message=err.message)
+    )
+    # We never expect StatusNotOk to be OK.
+    assert ret is not None, err
+    return ret
diff --git a/ortools/math_opt/python/solve_test.py b/ortools/math_opt/python/solve_test.py
index 6fcefa2759..3d432d0f85 100644
--- a/ortools/math_opt/python/solve_test.py
+++ b/ortools/math_opt/python/solve_test.py
@@ -68,9 +68,7 @@ class SolveTest(absltest.TestCase):
     def test_solve_error(self) -> None:
         mod = model.Model(name="test_model")
         mod.add_variable(lb=1.0, ub=-1.0, name="x1")
-        with self.assertRaisesRegex(
-            RuntimeError, "variables.*lower_bound > upper_bound"
-        ):
+        with self.assertRaisesRegex(ValueError, "variables.*lower_bound > upper_bound"):
             solve.solve(mod, parameters.SolverType.GLOP)
 
     def test_lp_solve(self) -> None:
@@ -213,16 +211,14 @@ class SolveTest(absltest.TestCase):
         mod = model.Model(name="test_model")
         mod.add_variable(lb=1.0, ub=1.0, name="x1")
         mod.add_variable(lb=1.0, ub=1.0, name="x1")
-        with self.assertRaisesRegex(RuntimeError, "duplicate name*"):
+        with self.assertRaisesRegex(ValueError, "duplicate name*"):
             solve.IncrementalSolver(mod, parameters.SolverType.GLOP)
 
     def test_incremental_solve_error(self) -> None:
         mod = model.Model(name="test_model")
         mod.add_variable(lb=1.0, ub=-1.0, name="x1")
         solver = solve.IncrementalSolver(mod, parameters.SolverType.GLOP)
-        with self.assertRaisesRegex(
-            RuntimeError, "variables.*lower_bound > upper_bound"
-        ):
+        with self.assertRaisesRegex(ValueError, "variables.*lower_bound > upper_bound"):
             solver.solve()
 
     def test_incremental_solve_error_on_reject(self) -> None:
@@ -250,7 +246,7 @@ class SolveTest(absltest.TestCase):
         )
 
         opt_model.add_binary_variable(name="x")
-        with self.assertRaisesRegex(RuntimeError, "duplicate name*"):
+        with self.assertRaisesRegex(ValueError, "duplicate name*"):
             solver.solve(
                 msg_cb=message_callback.printer_message_callback(prefix="[solve 2] ")
             )
diff --git a/ortools/math_opt/rpc.proto b/ortools/math_opt/rpc.proto
index 7ac5e8505d..fcbc425f4e 100644
--- a/ortools/math_opt/rpc.proto
+++ b/ortools/math_opt/rpc.proto
@@ -98,8 +98,8 @@ message SolveResponse {
   // Either `result` or `status` must be set. This is equivalent to C++
   // StatusOr<SolveResult>.
   oneof status_or {
-  // Description of the output of solving the model in the request.
-  SolveResultProto result = 1;
+    // Description of the output of solving the model in the request.
+    SolveResultProto result = 1;
 
     // The absl::Status returned by the solver. It should never be OK when set.
     StatusProto status = 3;
diff --git a/ortools/math_opt/samples/cpp/cutting_stock.cc b/ortools/math_opt/samples/cpp/cutting_stock.cc
index 41bf833c24..2e39111d5c 100644
--- a/ortools/math_opt/samples/cpp/cutting_stock.cc
+++ b/ortools/math_opt/samples/cpp/cutting_stock.cc
@@ -189,7 +189,7 @@ absl::StatusOr<CuttingStockSolution> SolveCuttingStock(
     add_config(Configuration{.pieces = {i}, .quantity = {1}});
   }
 
-  ASSIGN_OR_RETURN(auto solver, math_opt::IncrementalSolver::New(
+  ASSIGN_OR_RETURN(auto solver, math_opt::NewIncrementalSolver(
                                     &model, math_opt::SolverType::kGlop));
   int pricing_round = 0;
   while (true) {
diff --git a/ortools/math_opt/samples/cpp/facility_lp_benders.cc b/ortools/math_opt/samples/cpp/facility_lp_benders.cc
index 04316619b0..8cb0a89997 100644
--- a/ortools/math_opt/samples/cpp/facility_lp_benders.cc
+++ b/ortools/math_opt/samples/cpp/facility_lp_benders.cc
@@ -376,7 +376,7 @@ absl::StatusOr<std::unique_ptr<SecondStageSolver>> SecondStageSolver::New(
       absl::WrapUnique<SecondStageSolver>(
           new SecondStageSolver(std::move(instance), parameters));
   ASSIGN_OR_RETURN(std::unique_ptr<math_opt::IncrementalSolver> solver,
-                   math_opt::IncrementalSolver::New(
+                   math_opt::NewIncrementalSolver(
                        &second_stage_solver->second_stage_model_, solver_type));
   second_stage_solver->solver_ = std::move(solver);
   return std::move(second_stage_solver);
@@ -593,7 +593,7 @@ absl::Status Benders(const FacilityLocationInstance& instance,
   FirstStageProblem first_stage(instance.network, instance.facility_cost);
   ASSIGN_OR_RETURN(
       const std::unique_ptr<math_opt::IncrementalSolver> first_stage_solver,
-      math_opt::IncrementalSolver::New(&first_stage.model, solver_type));
+      math_opt::NewIncrementalSolver(&first_stage.model, solver_type));
   // Setup second stage solver.
   ASSIGN_OR_RETURN(std::unique_ptr<SecondStageSolver> second_stage_solver,
                    SecondStageSolver::New(instance, solver_type));
diff --git a/ortools/math_opt/samples/python/BUILD.bazel b/ortools/math_opt/samples/python/BUILD.bazel
index 80ddaa9ad2..838353a7e2 100644
--- a/ortools/math_opt/samples/python/BUILD.bazel
+++ b/ortools/math_opt/samples/python/BUILD.bazel
@@ -113,3 +113,13 @@ py_binary(
         "//ortools/math_opt/python:mathopt",
     ],
 )
+
+py_binary(
+    name = "writing_seminar",
+    srcs = ["writing_seminar.py"],
+    deps = [
+        requirement("absl-py"),
+        "//ortools/math_opt/python:mathopt",
+        "//ortools/math_opt/solvers:highs_solver",
+    ],
+)
diff --git a/ortools/math_opt/samples/python/writing_seminar.py b/ortools/math_opt/samples/python/writing_seminar.py
new file mode 100644
index 0000000000..a62ce110fb
--- /dev/null
+++ b/ortools/math_opt/samples/python/writing_seminar.py
@@ -0,0 +1,266 @@
+#!/usr/bin/env python3
+# Copyright 2010-2024 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Solves the problem of assigning students to classes to maximize welfare.
+
+Each student provides a ranking of their top five preferred writing seminars to
+be in, and can be assigned to at most one seminar. Each writing seminar has a
+capacity of how many students can attend. There is an increasing penalty for
+assigning students to their less preferred choice, e.g.
+
+assignment    | penalty
+--------------|--------
+first choice  | 0
+second choice | 1
+third choice  | 5
+fourth choice | 20
+fifth choice  | 100
+unassigned    | 1000
+
+(Unassigned can mean either they are in a class they did not choose, or they are
+not in any class). The goal is to find a feasible assignment that minimizes the
+sum of the penalties over all students.
+
+This method was used as of 2005 at Cornell to assign first year students to
+freshman writing seminars.
+
+We model the problem below with MIP, but it is in fact just a min-cost flow
+problem. Solving the LP relaxation of this MIP with a simplex based solver
+gives integer solutions.
+
+Warning: this method is not strategy proof, students can get a better assignment
+by not submitting their true preferences over the seminars.
+"""
+
+import dataclasses
+import datetime
+import random
+from typing import Dict, Sequence, Tuple
+
+from absl import app
+from absl import flags
+
+from ortools.math_opt.python import mathopt
+
+
+_SOLVER_TYPE = flags.DEFINE_enum_class(
+    "solver_type",
+    mathopt.SolverType.CP_SAT,
+    mathopt.SolverType,
+    "The solver to use, must be an LP solver if lp_relax=True (e.g. GLOP) and"
+    " must be an IP solver (e.g. SCIP, CP-SAT) otherwise.",
+)
+
+_NUM_CLASSES = flags.DEFINE_integer(
+    "num_classes", 50, "How many random classes to create."
+)
+
+_LP_RELAX = flags.DEFINE_boolean(
+    "lp_relax", False, "Solve the LP relaxation of the problem."
+)
+
+_TEST_DATA = flags.DEFINE_boolean("test_data", False, "Use the small test instance.")
+
+
+@dataclasses.dataclass(frozen=True)
+class Student:
+    preferred_classes: Tuple[int, ...]
+    name: str = ""
+
+
+@dataclasses.dataclass(frozen=True)
+class Seminar:
+    capacity: int
+    name: str = ""
+
+
+@dataclasses.dataclass(frozen=True)
+class WritingSeminarAssignmentProblem:
+    seminars: Tuple[Seminar, ...]
+    students: Tuple[Student, ...]
+    rank_penalty: Tuple[float, ...]
+    unmatched_penalty: float
+
+
+def _test_problem() -> WritingSeminarAssignmentProblem:
+    """A small deterministic instance for testing only."""
+    return WritingSeminarAssignmentProblem(
+        seminars=(Seminar(capacity=1, name="c1"), Seminar(capacity=1, name="c2")),
+        students=(
+            Student(preferred_classes=(1, 0), name="s1"),
+            Student(preferred_classes=(0, 1), name="s2"),
+        ),
+        rank_penalty=(0, 10),
+        unmatched_penalty=100,
+    )
+
+
+def _random_writing_seminar_assignment_problem(
+    *,
+    num_classes: int,
+    class_capacity: int,
+    num_students: int,
+    selections_per_student: int,
+    unmatched_penalty: float,
+    rank_penalty: Tuple[float, ...],
+) -> WritingSeminarAssignmentProblem:
+    """Generates a random instance of the WritingSeminarAssignmentProblem."""
+    if len(rank_penalty) != selections_per_student:
+        raise ValueError(
+            f"len(rank_penalty): {len(rank_penalty)} must equal"
+            f" selections_per_student: {selections_per_student}"
+        )
+    seminars = tuple(
+        Seminar(capacity=class_capacity, name=f"c_{i}") for i in range(num_classes)
+    )
+    students = []
+    all_class_ids = list(range(num_classes))
+    for s in range(num_students):
+        preferred = tuple(random.sample(all_class_ids, selections_per_student))
+        students.append(Student(preferred_classes=preferred, name=f"s_{s}"))
+    return WritingSeminarAssignmentProblem(
+        seminars=seminars,
+        students=tuple(students),
+        rank_penalty=rank_penalty,
+        unmatched_penalty=unmatched_penalty,
+    )
+
+
+def _assign_students(
+    problem: WritingSeminarAssignmentProblem,
+    solver_type: mathopt.SolverType,
+    time_limit: datetime.timedelta,
+) -> Dict[Student, Seminar]:
+    """Optimally assigns students to classes by solving an IP."""
+    # Problem data:
+    #   * i in S: the students.
+    #   * j in C: the classes (writing seminars).
+    #   * c_j: the capacity of class j.
+    #   * K: how many classes each student ranks.
+    #   * R_i for i in S, an ordered list K classes ranked for student i.
+    #   * p_k for k = 1,...,K the penalty for giving a student their kth choice.
+    #   * P: the penalty for not assigning a student.
+    #
+    # Decision variables:
+    #   x_i_j: student i takes seminar j
+    #   y_i: student i is not assigned
+    #
+    # Problem formulation:
+    #   min     sum_i P y_i +  sum_{k, j in enumerate(R_i)} p_k x_i_j
+    #   s.t.    y_i + sum_{j in R_i} x_i_j = 1    for all i in S
+    #           sum_{i : j in R_i} x_i_j <= c_j   for all j in C
+    #           x_i_j in {0, 1}                   for all i in S, for all j in R_i
+    #           y_i in {0, 1}                     for all i in S
+    model = mathopt.Model()
+    use_int_vars = not _LP_RELAX.value
+    # The y_i variables.
+    unassigned = [
+        model.add_variable(lb=0, ub=1, is_integer=use_int_vars, name=f"y_{i}")
+        for i, _ in enumerate(problem.students)
+    ]
+
+    # The x_i_j variables.
+    assignment_vars = [{} for _ in range(len(problem.students))]
+    for i, student in enumerate(problem.students):
+        for rank, j in enumerate(student.preferred_classes):
+            assignment_vars[i][j] = model.add_variable(
+                lb=0, ub=1, is_integer=use_int_vars, name=f"x_{i}_{j}"
+            )
+    # Transpose the variables in x. The first index of students_in_seminar
+    # is the class (j). The value is an unordered list of the variables for
+    # assigning students into this class.
+    students_in_seminar = [[] for _ in range(len(problem.seminars))]
+    for seminar_to_x in assignment_vars:
+        for j, x in seminar_to_x.items():
+            students_in_seminar[j].append(x)
+
+    # Create the objective
+    penalties = mathopt.fast_sum(unassigned) * problem.unmatched_penalty
+    for i, student in enumerate(problem.students):
+        penalties += mathopt.fast_sum(
+            problem.rank_penalty[rank] * assignment_vars[i][j]
+            for rank, j in enumerate(student.preferred_classes)
+        )
+    model.minimize(penalties)
+
+    # Each student is in at most one class
+    for i, student in enumerate(problem.students):
+        model.add_linear_constraint(
+            unassigned[i] + mathopt.fast_sum(assignment_vars[i].values()) == 1.0
+        )
+
+    # Each class does not exceed its capacity
+    for j, seminar in enumerate(problem.seminars):
+        model.add_linear_constraint(
+            mathopt.fast_sum(students_in_seminar[j]) <= seminar.capacity
+        )
+
+    solve_result = mathopt.solve(
+        model,
+        solver_type,
+        params=mathopt.SolveParameters(enable_output=True, time_limit=time_limit),
+    )
+    if solve_result.termination.reason not in {
+        mathopt.TerminationReason.OPTIMAL,
+        mathopt.TerminationReason.FEASIBLE,
+    }:
+        raise RuntimeError(
+            f"failed to find a feasible solution: {solve_result.termination}"
+        )
+
+    assignment = {}
+    for i, student in enumerate(problem.students):
+        for sem, x in assignment_vars[i].items():
+            x_val = solve_result.variable_values(x)
+            int_err = min(abs(x_val), abs(1 - x_val))
+            if int_err > 1e-3:
+                raise RuntimeError(
+                    "all variables should be within 1e-3 of either 0 or 1, but found"
+                    f" value: {x_val}"
+                )
+            if solve_result.variable_values(x) > 0.5:
+                assignment[student] = problem.seminars[sem]
+    return assignment
+
+
+def main(args: Sequence[str]) -> None:
+    del args  # Unused.
+    if _TEST_DATA.value:
+        problem = _test_problem()
+    else:
+        num_classes = _NUM_CLASSES.value
+        class_capacity = 15
+        num_students = num_classes * 12
+        selections_per_student = 5
+        unmatched_penalty = 1000
+        rank_penalty = (0, 1, 5, 20, 100)
+        problem = _random_writing_seminar_assignment_problem(
+            num_classes=num_classes,
+            class_capacity=class_capacity,
+            num_students=num_students,
+            selections_per_student=selections_per_student,
+            unmatched_penalty=unmatched_penalty,
+            rank_penalty=rank_penalty,
+        )
+    assignment = _assign_students(
+        problem, _SOLVER_TYPE.value, datetime.timedelta(minutes=1)
+    )
+    for student, seminar in assignment.items():
+        print(f"{student.name}: {seminar.name}")
+    num_unassigned = len(problem.students) - len(assignment)
+    print(f"Unassigned students: {num_unassigned}")
+
+
+if __name__ == "__main__":
+    app.run(main)
diff --git a/ortools/math_opt/solution.proto b/ortools/math_opt/solution.proto
index 31b9682f78..c34a352832 100644
--- a/ortools/math_opt/solution.proto
+++ b/ortools/math_opt/solution.proto
@@ -104,6 +104,14 @@ message DualSolutionProto {
   //  * dual_values.values must all be finite.
   SparseDoubleVectorProto dual_values = 1;
 
+  // Requirements:
+  //  * quadratic_dual_values.ids are keys of ModelProto.quadratic_constraints.
+  //  * quadratic_dual_values.values must all be finite.
+  // Note: Some solvers only return quadratic constraint duals when a
+  // solver-specific parameter is set
+  // (see go/mathopt-qcqp-dual#solver-specific).
+  SparseDoubleVectorProto quadratic_dual_values = 5;
+
   // Requirements:
   //  * reduced_costs.ids are elements of VariablesProto.ids.
   //  * reduced_costs.values must all be finite.
diff --git a/ortools/math_opt/solver_tests/callback_tests.cc b/ortools/math_opt/solver_tests/callback_tests.cc
index e9cbf801f2..a688549d49 100644
--- a/ortools/math_opt/solver_tests/callback_tests.cc
+++ b/ortools/math_opt/solver_tests/callback_tests.cc
@@ -335,7 +335,7 @@ TEST_P(CallbackTest, EventSimplex) {
   // solve, we know the starting basis. It would be simpler to set the starting
   // basis, once this is supported.
   ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
-                       IncrementalSolver::New(&model, GetParam().solver_type));
+                       NewIncrementalSolver(&model, GetParam().solver_type));
   {
     ASSERT_OK_AND_ASSIGN(const SolveResult result, solver->Solve(args));
     ASSERT_THAT(result, IsOptimal(6.0));
diff --git a/ortools/math_opt/solver_tests/generic_tests.cc b/ortools/math_opt/solver_tests/generic_tests.cc
index deccf73769..675eb8c4af 100644
--- a/ortools/math_opt/solver_tests/generic_tests.cc
+++ b/ortools/math_opt/solver_tests/generic_tests.cc
@@ -20,6 +20,7 @@
 #include <ostream>
 #include <sstream>
 #include <string>
+#include <tuple>
 #include <utility>
 #include <vector>
 
@@ -442,9 +443,8 @@ TEST_P(GenericTest, InvertedVariableBounds) {
     model.Maximize(3.0 * x);
 
     // The instantiation should not fail, even if the bounds are reversed.
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     // Solving should fail because of the inverted bounds.
     EXPECT_THAT(solver->Solve(solve_args),
@@ -482,9 +482,8 @@ TEST_P(GenericTest, InvertedVariableBounds) {
 
     model.Maximize(3.0 * x);
 
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     // As of 2022-11-17 the glp_interior() algorithm returns GLP_EFAIL when the
     // model is "empty" (no rows or columns). The issue is that the emptiness is
@@ -531,15 +530,14 @@ TEST_P(GenericTest, InvertedVariableBounds) {
 
     model.Maximize(3.0 * x);
 
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     EXPECT_THAT(solver->SolveWithoutUpdate(solve_args),
                 IsOkAndHolds(IsOptimal(3.0 * 4.0)));
 
     // Test the update using a new variable with inverted bounds (in case the
-    // update code path is not identical to the IncrementalSolver::New() one).
+    // update code path is not identical to the NewIncrementalSolver() one).
     const Variable y = model.AddVariable(/*lower_bound=*/lb, /*upper_bound=*/ub,
                                          GetParam().integer_variables, "y");
     model.Maximize(3.0 * x + y);
@@ -573,9 +571,8 @@ TEST_P(GenericTest, InvertedLinearConstraintBounds) {
     model.Maximize(3.0 * x);
 
     // The instantiation should not fail, even if the bounds are reversed.
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     // Solving should fail because of the inverted bounds.
     EXPECT_THAT(solver->Solve(solve_args),
@@ -602,9 +599,8 @@ TEST_P(GenericTest, InvertedLinearConstraintBounds) {
 
     model.Maximize(3.0 * x);
 
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     EXPECT_THAT(solver->SolveWithoutUpdate(solve_args),
                 IsOkAndHolds(IsOptimal(3.0 * 4.0)));
@@ -638,15 +634,14 @@ TEST_P(GenericTest, InvertedLinearConstraintBounds) {
 
     model.Maximize(3.0 * x);
 
-    ASSERT_OK_AND_ASSIGN(
-        const std::unique_ptr<IncrementalSolver> solver,
-        IncrementalSolver::New(&model, GetParam().solver_type));
+    ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                         NewIncrementalSolver(&model, GetParam().solver_type));
 
     EXPECT_THAT(solver->SolveWithoutUpdate(solve_args),
                 IsOkAndHolds(IsOptimal(3.0 * 4.0)));
 
     // Test the update with a new constraint with inverted bounds (in case the
-    // update code path is not identical to the IncrementalSolver::New() one).
+    // update code path is not identical to the NewIncrementalSolver() one).
     const LinearConstraint v = model.AddLinearConstraint(5.0 <= x <= 3.0, "v");
 
     ASSERT_OK(solver->Update());
diff --git a/ortools/math_opt/solver_tests/ip_model_solve_parameters_tests.cc b/ortools/math_opt/solver_tests/ip_model_solve_parameters_tests.cc
index df7b5a3628..ce434aecf9 100644
--- a/ortools/math_opt/solver_tests/ip_model_solve_parameters_tests.cc
+++ b/ortools/math_opt/solver_tests/ip_model_solve_parameters_tests.cc
@@ -310,7 +310,7 @@ TEST_P(BranchPrioritiesTest, PrioritiesClearedAfterIncrementalSolve) {
   // tight node limit. We expect the solver to load the priorities, but not to
   // make any progress towards the optimal solution.
   ASSERT_OK_AND_ASSIGN(const auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
   {
     SolveParameters params = SolveParams();
     params.node_limit = 0;
@@ -407,7 +407,7 @@ TEST_P(LazyConstraintsTest, AnnotationsAreClearedAfterSolve) {
   const LinearConstraint d = model.AddLinearConstraint(y >= -x);
   model.Minimize(y);
   ASSERT_OK_AND_ASSIGN(const auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
 
   SolveArguments args = {
       .parameters = NerfedSolveParams(),
diff --git a/ortools/math_opt/solver_tests/ip_parameter_tests.cc b/ortools/math_opt/solver_tests/ip_parameter_tests.cc
index bcfc4f3eea..883b4584ca 100644
--- a/ortools/math_opt/solver_tests/ip_parameter_tests.cc
+++ b/ortools/math_opt/solver_tests/ip_parameter_tests.cc
@@ -1099,9 +1099,6 @@ TEST_P(LargeInstanceIpParameterTest, IterationLimit) {
 }
 
 TEST_P(LargeInstanceIpParameterTest, NodeLimit) {
-  if (GetParam().solver_type == SolverType::kHighs) {
-    GTEST_SKIP() << "Ignoring this test as Highs 1.7+ returns unimplemented";
-  }
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<Model> model, Load23588());
   SolveParameters params = GetParam().base_parameters;
   params.node_limit = 1;
@@ -1222,9 +1219,6 @@ TEST_P(LargeInstanceIpParameterTest, BestBoundLimit) {
 }
 
 TEST_P(LargeInstanceIpParameterTest, SolutionLimit) {
-  if (GetParam().solver_type == SolverType::kHighs) {
-    GTEST_SKIP() << "Ignoring this test as Highs 1.7+ returns unimplemented";
-  }
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<Model> model, Load23588());
   SolveParameters params = GetParam().base_parameters;
   params.solution_limit = 1;
diff --git a/ortools/math_opt/solver_tests/logical_constraint_tests.cc b/ortools/math_opt/solver_tests/logical_constraint_tests.cc
index 75577065ff..c4109b347e 100644
--- a/ortools/math_opt/solver_tests/logical_constraint_tests.cc
+++ b/ortools/math_opt/solver_tests/logical_constraint_tests.cc
@@ -13,11 +13,8 @@
 
 #include "ortools/math_opt/solver_tests/logical_constraint_tests.h"
 
-#include <limits>
 #include <memory>
 #include <ostream>
-#include <string>
-#include <type_traits>
 #include <utility>
 
 #include "absl/status/status.h"
@@ -97,9 +94,9 @@ TEST_P(SimpleLogicalConstraintTest, CanBuildSos1Model) {
   model.AddSos1Constraint({3.0 * x + 2.0}, {3.0});
   model.AddSos1Constraint({2.0 * x + 1.0}, {});
   if (GetParam().supports_sos1) {
-    EXPECT_OK(IncrementalSolver::New(&model, GetParam().solver_type, {}));
+    EXPECT_OK(NewIncrementalSolver(&model, GetParam().solver_type, {}));
   } else {
-    EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+    EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
                 StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                                absl::StatusCode::kUnimplemented),
                          HasSubstr("sos1 constraints")));
@@ -113,9 +110,9 @@ TEST_P(SimpleLogicalConstraintTest, CanBuildSos2Model) {
   model.AddSos2Constraint({3.0 * x + 2.0}, {3.0});
   model.AddSos2Constraint({2.0 * x + 1.0}, {});
   if (GetParam().supports_sos2) {
-    EXPECT_OK(IncrementalSolver::New(&model, GetParam().solver_type, {}));
+    EXPECT_OK(NewIncrementalSolver(&model, GetParam().solver_type, {}));
   } else {
-    EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+    EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
                 StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                                absl::StatusCode::kUnimplemented),
                          HasSubstr("sos2 constraints")));
@@ -294,7 +291,7 @@ TEST_P(IncrementalLogicalConstraintTest, LinearToSos1Update) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(3.0, {{x, 1.0}, {y, 1.0}})));
 
@@ -346,7 +343,7 @@ TEST_P(IncrementalLogicalConstraintTest, LinearToSos2Update) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(
       solver->Solve({.parameters = GetParam().parameters}),
       IsOkAndHolds(IsOptimalWithSolution(6.0, {{x, 1.0}, {y, 1.0}, {z, 1.0}})));
@@ -409,7 +406,7 @@ TEST_P(IncrementalLogicalConstraintTest, UpdateDeletesSos1Constraint) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(2.5, {{x, 0.25}, {y, 0.75}})));
 
@@ -455,7 +452,7 @@ TEST_P(IncrementalLogicalConstraintTest, UpdateDeletesSos2Constraint) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(
       solver->Solve({.parameters = GetParam().parameters}),
       IsOkAndHolds(IsOptimalWithSolution(4.5, {{x, 0.5}, {y, 1.0}, {z, 0.5}})));
@@ -503,7 +500,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(
                   3.0, {{x, 0.0}, {y, 1.0}, {z, 0.0}, {w, 1.0}})));
@@ -551,7 +548,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(
                   5.0, {{x, 0.0}, {y, 1.0}, {z, 1.0}, {w, 1.0}})));
@@ -600,7 +597,7 @@ TEST_P(IncrementalLogicalConstraintTest, InstanceWithSos1AndSos2AndDeletion) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(
       solver->Solve({.parameters = GetParam().parameters}),
       IsOkAndHolds(IsOptimalWithSolution(3.5, {{x, 1.0}, {y, 1.0}, {z, 0.0}})));
@@ -629,9 +626,9 @@ TEST_P(SimpleLogicalConstraintTest, CanBuildIndicatorModel) {
   model.AddIndicatorConstraint(x, y <= 0.5);
 
   if (GetParam().supports_indicator_constraints) {
-    EXPECT_OK(IncrementalSolver::New(&model, GetParam().solver_type, {}));
+    EXPECT_OK(NewIncrementalSolver(&model, GetParam().solver_type, {}));
   } else {
-    EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+    EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
                 StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                                absl::StatusCode::kUnimplemented),
                          HasSubstr("indicator constraints")));
@@ -712,7 +709,7 @@ TEST_P(SimpleLogicalConstraintTest, IndicatorWithRangedImpliedConstraint) {
   const Variable y = model.AddContinuousVariable(0.0, 1.0, "y");
   model.AddIndicatorConstraint(x, 0.25 <= y <= 0.75);
 
-  EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+  EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
               StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                              absl::StatusCode::kUnimplemented),
                        HasSubstr("ranged")));
@@ -805,7 +802,7 @@ TEST_P(IncrementalLogicalConstraintTest, LinearToIndicatorUpdate) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(2.0, {{x, 1.0}, {y, 1.0}})));
 
@@ -865,7 +862,7 @@ TEST_P(IncrementalLogicalConstraintTest, UpdateDeletesIndicatorConstraint) {
   const IndicatorConstraint c = model.AddIndicatorConstraint(x, y <= 0.5);
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.5, {{x, 1.0}, {y, 0.5}})));
 
@@ -911,7 +908,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.0, {{x, 1.0}})));
 
@@ -953,7 +950,7 @@ TEST_P(IncrementalLogicalConstraintTest, UpdateDeletesIndicatorVariable) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.5, {{x, 1.0}, {y, 0.5}})));
 
@@ -1002,7 +999,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(
       solver->Solve({.parameters = GetParam().parameters}),
       IsOkAndHolds(IsOptimalWithSolution(3.0, {{x, 0.0}, {y, 1.0}, {z, 1.0}})));
@@ -1032,7 +1029,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_OK(solver->Solve({.parameters = GetParam().parameters}));
 
   model.set_continuous(x);
@@ -1084,7 +1081,7 @@ TEST_P(IncrementalLogicalConstraintTest, UpdateChangesIndicatorVariableBound) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   EXPECT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(
                   1.6, {{x, 1.0}, {y, 0.0}, {u, 0.6}, {v, 1.0}})));
@@ -1143,7 +1140,7 @@ TEST_P(IncrementalLogicalConstraintTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_OK(solver->Solve({.parameters = GetParam().parameters}));
 
   model.set_upper_bound(x, 2.0);
diff --git a/ortools/math_opt/solver_tests/lp_incomplete_solve_tests.cc b/ortools/math_opt/solver_tests/lp_incomplete_solve_tests.cc
index f064abbc0d..5e215c55c3 100644
--- a/ortools/math_opt/solver_tests/lp_incomplete_solve_tests.cc
+++ b/ortools/math_opt/solver_tests/lp_incomplete_solve_tests.cc
@@ -19,7 +19,6 @@
 #include <memory>
 #include <optional>
 #include <ostream>
-#include <type_traits>
 #include <vector>
 
 #include "absl/strings/str_cat.h"
@@ -536,7 +535,7 @@ TEST_P(LpIncompleteSolveTest, PrimalSimplexAlgorithm) {
 
   ASSERT_OK_AND_ASSIGN(
       const std::unique_ptr<IncrementalSolver> incremental_solver,
-      IncrementalSolver::New(&model, TestedSolver()));
+      NewIncrementalSolver(&model, TestedSolver()));
   SolveArguments args;
   args.parameters.lp_algorithm = GetParam().lp_algorithm;
   if (GetParam().supports_presolve) {
@@ -705,7 +704,7 @@ TEST_P(LpIncompleteSolveTest, PrimalSimplexAlgorithmRanged) {
   LinearConstraint c = model.AddLinearConstraint(Sum(x) >= 1);
   ASSERT_OK_AND_ASSIGN(
       const std::unique_ptr<IncrementalSolver> incremental_solver,
-      IncrementalSolver::New(&model, TestedSolver()));
+      NewIncrementalSolver(&model, TestedSolver()));
 
   SolveArguments args;
   args.parameters.lp_algorithm = LPAlgorithm::kPrimalSimplex;
@@ -963,7 +962,7 @@ TEST_P(LpIncompleteSolveTest, DualSimplexAlgorithmIncrementalCut) {
   model.Maximize(Sum(x));
   ASSERT_OK_AND_ASSIGN(
       const std::unique_ptr<IncrementalSolver> incremental_solver,
-      IncrementalSolver::New(&model, TestedSolver()));
+      NewIncrementalSolver(&model, TestedSolver()));
 
   ASSERT_OK(incremental_solver->Solve());
 
@@ -1109,7 +1108,7 @@ TEST_P(LpIncompleteSolveTest, PhaseIDualSimplexAlgorithm) {
 
   ASSERT_OK_AND_ASSIGN(
       const std::unique_ptr<IncrementalSolver> incremental_solver,
-      IncrementalSolver::New(&model, TestedSolver()));
+      NewIncrementalSolver(&model, TestedSolver()));
   SolveArguments args;
   args.parameters.lp_algorithm = GetParam().lp_algorithm;
   if (GetParam().supports_presolve) {
diff --git a/ortools/math_opt/solver_tests/lp_initial_basis_tests.cc b/ortools/math_opt/solver_tests/lp_initial_basis_tests.cc
index e4d0c29698..1f93a9abbb 100644
--- a/ortools/math_opt/solver_tests/lp_initial_basis_tests.cc
+++ b/ortools/math_opt/solver_tests/lp_initial_basis_tests.cc
@@ -16,7 +16,6 @@
 #include <limits>
 #include <memory>
 #include <optional>
-#include <string>
 #include <vector>
 
 #include "absl/log/check.h"
@@ -53,7 +52,7 @@ SolveStats LpBasisStartTest::SolveWithWarmStart(
 SolveStats LpBasisStartTest::RoundTripSolve() {
   model_.Maximize(objective_expression_);
   const std::unique_ptr<IncrementalSolver> solver =
-      IncrementalSolver::New(&model_, TestedSolver()).value();
+      NewIncrementalSolver(&model_, TestedSolver()).value();
   const SolveResult max_result = solver->Solve({.parameters = params_}).value();
   CHECK_OK(max_result.termination.EnsureIsOptimal());
   ModelSolveParameters max_model_parameters;
diff --git a/ortools/math_opt/solver_tests/lp_tests.cc b/ortools/math_opt/solver_tests/lp_tests.cc
index dfae7371f5..7d0fda5f2b 100644
--- a/ortools/math_opt/solver_tests/lp_tests.cc
+++ b/ortools/math_opt/solver_tests/lp_tests.cc
@@ -81,7 +81,7 @@ IncrementalLpTest::IncrementalLpTest()
       y_3_(model_.AddContinuousVariable(0, 1, "y_3")),
       c_3_(model_.AddLinearConstraint(x_3_ + y_3_ <= 1.5, "c_3")) {
   model_.Maximize(0.1 + 3 * (x_1_ + x_2_ + x_3_) + 2 * (y_1_ + y_2_ + y_3_));
-  solver_ = IncrementalSolver::New(&model_, TestedSolver()).value();
+  solver_ = NewIncrementalSolver(&model_, TestedSolver()).value();
   const SolveResult first_solve = solver_->Solve().value();
   CHECK_OK(first_solve.termination.EnsureIsOptimal());
   CHECK_LE(std::abs(first_solve.objective_value() - 12.1), kTolerance);
@@ -953,7 +953,7 @@ TEST_P(SimpleLpTest, OptimalAfterInfeasible) {
   const SolveArguments arguments{.parameters = GetParam().parameters};
 
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
   EXPECT_THAT(solver->Solve(arguments),
               IsOkAndHolds(TerminatesWithOneOf(
                   {TerminationReason::kInfeasible,
@@ -975,7 +975,7 @@ TEST_P(SimpleLpTest, OptimalAfterUnbounded) {
   const SolveArguments arguments{.parameters = GetParam().parameters};
 
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
   EXPECT_THAT(solver->Solve(arguments),
               IsOkAndHolds(TerminatesWithOneOf(
                   {TerminationReason::kUnbounded,
diff --git a/ortools/math_opt/solver_tests/mip_tests.cc b/ortools/math_opt/solver_tests/mip_tests.cc
index 8c9cbd36e8..4d846f1d20 100644
--- a/ortools/math_opt/solver_tests/mip_tests.cc
+++ b/ortools/math_opt/solver_tests/mip_tests.cc
@@ -17,9 +17,6 @@
 #include <initializer_list>
 #include <limits>
 #include <ostream>
-#include <string>
-#include <type_traits>
-#include <utility>
 
 #include "absl/log/check.h"
 #include "absl/status/statusor.h"
@@ -65,7 +62,7 @@ IncrementalMipTest::IncrementalMipTest()
       y_(model_.AddIntegerVariable(0.0, 2.0, "y")),
       c_(model_.AddLinearConstraint(0 <= x_ + y_ <= 1.5, "c")) {
   model_.Maximize(3.0 * x_ + 2.0 * y_ + 0.1);
-  solver_ = IncrementalSolver::New(&model_, TestedSolver()).value();
+  solver_ = NewIncrementalSolver(&model_, TestedSolver()).value();
   const SolveResult first_solve = solver_->Solve().value();
   CHECK(first_solve.has_primal_feasible_solution());
   CHECK_LE(std::abs(first_solve.objective_value() - 3.6), kTolerance)
@@ -179,7 +176,7 @@ TEST_P(IncrementalMipTest, DISABLED_MakeContinuousWithNonIntegralBounds) {
   model.Maximize(x);
 
   ASSERT_OK_AND_ASSIGN(const auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
   ASSERT_THAT(solver->Solve(), IsOkAndHolds(IsOptimal(1.0)));
 
   // Switching to continuous should use the fractional bound. For solvers that
@@ -203,7 +200,7 @@ TEST_P(IncrementalMipTest, MakeIntegralWithNonIntegralBounds) {
   model.Maximize(x);
 
   ASSERT_OK_AND_ASSIGN(const auto solver,
-                       IncrementalSolver::New(&model, TestedSolver()));
+                       NewIncrementalSolver(&model, TestedSolver()));
   ASSERT_THAT(solver->Solve(), IsOkAndHolds(IsOptimal(1.5)));
 
   model.set_integer(x);
diff --git a/ortools/math_opt/solver_tests/multi_objective_tests.cc b/ortools/math_opt/solver_tests/multi_objective_tests.cc
index 88f9009e36..5b8f1394bb 100644
--- a/ortools/math_opt/solver_tests/multi_objective_tests.cc
+++ b/ortools/math_opt/solver_tests/multi_objective_tests.cc
@@ -13,7 +13,6 @@
 
 #include "ortools/math_opt/solver_tests/multi_objective_tests.h"
 
-#include <limits>
 #include <memory>
 #include <ostream>
 
@@ -81,9 +80,9 @@ TEST_P(SimpleMultiObjectiveTest, CanBuildMultiObjectiveModel) {
   model.AddMinimizationObjective(-3.0 * x + 2.0, /*priority=*/1);
 
   if (GetParam().supports_auxiliary_objectives) {
-    EXPECT_OK(IncrementalSolver::New(&model, GetParam().solver_type, {}));
+    EXPECT_OK(NewIncrementalSolver(&model, GetParam().solver_type, {}));
   } else {
-    EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+    EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
                 StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                                absl::StatusCode::kUnimplemented),
                          HasSubstr("multiple objectives")));
@@ -145,7 +144,7 @@ TEST_P(SimpleMultiObjectiveTest, PrimaryAndAuxiliaryObjectiveSharePriority) {
   Model model;
   model.set_objective_priority(model.primary_objective(), 1);
   model.AddAuxiliaryObjective(1);
-  EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+  EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
               StatusIs(absl::StatusCode::kInvalidArgument,
                        HasSubstr("repeated objective priority: 1")));
 }
@@ -157,7 +156,7 @@ TEST_P(SimpleMultiObjectiveTest, AuxiliaryObjectivesSharePriority) {
   Model model;
   model.AddAuxiliaryObjective(1);
   model.AddAuxiliaryObjective(1);
-  EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+  EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
               StatusIs(absl::StatusCode::kInvalidArgument,
                        HasSubstr("repeated objective priority: 1")));
 }
@@ -369,7 +368,7 @@ TEST_P(IncrementalMultiObjectiveTest, SingleToMultiObjectiveModel) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   // Since there are multiple optimal solutions we do not match against the
   // solution value.
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
@@ -438,7 +437,7 @@ TEST_P(IncrementalMultiObjectiveTest, AddObjectiveToMultiObjectiveModel) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   {
     ASSERT_OK_AND_ASSIGN(const SolveResult result,
                          solver->Solve({.parameters = GetParam().parameters}));
@@ -493,7 +492,7 @@ TEST_P(IncrementalMultiObjectiveTest, DeleteObjectiveFromMultiObjectiveModel) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   {
     ASSERT_OK_AND_ASSIGN(const SolveResult result,
                          solver->Solve({.parameters = GetParam().parameters}));
@@ -547,7 +546,7 @@ TEST_P(IncrementalMultiObjectiveTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   {
     ASSERT_OK_AND_ASSIGN(const SolveResult result,
                          solver->Solve({.parameters = GetParam().parameters}));
@@ -587,7 +586,7 @@ TEST_P(IncrementalMultiObjectiveTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   {
     ASSERT_OK_AND_ASSIGN(const SolveResult result,
                          solver->Solve({.parameters = GetParam().parameters}));
diff --git a/ortools/math_opt/solver_tests/qc_tests.cc b/ortools/math_opt/solver_tests/qc_tests.cc
index 7aed4374a4..8ebab7eec4 100644
--- a/ortools/math_opt/solver_tests/qc_tests.cc
+++ b/ortools/math_opt/solver_tests/qc_tests.cc
@@ -16,8 +16,6 @@
 #include <cmath>
 #include <memory>
 #include <ostream>
-#include <string>
-#include <type_traits>
 #include <utility>
 
 #include "absl/status/status.h"
@@ -94,10 +92,10 @@ TEST_P(SimpleQcTest, CanBuildQcModel) {
   UnivariateQcProblem qc_problem(GetParam().use_integer_variables);
   if (GetParam().supports_qc) {
     EXPECT_OK(
-        IncrementalSolver::New(&qc_problem.model, GetParam().solver_type, {}));
+        NewIncrementalSolver(&qc_problem.model, GetParam().solver_type, {}));
   } else {
     EXPECT_THAT(
-        IncrementalSolver::New(&qc_problem.model, GetParam().solver_type, {}),
+        NewIncrementalSolver(&qc_problem.model, GetParam().solver_type, {}),
         StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                        absl::StatusCode::kUnimplemented),
                  HasSubstr("quadratic constraints")));
@@ -180,7 +178,7 @@ TEST_P(IncrementalQcTest, LinearToQuadraticUpdate) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(2.0, {{x, 1.0}})));
 
@@ -237,7 +235,7 @@ TEST_P(IncrementalQcTest, UpdateDeletesQuadraticConstraint) {
   HalfEllipseProblem qc_problem(GetParam().use_integer_variables);
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&qc_problem.model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&qc_problem.model, GetParam().solver_type, {}));
   // We test that the solution is correct elsewhere.
   ASSERT_OK(solver->Solve({.parameters = GetParam().parameters}));
 
@@ -272,7 +270,7 @@ TEST_P(IncrementalQcTest, UpdateDeletesVariableInQuadraticConstraint) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&qc_problem.model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&qc_problem.model, GetParam().solver_type, {}));
   // We test that the solution is correct elsewhere.
   ASSERT_OK(solver->Solve({.parameters = GetParam().parameters}));
 
@@ -287,5 +285,162 @@ TEST_P(IncrementalQcTest, UpdateDeletesVariableInQuadraticConstraint) {
                                                  kTolerance)));
 }
 
+// Primal:
+//   min_{x} x
+//   s.t.
+//   Quadratic constraint:
+//       x^2 <= 1
+//
+// Optimal solution: x* = -1.
+//
+// Dual (go/mathopt-qcqp-dual):
+//   max_{mu, x, r}  mu + mu*x^2
+//   s.t.   mu*2*x + r  = 1
+//                  mu <= 0
+//                   r  = 0
+//
+// Optimal solution: x* = -1, mu* = -0.5.
+TEST_P(QcDualsTest, OnlyQuadraticConstraintLess) {
+  if (!GetParam().supports_qc) {
+    return;
+  }
+  Model model;
+  const Variable x = model.AddVariable();
+  const QuadraticConstraint mu = model.AddQuadraticConstraint(x * x <= 1);
+  model.Minimize(x);
+
+  ASSERT_OK_AND_ASSIGN(const SolveResult solve_result, SimpleSolve(model));
+  const double expected_objective_value = -1.0;
+  EXPECT_THAT(solve_result,
+              IsOptimalWithSolution(expected_objective_value, {{x, -1.0}}));
+  EXPECT_THAT(solve_result,
+              IsOptimalWithDualSolution(expected_objective_value, {},
+                                        {{mu, -0.5}}, {{x, 0.0}}));
+}
+
+// Primal:
+//   min_{x} x
+//   s.t.
+//   Quadratic constraint:
+//       -x^2 >= -1
+//
+// Optimal solution: x* = -1.
+//
+// Dual (go/mathopt-qcqp-dual):
+//   max_{mu, x, r}  -mu - mu*x^2
+//   s.t.  -mu*2*x + r  = 1
+//                  mu >= 0
+//                   r  = 0
+//
+// Optimal solution: x* = -1, mu* = 0.5.
+TEST_P(QcDualsTest, OnlyQuadraticConstraintGreater) {
+  if (!GetParam().supports_qc) {
+    return;
+  }
+  Model model;
+  const Variable x = model.AddVariable();
+  const QuadraticConstraint mu = model.AddQuadraticConstraint(-x * x >= -1);
+  model.Minimize(x);
+
+  ASSERT_OK_AND_ASSIGN(const SolveResult solve_result, SimpleSolve(model));
+  const double expected_objective_value = -1.0;
+  EXPECT_THAT(solve_result,
+              IsOptimalWithSolution(expected_objective_value, {{x, -1.0}}));
+  EXPECT_THAT(solve_result,
+              IsOptimalWithDualSolution(expected_objective_value, {},
+                                        {{mu, 0.5}}, {{x, 0.0}}));
+}
+
+// Primal:
+//   min_{x} x1^2 - 10 x1
+//   s.t.
+//   Quadratic constraints:
+//        x1^2 + x0 <= 2
+//   Linear constraints:
+//          x1 - x0 <= 0
+//         -x1 - x0 <= 0
+//
+// Optimal solution: x* = (1, 1).
+//
+// Dual (go/mathopt-qcqp-dual):
+//   max_{mu, x, y, r}  2*mu + mu*x1^2 - x1^2
+//   s.t.   -y0 - y1 + r0 + mu       = 0
+//           y0 - y1 + r1 + mu*2*x1  = 2*x1 - 10
+//                               mu <= 0
+//                               y0 <= 0
+//                               y1 <= 0
+//                               r0  = 0
+//                               r1  = 0
+//
+// Optimal solution: x* = (1, 1), mu* = -8/3, y = (-8/3, 0), r = (0, 0).
+TEST_P(QcDualsTest, QuadraticObjectiveAndLinearAndQuadraticConstraints) {
+  if (!GetParam().supports_qc) {
+    return;
+  }
+  Model model;
+  const Variable x0 = model.AddVariable();
+  const Variable x1 = model.AddVariable();
+  const LinearConstraint y0 = model.AddLinearConstraint(x1 - x0 <= 0.0);
+  const LinearConstraint y1 = model.AddLinearConstraint(-x1 - x0 <= 0.0);
+  const QuadraticConstraint mu =
+      model.AddQuadraticConstraint(x1 * x1 + x0 <= 2);
+  model.Minimize(x1 * x1 - 10.0 * x1);
+
+  ASSERT_OK_AND_ASSIGN(const SolveResult solve_result, SimpleSolve(model));
+  const double expected_objective_value = -9.0;
+  EXPECT_THAT(solve_result, IsOptimalWithSolution(expected_objective_value,
+                                                  {{x0, 1.0}, {x1, 1.0}}));
+  EXPECT_THAT(solve_result,
+              IsOptimalWithDualSolution(
+                  expected_objective_value, {{y0, -8.0 / 3.0}, {y1, 0.0}},
+                  {{mu, -8.0 / 3.0}}, {{x0, 0.0}, {x1, 0.0}}));
+}
+
+// Primal:
+//   max_{x} -x0^2 + 4x0
+//   s.t.
+//   Quadratic constraints:
+//        x0^2 + x1^2 + x2^2 <= 3
+//   Linear constraints:
+//          x1 = 1
+//   Variable bounds:
+//          x2 = 1
+//
+// Optimal solution: x* = (1, 1, 1).
+//
+// Dual (go/mathopt-qcqp-dual):
+//   min_{mu, x, y, r}  y + r2 + 3*mu + mu*(x0^2 + x1^2 + x2^2) + x0^2
+//   s.t.       r0     + mu*2*x0   = -2x0 + 4
+//              r1 + y + mu*2*x1   = 0
+//              r2 + mu*2*x2       = 0
+//                             mu >= 0
+//                             r0  = 0
+//                             r1  = 0
+//
+// Optimal solution: x* = (1, 1, 1), mu* = 1, y = -2, r = (0, 0, -2).
+TEST_P(QcDualsTest, MaxAndVariableBounds) {
+  if (!GetParam().supports_qc) {
+    return;
+  }
+  Model model;
+  const Variable x0 = model.AddVariable();
+  const Variable x1 = model.AddVariable();
+  const Variable x2 = model.AddContinuousVariable(1.0, 1.0);
+  const LinearConstraint y = model.AddLinearConstraint(x1 == 1.0);
+  const QuadraticConstraint mu =
+      model.AddQuadraticConstraint(x0 * x0 + x1 * x1 + x2 * x2 <= 3.0);
+  model.Maximize(-x0 * x0 + 4.0 * x0);
+
+  ASSERT_OK_AND_ASSIGN(const SolveResult solve_result, SimpleSolve(model));
+  const double expected_objective_value = 3.0;
+  EXPECT_THAT(solve_result,
+              IsOptimalWithSolution(expected_objective_value,
+                                    {{x0, 1.0}, {x1, 1.0}, {x2, 1.0}}));
+  EXPECT_THAT(solve_result,
+              IsOptimalWithDualSolution(expected_objective_value, {{y, -2.0}},
+                                        {{mu, 1.0}},
+                                        {{x0, 0.0}, {x1, 0.0}, {x2, -2.0}}));
+}
+
 }  // namespace
 }  // namespace operations_research::math_opt
diff --git a/ortools/math_opt/solver_tests/qc_tests.h b/ortools/math_opt/solver_tests/qc_tests.h
index b2d4fae4ce..4c04a8b5f9 100644
--- a/ortools/math_opt/solver_tests/qc_tests.h
+++ b/ortools/math_opt/solver_tests/qc_tests.h
@@ -83,6 +83,23 @@ class SimpleQcTest : public testing::TestWithParam<QcTestParameters> {
 //                                        use_integer_variables)));
 class IncrementalQcTest : public testing::TestWithParam<QcTestParameters> {};
 
+// A suite of unit tests focused on testing dual solutions from QC solvers.
+//
+// To use these tests, in file <solver>_test.cc, write:
+//   INSTANTIATE_TEST_SUITE_P(
+//       <Solver>QcDualsTest, QcDualsTest,
+//       testing::Values(QcTestParameters(SolverType::k<Solver>,
+//                                        parameters, supports_qc,
+//                                        supports_incremental_add_and_deletes,
+//                                        supports_incremental_variable_deletions,
+//                                        use_integer_variables)));
+class QcDualsTest : public testing::TestWithParam<QcTestParameters> {
+ protected:
+  absl::StatusOr<SolveResult> SimpleSolve(const Model& model) {
+    return Solve(model, GetParam().solver_type,
+                 {.parameters = GetParam().parameters});
+  }
+};
 }  // namespace operations_research::math_opt
 
 #endif  // OR_TOOLS_MATH_OPT_SOLVER_TESTS_QC_TESTS_H_
diff --git a/ortools/math_opt/solver_tests/qp_tests.cc b/ortools/math_opt/solver_tests/qp_tests.cc
index a161800118..74e13c9e85 100644
--- a/ortools/math_opt/solver_tests/qp_tests.cc
+++ b/ortools/math_opt/solver_tests/qp_tests.cc
@@ -202,9 +202,8 @@ TEST_P(IncrementalQpTest, EmptyUpdate) {
 
   UnivariateQpProblem qp_problem(GetParam().use_integer_variables);
 
-  ASSERT_OK_AND_ASSIGN(
-      const std::unique_ptr<IncrementalSolver> solver,
-      IncrementalSolver::New(&qp_problem.model, TestedSolver()));
+  ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                       NewIncrementalSolver(&qp_problem.model, TestedSolver()));
   ASSERT_OK_AND_ASSIGN(const SolveResult first_result,
                        solver->Solve({.parameters = GetParam().parameters}));
   ASSERT_THAT(first_result,
@@ -234,9 +233,8 @@ TEST_P(IncrementalQpTest, LinearToQuadraticUpdate) {
   // We remove the quadratic coefficient x * x from the objective, leaving an LP
   UnivariateQpProblem qp_problem(GetParam().use_integer_variables);
   qp_problem.model.set_objective_coefficient(qp_problem.x, qp_problem.x, 0);
-  ASSERT_OK_AND_ASSIGN(
-      const std::unique_ptr<IncrementalSolver> solver,
-      IncrementalSolver::New(&qp_problem.model, TestedSolver()));
+  ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                       NewIncrementalSolver(&qp_problem.model, TestedSolver()));
   ASSERT_OK_AND_ASSIGN(const SolveResult first_result,
                        solver->Solve({.parameters = GetParam().parameters}));
   ASSERT_THAT(first_result, IsOptimal(0.0625 - 0.5));
@@ -289,9 +287,8 @@ TEST_P(IncrementalQpTest, ModifyQuadraticObjective) {
 
   UnivariateQpProblem qp_problem(GetParam().use_integer_variables);
 
-  ASSERT_OK_AND_ASSIGN(
-      const std::unique_ptr<IncrementalSolver> solver,
-      IncrementalSolver::New(&qp_problem.model, TestedSolver()));
+  ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                       NewIncrementalSolver(&qp_problem.model, TestedSolver()));
   ASSERT_OK_AND_ASSIGN(const SolveResult first_result,
                        solver->Solve({.parameters = GetParam().parameters}));
   ASSERT_THAT(first_result,
@@ -327,9 +324,8 @@ TEST_P(IncrementalQpTest, DeleteVariable) {
 
   SimplexConstrainedQpProblem qp_problem(GetParam().use_integer_variables);
 
-  ASSERT_OK_AND_ASSIGN(
-      const std::unique_ptr<IncrementalSolver> solver,
-      IncrementalSolver::New(&qp_problem.model, TestedSolver()));
+  ASSERT_OK_AND_ASSIGN(const std::unique_ptr<IncrementalSolver> solver,
+                       NewIncrementalSolver(&qp_problem.model, TestedSolver()));
   ASSERT_OK_AND_ASSIGN(const SolveResult first_result,
                        solver->Solve({.parameters = GetParam().parameters}));
   ASSERT_THAT(first_result,
diff --git a/ortools/math_opt/solver_tests/second_order_cone_tests.cc b/ortools/math_opt/solver_tests/second_order_cone_tests.cc
index a11a5dab6e..1df7db1867 100644
--- a/ortools/math_opt/solver_tests/second_order_cone_tests.cc
+++ b/ortools/math_opt/solver_tests/second_order_cone_tests.cc
@@ -16,8 +16,6 @@
 #include <cmath>
 #include <memory>
 #include <ostream>
-#include <string>
-#include <type_traits>
 #include <utility>
 
 #include "absl/status/status.h"
@@ -75,9 +73,9 @@ TEST_P(SimpleSecondOrderConeTest, CanBuildSecondOrderConeModel) {
   const Variable x = model.AddContinuousVariable(0.0, 1.0, "x");
   model.AddSecondOrderConeConstraint({x}, 2 * x);
   if (GetParam().supports_soc_constraints) {
-    EXPECT_OK(IncrementalSolver::New(&model, GetParam().solver_type, {}));
+    EXPECT_OK(NewIncrementalSolver(&model, GetParam().solver_type, {}));
   } else {
-    EXPECT_THAT(IncrementalSolver::New(&model, GetParam().solver_type, {}),
+    EXPECT_THAT(NewIncrementalSolver(&model, GetParam().solver_type, {}),
                 StatusIs(AnyOf(absl::StatusCode::kInvalidArgument,
                                absl::StatusCode::kUnimplemented),
                          HasSubstr("second-order cone constraints")));
@@ -179,7 +177,7 @@ TEST_P(IncrementalSecondOrderConeTest, LinearToSecondOrderConeUpdate) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.5, {{x, 0.5}, {y, 1.0}})));
 
@@ -244,7 +242,7 @@ TEST_P(IncrementalSecondOrderConeTest, UpdateDeletesSecondOrderConeConstraint) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.0, {{x, 0.5}, {y, 0.5}},
                                                  kTolerance)));
@@ -286,7 +284,7 @@ TEST_P(IncrementalSecondOrderConeTest, UpdateDeletesUpperBoundingVariable) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(2.0, {{x, 1.0}, {y, 1.0}})));
 
@@ -329,7 +327,7 @@ TEST_P(IncrementalSecondOrderConeTest,
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(3.0, {{x, 2.0}, {y, 1.0}})));
 
@@ -371,7 +369,7 @@ TEST_P(IncrementalSecondOrderConeTest, UpdateDeletesVariableThatIsAnArgument) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(1.0, {{x, 1.0}, {y, 1.0}})));
 
@@ -413,7 +411,7 @@ TEST_P(IncrementalSecondOrderConeTest, UpdateDeletesVariableInAnArgument) {
 
   ASSERT_OK_AND_ASSIGN(
       const auto solver,
-      IncrementalSolver::New(&model, GetParam().solver_type, {}));
+      NewIncrementalSolver(&model, GetParam().solver_type, {}));
   ASSERT_THAT(solver->Solve({.parameters = GetParam().parameters}),
               IsOkAndHolds(IsOptimalWithSolution(2.0, {{x, 1.0}, {y, 2.0}})));
 
diff --git a/ortools/math_opt/solver_tests/status_tests.cc b/ortools/math_opt/solver_tests/status_tests.cc
index b973db956f..a3a9e660b4 100644
--- a/ortools/math_opt/solver_tests/status_tests.cc
+++ b/ortools/math_opt/solver_tests/status_tests.cc
@@ -99,9 +99,6 @@ TEST_P(StatusTest, PrimalAndDualInfeasible) {
         << "Ignoring this test as GLPK gets stuck in presolve for IP's "
            "with a primal-dual infeasible LP relaxation.";
   }
-  if (GetParam().solver_type == SolverType::kHighs) {
-    GTEST_SKIP() << "Ignoring this test as Highs 1.7+ returns error.";
-  }
 
   Model model;
   const Variable x1 =
@@ -337,7 +334,7 @@ TEST_P(StatusTest, InfeasibleIpWithPrimalDualFeasibleRelaxation2) {
     EXPECT_THAT(result, TerminatesWith(TerminationReason::kInfeasible));
     // Result validators imply primal problem status is kInfeasible.
     EXPECT_THAT(result.termination.problem_status,
-        Not(DualStatusIs(FeasibilityStatus::kInfeasible)));
+                Not(DualStatusIs(FeasibilityStatus::kInfeasible)));
   }
 }
 
@@ -376,9 +373,6 @@ TEST_P(StatusTest, IncompleteIpSolve) {
     GTEST_SKIP() << "Ignoring this test as it is an IP-only test and requires "
                     "support for node_limit.";
   }
-  if (GetParam().solver_type == SolverType::kHighs) {
-    GTEST_SKIP() << "Ignoring this test as Highs 1.7+ returns MODEL_INVALID";
-  }
   ASSERT_OK_AND_ASSIGN(const std::unique_ptr<const Model> model, Load23588());
   SolveArguments args = {
       .parameters = {.enable_output = true, .node_limit = 1}};
diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel
index e7e8054a97..e5ee020d12 100644
--- a/ortools/math_opt/solvers/BUILD.bazel
+++ b/ortools/math_opt/solvers/BUILD.bazel
@@ -313,7 +313,6 @@ cc_test(
     name = "cp_sat_solver_test",
     srcs = ["cp_sat_solver_test.cc"],
     shard_count = 10,
-    timeout = "eternal",
     deps = [
         ":cp_sat_solver",
         "//ortools/base:gmock_main",
diff --git a/ortools/math_opt/solvers/cp_sat_solver_test.cc b/ortools/math_opt/solvers/cp_sat_solver_test.cc
index 03a90d427e..d701648f52 100644
--- a/ortools/math_opt/solvers/cp_sat_solver_test.cc
+++ b/ortools/math_opt/solvers/cp_sat_solver_test.cc
@@ -112,6 +112,7 @@ INSTANTIATE_TEST_SUITE_P(CpSatSimpleQcTest, SimpleQcTest,
                          ValuesIn(GetCpSatQcTestParameters()));
 INSTANTIATE_TEST_SUITE_P(CpSatIncrementalQcTest, IncrementalQcTest,
                          ValuesIn(GetCpSatQcTestParameters()));
+GTEST_ALLOW_UNINSTANTIATED_PARAMETERIZED_TEST(QcDualsTest);
 
 SecondOrderConeTestParameters GetCpSatSecondOrderConeTestParameters() {
   return SecondOrderConeTestParameters(
diff --git a/ortools/math_opt/solvers/glop_solver_test.cc b/ortools/math_opt/solvers/glop_solver_test.cc
index 5c75f795e7..1957e8aa6b 100644
--- a/ortools/math_opt/solvers/glop_solver_test.cc
+++ b/ortools/math_opt/solvers/glop_solver_test.cc
@@ -145,6 +145,7 @@ INSTANTIATE_TEST_SUITE_P(GlopSimpleQcTest, SimpleQcTest,
                          testing::Values(GetGlopQcTestParameters()));
 INSTANTIATE_TEST_SUITE_P(GlopIncrementalQcTest, IncrementalQcTest,
                          testing::Values(GetGlopQcTestParameters()));
+GTEST_ALLOW_UNINSTANTIATED_PARAMETERIZED_TEST(QcDualsTest);
 
 SecondOrderConeTestParameters GetGlopSecondOrderConeTestParameters() {
   return SecondOrderConeTestParameters(
diff --git a/ortools/math_opt/solvers/glpk_solver_test.cc b/ortools/math_opt/solvers/glpk_solver_test.cc
index 3e6cffacbe..07a2c770ff 100644
--- a/ortools/math_opt/solvers/glpk_solver_test.cc
+++ b/ortools/math_opt/solvers/glpk_solver_test.cc
@@ -259,6 +259,7 @@ INSTANTIATE_TEST_SUITE_P(GlpkSimpleQcTest, SimpleQcTest,
                          ValuesIn(GetGlpkQcTestParameters()));
 INSTANTIATE_TEST_SUITE_P(GlpkIncrementalQcTest, IncrementalQcTest,
                          ValuesIn(GetGlpkQcTestParameters()));
+GTEST_ALLOW_UNINSTANTIATED_PARAMETERIZED_TEST(QcDualsTest);
 
 SecondOrderConeTestParameters GetGlpkSecondOrderConeTestParameters() {
   return SecondOrderConeTestParameters(
@@ -423,7 +424,7 @@ TEST(GlpkSolverDeathTest, DestroySolverFromAnotherThread) {
   Model model("model");
 
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<IncrementalSolver> incremental_solver,
-                       IncrementalSolver::New(&model, SolverType::kGlpk));
+                       NewIncrementalSolver(&model, SolverType::kGlpk));
 
 #if 0
   EXPECT_DEATH_IF_SUPPORTED(
@@ -450,7 +451,7 @@ TEST(GlpkSolverTest, SolveFromAnotherThread) {
   model.Maximize(x + y);
 
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<IncrementalSolver> incremental_solver,
-                       IncrementalSolver::New(&model, SolverType::kGlpk));
+                       NewIncrementalSolver(&model, SolverType::kGlpk));
 
   absl::StatusOr<SolveResult> solve_result_or;
   std::thread([&]() { solve_result_or = incremental_solver->Solve(); }).join();
@@ -466,7 +467,7 @@ TEST(GlpkSolverTest, UpdateFromAnotherThread) {
   model.Maximize(x + y);
 
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<IncrementalSolver> incremental_solver,
-                       IncrementalSolver::New(&model, SolverType::kGlpk));
+                       NewIncrementalSolver(&model, SolverType::kGlpk));
 
   model.set_lower_bound(x, 1.2);
 
@@ -486,7 +487,7 @@ TEST(GlpkSolverTest, FailedUpdateFromAnotherThread) {
   model.Maximize(x + y);
 
   ASSERT_OK_AND_ASSIGN(std::unique_ptr<IncrementalSolver> incremental_solver,
-                       IncrementalSolver::New(&model, SolverType::kGlpk));
+                       NewIncrementalSolver(&model, SolverType::kGlpk));
 
   // Quadratic objectives are not supported by GLPK.
   model.Maximize(x * x);
diff --git a/ortools/math_opt/solvers/gscip_solver_test.cc b/ortools/math_opt/solvers/gscip_solver_test.cc
index 4d150f4cb4..4deb0a3815 100644
--- a/ortools/math_opt/solvers/gscip_solver_test.cc
+++ b/ortools/math_opt/solvers/gscip_solver_test.cc
@@ -119,6 +119,7 @@ INSTANTIATE_TEST_SUITE_P(GscipSimpleQcTest, SimpleQcTest,
                          ValuesIn(GetGscipQcTestParameters()));
 INSTANTIATE_TEST_SUITE_P(GscipIncrementalQcTest, IncrementalQcTest,
                          ValuesIn(GetGscipQcTestParameters()));
+GTEST_ALLOW_UNINSTANTIATED_PARAMETERIZED_TEST(QcDualsTest);
 
 SecondOrderConeTestParameters GetGscipSecondOrderConeTestParameters() {
   return SecondOrderConeTestParameters(
@@ -343,7 +344,7 @@ TEST(GScipSolverTest, UpdatingLowerBoundNotAllowedOnBinaryVariables) {
   Model model;
   const Variable x = model.AddBinaryVariable("x");
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, SolverType::kGscip, {}));
+                       NewIncrementalSolver(&model, SolverType::kGscip, {}));
 
   model.set_lower_bound(x, -1.0);
   EXPECT_THAT(solver->Update(), IsOkAndHolds(Not(DidUpdate())));
@@ -353,7 +354,7 @@ TEST(GScipSolverTest, UpdatingUpperBoundNotAllowedOnBinaryVariables) {
   Model model;
   const Variable x = model.AddBinaryVariable("x");
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, SolverType::kGscip, {}));
+                       NewIncrementalSolver(&model, SolverType::kGscip, {}));
 
   model.set_upper_bound(x, 2.0);
   EXPECT_THAT(solver->Update(), IsOkAndHolds(Not(DidUpdate())));
@@ -364,7 +365,7 @@ TEST(GScipSolverTest, UpdatingLowerBoundNotAllowedOnImplicitBinaryVariables) {
   // This will be silently converted to a binary variable in SCIP.
   const Variable y = model.AddIntegerVariable(0.0, 1.0, "y");
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, SolverType::kGscip, {}));
+                       NewIncrementalSolver(&model, SolverType::kGscip, {}));
 
   model.set_lower_bound(y, -1.0);
   EXPECT_THAT(solver->Update(), IsOkAndHolds(Not(DidUpdate())));
@@ -375,7 +376,7 @@ TEST(GScipSolverTest, UpdatingUpperBoundNotAllowedOnImplicitBinaryVariables) {
   // This will be silently converted to a binary variable in SCIP.
   const Variable y = model.AddIntegerVariable(0.0, 1.0, "y");
   ASSERT_OK_AND_ASSIGN(auto solver,
-                       IncrementalSolver::New(&model, SolverType::kGscip, {}));
+                       NewIncrementalSolver(&model, SolverType::kGscip, {}));
 
   model.set_upper_bound(y, 2.0);
   EXPECT_THAT(solver->Update(), IsOkAndHolds(Not(DidUpdate())));
diff --git a/ortools/math_opt/solvers/gurobi/BUILD.bazel b/ortools/math_opt/solvers/gurobi/BUILD.bazel
index 75b52fb180..32f8f39ffd 100644
--- a/ortools/math_opt/solvers/gurobi/BUILD.bazel
+++ b/ortools/math_opt/solvers/gurobi/BUILD.bazel
@@ -11,12 +11,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-package(
-    default_visibility = [
-        "//ortools/math_opt:__subpackages__",
-    ],
-)
-
 cc_library(
     name = "g_gurobi",
     srcs = [
diff --git a/ortools/math_opt/solvers/gurobi_init_arguments.h b/ortools/math_opt/solvers/gurobi_init_arguments.h
index 2353decb16..2b9adacdb0 100644
--- a/ortools/math_opt/solvers/gurobi_init_arguments.h
+++ b/ortools/math_opt/solvers/gurobi_init_arguments.h
@@ -63,11 +63,11 @@ namespace math_opt {
 //
 //   ASSIGN_OR_RETURN(
 //       const std::unique_ptr<IncrementalSolver> incremental_solve_1,
-//       IncrementalSolver::New(model, SOLVER_TYPE_GUROBI,
+//       NewIncrementalSolver(model, SOLVER_TYPE_GUROBI,
 //                              SolverInitArguments(gurobi_args)));
 //   ASSIGN_OR_RETURN(
 //       const std::unique_ptr<IncrementalSolver> incremental_solve_2,
-//       IncrementalSolver::New(model, SOLVER_TYPE_GUROBI,
+//       NewIncrementalSolver(model, SOLVER_TYPE_GUROBI,
 //                              SolverInitArguments(gurobi_args)));
 //
 //   ASSIGN_OR_RETURN(const SolveResult result_1, incremental_solve_1->Solve());
diff --git a/ortools/math_opt/solvers/gurobi_solver.cc b/ortools/math_opt/solvers/gurobi_solver.cc
index b4b4525f2b..d83743c509 100644
--- a/ortools/math_opt/solvers/gurobi_solver.cc
+++ b/ortools/math_opt/solvers/gurobi_solver.cc
@@ -1393,7 +1393,7 @@ absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetLpSolution(
   ASSIGN_OR_RETURN(auto primal_solution_and_claim,
                    GetConvexPrimalSolutionIfAvailable(model_parameters));
   ASSIGN_OR_RETURN(auto dual_solution_and_claim,
-                   GetLpDualSolutionIfAvailable(model_parameters));
+                   GetConvexDualSolutionIfAvailable(model_parameters));
   ASSIGN_OR_RETURN(auto basis, GetBasisIfAvailable());
   const SolutionClaims solution_claims = {
       .primal_feasible_solution_exists =
@@ -1423,7 +1423,7 @@ absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetLpSolution(
 }
 
 absl::StatusOr<GurobiSolver::SolutionAndClaim<DualSolutionProto>>
-GurobiSolver::GetLpDualSolutionIfAvailable(
+GurobiSolver::GetConvexDualSolutionIfAvailable(
     const ModelSolveParametersProto& model_parameters) {
   if (!gurobi_->IsAttrAvailable(GRB_DBL_ATTR_PI) ||
       !gurobi_->IsAttrAvailable(GRB_DBL_ATTR_RC)) {
@@ -1450,6 +1450,17 @@ GurobiSolver::GetLpDualSolutionIfAvailable(
                                    *dual_solution.mutable_reduced_costs(),
                                    model_parameters.reduced_costs_filter());
 
+  if (!quadratic_constraints_map_.empty() &&
+      gurobi_->IsAttrAvailable(GRB_DBL_ATTR_QCPI)) {
+    ASSIGN_OR_RETURN(
+        const std::vector<double> grb_quad_constraint_duals,
+        gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_QCPI, num_gurobi_quad_cons_));
+    GurobiVectorToSparseDoubleVector(
+        grb_quad_constraint_duals, quadratic_constraints_map_,
+        *dual_solution.mutable_quadratic_dual_values(),
+        model_parameters.quadratic_dual_values_filter());
+  }
+
   ASSIGN_OR_RETURN(const int grb_termination,
                    gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
   if (grb_termination == GRB_OPTIMAL &&
@@ -1532,7 +1543,7 @@ absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetQpSolution(
   // TODO(b/225189115): Expand QpDualsTest to check maximization problems and
   // other edge cases.
   ASSIGN_OR_RETURN((auto [dual_solution, found_dual_feasible_solution]),
-                   GetLpDualSolutionIfAvailable(model_parameters));
+                   GetConvexDualSolutionIfAvailable(model_parameters));
   // TODO(b/280353996): we want to extract the basis here (when we solve via
   // simplex), but the existing code extracts a basis which fails our validator.
 
@@ -1559,11 +1570,8 @@ absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetQcpSolution(
     const ModelSolveParametersProto& model_parameters) {
   ASSIGN_OR_RETURN((auto [primal_solution, found_primal_feasible_solution]),
                    GetConvexPrimalSolutionIfAvailable(model_parameters));
-  // TODO(b/227217735): Expose duals. Note that the user must set the QCPDual
-  //                    parameter for Gurobi to return any dual values.
-  ASSIGN_OR_RETURN((auto [_, found_dual_feasible_solution]),
-                   GetLpDualSolutionIfAvailable(model_parameters));
-
+  ASSIGN_OR_RETURN((auto [dual_solution, found_dual_feasible_solution]),
+                   GetConvexDualSolutionIfAvailable(model_parameters));
   ASSIGN_OR_RETURN(const int grb_termination,
                    gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
   // By default, Gurobi will not return duals for optimally solved QCPs.
@@ -1574,9 +1582,15 @@ absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetQcpSolution(
           found_dual_feasible_solution || proven_feasible};
 
   SolutionsAndClaims solution_and_claims{.solution_claims = solution_claims};
-  if (primal_solution.has_value()) {
-    *solution_and_claims.solutions.emplace_back().mutable_primal_solution() =
-        *std::move(primal_solution);
+  if (primal_solution.has_value() || dual_solution.has_value()) {
+    SolutionProto& solution =
+        solution_and_claims.solutions.emplace_back(SolutionProto());
+    if (primal_solution.has_value()) {
+      *solution.mutable_primal_solution() = *std::move(primal_solution);
+    }
+    if (dual_solution.has_value()) {
+      *solution.mutable_dual_solution() = *std::move(dual_solution);
+    }
   }
   return solution_and_claims;
 }
diff --git a/ortools/math_opt/solvers/gurobi_solver.h b/ortools/math_opt/solvers/gurobi_solver.h
index 3f2eb76185..eb8a695bc2 100644
--- a/ortools/math_opt/solvers/gurobi_solver.h
+++ b/ortools/math_opt/solvers/gurobi_solver.h
@@ -256,7 +256,7 @@ class GurobiSolver : public SolverInterface {
   GetConvexPrimalSolutionIfAvailable(
       const ModelSolveParametersProto& model_parameters);
   absl::StatusOr<SolutionAndClaim<DualSolutionProto>>
-  GetLpDualSolutionIfAvailable(
+  GetConvexDualSolutionIfAvailable(
       const ModelSolveParametersProto& model_parameters);
   absl::StatusOr<std::optional<BasisProto>> GetBasisIfAvailable();
 
diff --git a/ortools/math_opt/solvers/pdlp_solver_test.cc b/ortools/math_opt/solvers/pdlp_solver_test.cc
index 4b5035fa9d..9b82eba3dd 100644
--- a/ortools/math_opt/solvers/pdlp_solver_test.cc
+++ b/ortools/math_opt/solvers/pdlp_solver_test.cc
@@ -111,6 +111,7 @@ INSTANTIATE_TEST_SUITE_P(PdlpSimpleQcTest, SimpleQcTest,
                          testing::Values(GetPdlpQcTestParameters()));
 INSTANTIATE_TEST_SUITE_P(PdlpIncrementalQcTest, IncrementalQcTest,
                          testing::Values(GetPdlpQcTestParameters()));
+GTEST_ALLOW_UNINSTANTIATED_PARAMETERIZED_TEST(QcDualsTest);
 
 SecondOrderConeTestParameters GetPdlpSecondOrderConeTestParameters() {
   return SecondOrderConeTestParameters(
diff --git a/ortools/math_opt/storage/BUILD.bazel b/ortools/math_opt/storage/BUILD.bazel
index 01f900655b..cb85a81422 100644
--- a/ortools/math_opt/storage/BUILD.bazel
+++ b/ortools/math_opt/storage/BUILD.bazel
@@ -180,6 +180,7 @@ cc_library(
         ":atomic_constraint_storage",
         ":iterators",
         ":linear_constraint_storage",
+        ":model_storage_types",
         ":objective_storage",
         ":sparse_matrix",
         ":update_trackers",
@@ -202,7 +203,6 @@ cc_library(
         "@com_google_absl//absl/container:flat_hash_map",
         "@com_google_absl//absl/container:flat_hash_set",
         "@com_google_absl//absl/log:check",
-        "@com_google_absl//absl/meta:type_traits",
         "@com_google_absl//absl/status",
         "@com_google_absl//absl/status:statusor",
         "@com_google_absl//absl/strings",
diff --git a/ortools/math_opt/storage/model_storage.h b/ortools/math_opt/storage/model_storage.h
index c68e108ced..5574b6f5e1 100644
--- a/ortools/math_opt/storage/model_storage.h
+++ b/ortools/math_opt/storage/model_storage.h
@@ -20,18 +20,14 @@
 #include <optional>
 #include <string>
 #include <tuple>
-#include <utility>
 #include <vector>
 
 #include "absl/container/flat_hash_map.h"
 #include "absl/container/flat_hash_set.h"
 #include "absl/log/check.h"
-#include "absl/meta/type_traits.h"
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
 #include "absl/strings/string_view.h"
-#include "ortools/base/map_util.h"
-#include "ortools/base/strong_int.h"
 #include "ortools/math_opt/constraints/indicator/storage.h"  // IWYU pragma: export
 #include "ortools/math_opt/constraints/quadratic/storage.h"  // IWYU pragma: export
 #include "ortools/math_opt/constraints/second_order_cone/storage.h"
@@ -42,8 +38,8 @@
 #include "ortools/math_opt/storage/atomic_constraint_storage.h"  // IWYU pragma: export
 #include "ortools/math_opt/storage/iterators.h"
 #include "ortools/math_opt/storage/linear_constraint_storage.h"
+#include "ortools/math_opt/storage/model_storage_types.h"
 #include "ortools/math_opt/storage/objective_storage.h"
-#include "ortools/math_opt/storage/sparse_matrix.h"
 #include "ortools/math_opt/storage/update_trackers.h"
 #include "ortools/math_opt/storage/variable_storage.h"
 
@@ -395,6 +391,8 @@ class ModelStorage {
   inline const absl::flat_hash_map<VariableId, double>& linear_objective(
       ObjectiveId id) const;
 
+  inline int64_t num_linear_objective_terms(ObjectiveId id) const;
+
   inline int64_t num_quadratic_objective_terms(ObjectiveId id) const;
 
   // The variable pairs with nonzero quadratic objective coefficients. The keys
@@ -1012,6 +1010,10 @@ const absl::flat_hash_map<VariableId, double>& ModelStorage::linear_objective(
   return objectives_.linear_terms(id);
 }
 
+int64_t ModelStorage::num_linear_objective_terms(const ObjectiveId id) const {
+  return objectives_.linear_terms(id).size();
+}
+
 int64_t ModelStorage::num_quadratic_objective_terms(
     const ObjectiveId id) const {
   return objectives_.quadratic_terms(id).nonzeros();
diff --git a/ortools/math_opt/storage/objective_storage.h b/ortools/math_opt/storage/objective_storage.h
index acaaba5968..f20cbc388b 100644
--- a/ortools/math_opt/storage/objective_storage.h
+++ b/ortools/math_opt/storage/objective_storage.h
@@ -393,12 +393,12 @@ void ObjectiveStorage::Clear(const ObjectiveId id,
     if (!diff.objective_tracked(id)) {
       continue;
     }
-    for (const auto [var, _] : data.linear_terms.terms()) {
+    for (const auto& [var, _] : data.linear_terms.terms()) {
       if (var < diff.variable_checkpoint) {
         diff.objective_diffs[id].linear_coefficients.insert(var);
       }
     }
-    for (const auto [v1, v2, _] : data.quadratic_terms.Terms()) {
+    for (const auto& [v1, v2, _] : data.quadratic_terms.Terms()) {
       if (v2 < diff.variable_checkpoint) {  // v1 <= v2 is implied
         diff.objective_diffs[id].quadratic_coefficients.insert({v1, v2});
       }
diff --git a/ortools/math_opt/storage/sparse_matrix.cc b/ortools/math_opt/storage/sparse_matrix.cc
index e23a3a0c9a..efb94d1184 100644
--- a/ortools/math_opt/storage/sparse_matrix.cc
+++ b/ortools/math_opt/storage/sparse_matrix.cc
@@ -151,7 +151,7 @@ SparseDoubleMatrixProto SparseSymmetricMatrix::Proto() const {
     // supported.
     std::vector<std::pair<VariableId, double>> related = Terms(v);
     absl::c_sort(related);
-    for (const auto [other, coef] : related) {
+    for (const auto& [other, coef] : related) {
       if (v <= other) {
         result.add_row_ids(v.value());
         result.add_column_ids(other.value());
@@ -167,7 +167,7 @@ SparseDoubleMatrixProto SparseSymmetricMatrix::Update(
     const absl::Span<const VariableId> new_variables,
     const absl::flat_hash_set<std::pair<VariableId, VariableId>>& dirty) const {
   std::vector<std::tuple<VariableId, VariableId, double>> updates;
-  for (const std::pair<VariableId, VariableId> pair : dirty) {
+  for (const std::pair<VariableId, VariableId>& pair : dirty) {
     // If either variable has been deleted, don't add it.
     if (deleted_variables.contains(pair.first) ||
         deleted_variables.contains(pair.second)) {
@@ -179,7 +179,7 @@ SparseDoubleMatrixProto SparseSymmetricMatrix::Update(
   for (const VariableId v : new_variables) {
     if (related_variables_.contains(v)) {
       // TODO(b/233630053): do not allocate here.
-      for (const auto [other, coef] : Terms(v)) {
+      for (const auto& [other, coef] : Terms(v)) {
         if (v <= other) {
           updates.push_back({v, other, coef});
         } else if (new_variables.empty() || other < new_variables[0]) {
diff --git a/ortools/math_opt/tools/BUILD.bazel b/ortools/math_opt/tools/BUILD.bazel
index f3aa9c4afc..55c1f2f460 100644
--- a/ortools/math_opt/tools/BUILD.bazel
+++ b/ortools/math_opt/tools/BUILD.bazel
@@ -74,6 +74,7 @@ cc_library(
         "//ortools/math_opt:model_cc_proto",
         "//ortools/math_opt:model_parameters_cc_proto",
         "//ortools/math_opt/io:lp_converter",
+        "//ortools/math_opt/io:lp_parser",
         "//ortools/math_opt/io:mps_converter",
         "//ortools/math_opt/io:proto_converter",
         "@com_google_absl//absl/algorithm:container",
diff --git a/ortools/math_opt/tools/file_format_flags.cc b/ortools/math_opt/tools/file_format_flags.cc
index c301e1ab53..c880371215 100644
--- a/ortools/math_opt/tools/file_format_flags.cc
+++ b/ortools/math_opt/tools/file_format_flags.cc
@@ -34,6 +34,7 @@
 #include "ortools/base/status_macros.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
 #include "ortools/math_opt/io/lp_converter.h"
+#include "ortools/math_opt/io/lp_parser.h"
 #include "ortools/math_opt/io/mps_converter.h"
 #include "ortools/math_opt/io/proto_converter.h"
 #include "ortools/math_opt/model.pb.h"
@@ -225,8 +226,10 @@ ReadModel(const absl::string_view file_path, const FileFormat format) {
       return std::make_pair(std::move(model), std::nullopt);
     }
     case FileFormat::kLP: {
-      return absl::UnimplementedError(
-          "MathOpt does not yet support reading .lp files; only writing them");
+      ASSIGN_OR_RETURN(const std::string lp_data,
+                       file::GetContents(file_path, file::Defaults()));
+      ASSIGN_OR_RETURN(ModelProto model, ModelProtoFromLp(lp_data));
+      return std::make_pair(std::move(model), std::nullopt);
     }
   }
 }