diff --git a/makefiles/Makefile.gen.mk b/makefiles/Makefile.gen.mk
index b893e4b69d..94f8d096ac 100644
--- a/makefiles/Makefile.gen.mk
+++ b/makefiles/Makefile.gen.mk
@@ -2637,11 +2637,14 @@ $(OBJ_DIR)/bop/bop_parameters.pb.$O: \
 
 LP_DEPS = \
  $(SRC_DIR)/ortools/linear_solver/glop_utils.h \
+ $(SRC_DIR)/ortools/linear_solver/gurobi_proto_solver.h \
  $(SRC_DIR)/ortools/linear_solver/linear_expr.h \
  $(SRC_DIR)/ortools/linear_solver/linear_solver.h \
  $(SRC_DIR)/ortools/linear_solver/model_exporter.h \
  $(SRC_DIR)/ortools/linear_solver/model_exporter_swig_helper.h \
  $(SRC_DIR)/ortools/linear_solver/model_validator.h \
+ $(SRC_DIR)/ortools/linear_solver/scip_helper_macros.h \
+ $(SRC_DIR)/ortools/linear_solver/scip_proto_solver.h \
  $(GEN_DIR)/ortools/linear_solver/linear_solver.pb.h
 
 LP_LIB_OBJS = \
@@ -2653,6 +2656,7 @@ LP_LIB_OBJS = \
  $(OBJ_DIR)/linear_solver/glop_utils.$O \
  $(OBJ_DIR)/linear_solver/glpk_interface.$O \
  $(OBJ_DIR)/linear_solver/gurobi_interface.$O \
+ $(OBJ_DIR)/linear_solver/gurobi_proto_solver.$O \
  $(OBJ_DIR)/linear_solver/linear_expr.$O \
  $(OBJ_DIR)/linear_solver/linear_solver.$O \
  $(OBJ_DIR)/linear_solver/model_exporter.$O \
@@ -2767,6 +2771,10 @@ objs/linear_solver/gurobi_interface.$O: \
  ortools/linear_solver/gurobi_interface.cc | $(OBJ_DIR)/linear_solver
 	$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Slinear_solver$Sgurobi_interface.cc $(OBJ_OUT)$(OBJ_DIR)$Slinear_solver$Sgurobi_interface.$O
 
+objs/linear_solver/gurobi_proto_solver.$O: \
+ ortools/linear_solver/gurobi_proto_solver.cc | $(OBJ_DIR)/linear_solver
+	$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Slinear_solver$Sgurobi_proto_solver.cc $(OBJ_OUT)$(OBJ_DIR)$Slinear_solver$Sgurobi_proto_solver.$O
+
 objs/linear_solver/linear_expr.$O: ortools/linear_solver/linear_expr.cc \
  ortools/linear_solver/linear_expr.h ortools/base/logging.h \
  ortools/base/integral_types.h ortools/base/macros.h \
diff --git a/ortools/linear_solver/gurobi_interface.cc b/ortools/linear_solver/gurobi_interface.cc
index f359b1155e..305f0f02d4 100644
--- a/ortools/linear_solver/gurobi_interface.cc
+++ b/ortools/linear_solver/gurobi_interface.cc
@@ -22,6 +22,7 @@
 #include <utility>
 #include <vector>
 
+#include "absl/strings/match.h"
 #include "absl/strings/str_format.h"
 #include "ortools/base/commandlineflags.h"
 #include "ortools/base/integral_types.h"
diff --git a/ortools/linear_solver/gurobi_proto_solver.cc b/ortools/linear_solver/gurobi_proto_solver.cc
new file mode 100644
index 0000000000..3030225e85
--- /dev/null
+++ b/ortools/linear_solver/gurobi_proto_solver.cc
@@ -0,0 +1,265 @@
+// Copyright 2010-2018 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.
+
+#if defined(USE_GUROBI)
+
+#include "ortools/linear_solver/gurobi_proto_solver.h"
+
+#include <limits>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "absl/strings/str_cat.h"
+#include "absl/strings/str_format.h"
+#include "absl/strings/str_split.h"
+#include "ortools/base/canonical_errors.h"
+#include "ortools/base/cleanup.h"
+#include "ortools/base/status.h"
+#include "ortools/base/status_macros.h"
+#include "ortools/linear_solver/linear_solver.pb.h"
+#include "ortools/linear_solver/model_validator.h"
+
+extern "C" {
+#include "gurobi_c.h"
+int __stdcall GRBisqp(GRBenv**, const char*, const char*, const char*, int,
+                      const char*);
+}
+
+namespace operations_research {
+
+namespace {
+constexpr int GRB_OK = 0;
+
+inline util::Status GurobiCodeToUtilStatus(int error_code,
+                                           const char* source_file,
+                                           int source_line,
+                                           const char* statement,
+                                           GRBenv* const env) {
+  if (error_code == GRB_OK) return util::OkStatus();
+  return util::InvalidArgumentError(absl::StrFormat(
+      "Gurobi error code %d (file '%s', line %d) on '%s': %s", error_code,
+      source_file, source_line, statement, GRBgeterrormsg(env)));
+}
+
+util::Status CreateEnvironment(GRBenv** gurobi) {
+  const char kGurobiEnvErrorMsg[] =
+      "Could not load Gurobi environment. Is gurobi correctly installed and "
+      "licensed on this machine?";
+
+  if (GRBloadenv(gurobi, nullptr) != 0 || *gurobi == nullptr) {
+    return util::FailedPreconditionError(
+        absl::StrFormat("%s %s", kGurobiEnvErrorMsg, GRBgeterrormsg(*gurobi)));
+  }
+  return util::OkStatus();
+}
+}  // namespace
+
+util::StatusOr<MPSolutionResponse> GurobiSolveProto(
+    const MPModelRequest& request) {
+  MPSolutionResponse response;
+  if (MPRequestIsEmptyOrInvalid(request, &response)) {
+    return response;
+  }
+
+  const MPModelProto& model = request.model();
+  if (request.has_solver_specific_parameters()) {
+    // TODO(user): Support solver-specific parameters without duplicating
+    // all of the write file / read file code in linear_solver.cc.
+    return util::UnimplementedError(
+        "Solver-specific parameters not supported.");
+  }
+  if (model.general_constraint_size() > 0) {
+    // TODO(user): Move indicator constraint logic from linear_solver.cc to
+    // this file.
+    return util::UnimplementedError("General constraints not supported.");
+  }
+  if (model.has_solution_hint()) {
+    // TODO(user): Support solution hints.
+    return util::UnimplementedError("Solution hint not supported.");
+  }
+
+  GRBenv* gurobi = nullptr;
+  GRBmodel* gurobi_model = nullptr;
+
+// `gurobi` references ther GRBenv variable defined above.
+#define RETURN_IF_GUROBI_ERROR(x) \
+  RETURN_IF_ERROR(GurobiCodeToUtilStatus(x, __FILE__, __LINE__, #x, gurobi));
+
+  auto delete_gurobi_objects = [&]() -> util::Status {
+    // Release all created pointers.
+    if (gurobi_model) RETURN_IF_GUROBI_ERROR(GRBfreemodel(gurobi_model));
+    if (gurobi) GRBfreeenv(gurobi);
+    return util::OkStatus();
+  };
+  auto gurobi_deleter = gtl::MakeCleanup([delete_gurobi_objects]() {
+    const util::Status deleter_status = delete_gurobi_objects();
+    LOG_IF(DFATAL, !deleter_status.ok()) << deleter_status;
+  });
+
+  RETURN_IF_ERROR(CreateEnvironment(&gurobi));
+  RETURN_IF_GUROBI_ERROR(GRBnewmodel(gurobi, &gurobi_model,
+                                     model.name().c_str(),
+                                     /*numvars=*/0,
+                                     /*obj=*/nullptr,
+                                     /*lb=*/nullptr,
+                                     /*ub=*/nullptr,
+                                     /*vtype=*/nullptr,
+                                     /*varnames=*/nullptr));
+
+  if (request.solver_time_limit_seconds() > 0) {
+    RETURN_IF_GUROBI_ERROR(GRBsetdblparam(GRBgetenv(gurobi_model),
+                                          GRB_DBL_PAR_TIMELIMIT,
+                                          request.solver_time_limit_seconds()));
+  }
+  RETURN_IF_GUROBI_ERROR(
+      GRBsetintparam(GRBgetenv(gurobi_model), GRB_INT_PAR_OUTPUTFLAG,
+                     request.enable_internal_solver_output()));
+
+  const int variable_size = model.variable_size();
+  bool has_integer_variables = false;
+  {
+    std::vector<double> obj_coeffs(variable_size, 0);
+    std::vector<double> lb(variable_size);
+    std::vector<double> ub(variable_size);
+    std::vector<char> ctype(variable_size);
+    std::vector<const char*> varnames(variable_size);
+    for (int v = 0; v < variable_size; ++v) {
+      const MPVariableProto& variable = model.variable(v);
+      obj_coeffs[v] = variable.objective_coefficient();
+      lb[v] = variable.lower_bound();
+      ub[v] = variable.upper_bound();
+      ctype[v] = variable.is_integer() ? GRB_INTEGER : GRB_CONTINUOUS;
+      if (variable.is_integer()) has_integer_variables = true;
+      if (!variable.name().empty()) varnames[v] = variable.name().c_str();
+    }
+
+    RETURN_IF_GUROBI_ERROR(
+        GRBaddvars(gurobi_model, variable_size, 0, nullptr, nullptr, nullptr,
+                   /*obj=*/obj_coeffs.data(),
+                   /*lb=*/lb.data(), /*ub=*/ub.data(), /*vtype=*/ctype.data(),
+                   /*varnames=*/const_cast<char**>(varnames.data())));
+  }
+
+  {
+    std::vector<int> ct_variables;
+    std::vector<double> ct_coefficients;
+    for (int c = 0; c < model.constraint_size(); ++c) {
+      const MPConstraintProto& constraint = model.constraint(c);
+      const int size = constraint.var_index_size();
+      ct_variables.resize(size, 0);
+      ct_coefficients.resize(size, 0);
+      for (int i = 0; i < size; ++i) {
+        ct_variables[i] = constraint.var_index(i);
+        ct_coefficients[i] = constraint.coefficient(i);
+      }
+      RETURN_IF_GUROBI_ERROR(GRBaddrangeconstr(
+          gurobi_model, /*numnz=*/size, /*cind=*/ct_variables.data(),
+          /*cval=*/ct_coefficients.data(), /*lower=*/constraint.lower_bound(),
+          /*upper=*/constraint.upper_bound(),
+          /*constrname=*/constraint.name().c_str()));
+    }
+  }
+
+  RETURN_IF_GUROBI_ERROR(GRBsetintattr(gurobi_model, GRB_INT_ATTR_MODELSENSE,
+                                       model.maximize() ? -1 : 1));
+  RETURN_IF_GUROBI_ERROR(GRBsetdblattr(gurobi_model, GRB_DBL_ATTR_OBJCON,
+                                       model.objective_offset()));
+  if (model.has_quadratic_objective()) {
+    MPQuadraticObjective qobj = model.quadratic_objective();
+    if (qobj.coefficient_size() > 0) {
+      RETURN_IF_GUROBI_ERROR(
+          GRBaddqpterms(gurobi_model, /*numqnz=*/qobj.coefficient_size(),
+                        /*qrow=*/qobj.mutable_qvar1_index()->mutable_data(),
+                        /*qcol=*/qobj.mutable_qvar2_index()->mutable_data(),
+                        /*qval=*/qobj.mutable_coefficient()->mutable_data()));
+    }
+  }
+
+  RETURN_IF_GUROBI_ERROR(GRBupdatemodel(gurobi_model));
+  RETURN_IF_GUROBI_ERROR(GRBoptimize(gurobi_model));
+
+  int optimization_status = 0;
+  RETURN_IF_GUROBI_ERROR(
+      GRBgetintattr(gurobi_model, GRB_INT_ATTR_STATUS, &optimization_status));
+  int solution_count = 0;
+  RETURN_IF_GUROBI_ERROR(
+      GRBgetintattr(gurobi_model, GRB_INT_ATTR_SOLCOUNT, &solution_count));
+  switch (optimization_status) {
+    case GRB_OPTIMAL:
+      response.set_status(MPSOLVER_OPTIMAL);
+      break;
+    case GRB_INF_OR_UNBD:
+      DLOG(INFO) << "Gurobi solve returned GRB_INF_OR_UNBD, which we treat as "
+                    "INFEASIBLE even though it may mean UNBOUNDED.";
+      response.set_status_str(
+          "The model may actually be unbounded: Gurobi returned "
+          "GRB_INF_OR_UNBD");
+      ABSL_FALLTHROUGH_INTENDED;
+    case GRB_INFEASIBLE:
+      response.set_status(MPSOLVER_INFEASIBLE);
+      break;
+    case GRB_UNBOUNDED:
+      response.set_status(MPSOLVER_UNBOUNDED);
+      break;
+    default: {
+      if (solution_count > 0) {
+        response.set_status(MPSOLVER_FEASIBLE);
+      } else {
+        response.set_status(MPSOLVER_NOT_SOLVED);
+        response.set_status_str(
+            absl::StrFormat("Gurobi status code %d", optimization_status));
+      }
+      break;
+    }
+  }
+
+  if (solution_count > 0 && (response.status() == MPSOLVER_FEASIBLE ||
+                             response.status() == MPSOLVER_OPTIMAL)) {
+    double objective_value = 0;
+    RETURN_IF_GUROBI_ERROR(
+        GRBgetdblattr(gurobi_model, GRB_DBL_ATTR_OBJVAL, &objective_value));
+    response.set_objective_value(objective_value);
+    if (has_integer_variables) {
+      double best_objective_bound = 0;
+      const int error = GRBgetdblattr(gurobi_model, GRB_DBL_ATTR_OBJBOUND,
+                                      &best_objective_bound);
+      if (response.status() == MPSOLVER_OPTIMAL &&
+          error == GRB_ERROR_DATA_NOT_AVAILABLE) {
+        // If the presolve deletes all variables, there's no best bound.
+        response.set_best_objective_bound(objective_value);
+      } else {
+        RETURN_IF_GUROBI_ERROR(error);
+        response.set_best_objective_bound(best_objective_bound);
+      }
+    }
+
+    response.mutable_variable_value()->Resize(variable_size, 0);
+    RETURN_IF_GUROBI_ERROR(
+        GRBgetdblattrarray(gurobi_model, GRB_DBL_ATTR_X, 0, variable_size,
+                           response.mutable_variable_value()->mutable_data()));
+    if (!has_integer_variables) {
+      response.mutable_dual_value()->Resize(model.constraint_size(), 0);
+      RETURN_IF_GUROBI_ERROR(GRBgetdblattrarray(
+          gurobi_model, GRB_DBL_ATTR_PI, 0, model.constraint_size(),
+          response.mutable_dual_value()->mutable_data()));
+    }
+  }
+#undef RETURN_IF_GUROBI_ERROR
+
+  return response;
+}
+
+}  // namespace operations_research
+
+#endif  //  #if defined(USE_GUROBI)
diff --git a/ortools/linear_solver/gurobi_proto_solver.h b/ortools/linear_solver/gurobi_proto_solver.h
new file mode 100644
index 0000000000..f27d09e851
--- /dev/null
+++ b/ortools/linear_solver/gurobi_proto_solver.h
@@ -0,0 +1,26 @@
+// Copyright 2010-2018 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_LINEAR_SOLVER_GUROBI_PROTO_SOLVER_H_
+#define OR_TOOLS_LINEAR_SOLVER_GUROBI_PROTO_SOLVER_H_
+
+#include "ortools/base/statusor.h"
+#include "ortools/linear_solver/linear_solver.pb.h"
+
+namespace operations_research {
+
+util::StatusOr<MPSolutionResponse> GurobiSolveProto(
+    const MPModelRequest& request);
+
+}  // namespace operations_research
+#endif  // OR_TOOLS_LINEAR_SOLVER_GUROBI_PROTO_SOLVER_H_
diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc
index 522708b4c6..6c6476d307 100644
--- a/ortools/linear_solver/linear_solver.cc
+++ b/ortools/linear_solver/linear_solver.cc
@@ -598,6 +598,14 @@ MPSolverResponseStatus MPSolver::LoadModelFromProtoInternal(
     }
   }
 
+  if (input_model.has_quadratic_objective()) {
+    *error_message =
+        "Optimizing a quadratic objective is only supported through direct "
+        "proto solves. Please use MPSolver::SolveWithProto, or the solver's "
+        "direct proto solve function.";
+    return MPSOLVER_MODEL_INVALID;
+  }
+
   MPObjective* const objective = MutableObjective();
   // Passing empty names makes the MPSolver generate unique names.
   const std::string empty;
@@ -631,41 +639,49 @@ MPSolverResponseStatus MPSolver::LoadModelFromProtoInternal(
 
   for (const MPGeneralConstraintProto& general_constraint :
        input_model.general_constraint()) {
-    if (general_constraint.has_indicator_constraint()) {
-      const auto& proto =
-          general_constraint.indicator_constraint().constraint();
-      if (proto.lower_bound() == -infinity() &&
-          proto.upper_bound() == infinity()) {
-        continue;
+    switch (general_constraint.general_constraint_case()) {
+      case MPGeneralConstraintProto::kIndicatorConstraint: {
+        const auto& proto =
+            general_constraint.indicator_constraint().constraint();
+        if (proto.lower_bound() == -infinity() &&
+            proto.upper_bound() == infinity()) {
+          continue;
+        }
+
+        const int constraint_index = NumConstraints();
+        MPConstraint* const constraint = new MPConstraint(
+            constraint_index, proto.lower_bound(), proto.upper_bound(),
+            clear_names ? "" : proto.name(), interface_.get());
+        if (constraint_name_to_index_) {
+          gtl::InsertOrDie(&*constraint_name_to_index_, constraint->name(),
+                           constraint_index);
+        }
+        constraints_.push_back(constraint);
+        constraint_is_extracted_.push_back(false);
+
+        constraint->set_is_lazy(proto.is_lazy());
+        for (int j = 0; j < proto.var_index_size(); ++j) {
+          constraint->SetCoefficient(variables_[proto.var_index(j)],
+                                     proto.coefficient(j));
+        }
+
+        MPVariable* const variable =
+            variables_[general_constraint.indicator_constraint().var_index()];
+        constraint->indicator_variable_ = variable;
+        constraint->indicator_value_ =
+            general_constraint.indicator_constraint().var_value();
+
+        if (!interface_->AddIndicatorConstraint(constraint)) {
+          *error_message = "Solver doesn't support indicator constraints";
+          return MPSOLVER_MODEL_INVALID;
+        }
+        break;
       }
-
-      const int constraint_index = NumConstraints();
-      MPConstraint* const constraint = new MPConstraint(
-          constraint_index, proto.lower_bound(), proto.upper_bound(),
-          clear_names ? "" : proto.name(), interface_.get());
-      if (constraint_name_to_index_) {
-        gtl::InsertOrDie(&*constraint_name_to_index_, constraint->name(),
-                         constraint_index);
-      }
-      constraints_.push_back(constraint);
-      constraint_is_extracted_.push_back(false);
-
-      constraint->set_is_lazy(proto.is_lazy());
-      for (int j = 0; j < proto.var_index_size(); ++j) {
-        constraint->SetCoefficient(variables_[proto.var_index(j)],
-                                   proto.coefficient(j));
-      }
-
-      MPVariable* const variable =
-          variables_[general_constraint.indicator_constraint().var_index()];
-      constraint->indicator_variable_ = variable;
-      constraint->indicator_value_ =
-          general_constraint.indicator_constraint().var_value();
-
-      if (!interface_->AddIndicatorConstraint(constraint)) {
-        *error_message = "Solver doesn't support indicator constraints.";
+      default:
+        *error_message =
+            absl::StrCat("Solver doesn't support general constraints of type ",
+                         general_constraint.general_constraint_case());
         return MPSOLVER_MODEL_INVALID;
-      }
     }
   }
 
@@ -744,14 +760,21 @@ void MPSolver::SolveWithProto(const MPModelRequest& model_request,
   if (model_request.enable_internal_solver_output()) {
     solver.EnableOutput();
   }
+
+  auto optional_response = solver.interface_->DirectlySolveProto(model_request);
+  if (optional_response) {
+    *response = std::move(optional_response).value();
+    return;
+  }
+
   std::string error_message;
   response->set_status(solver.LoadModelFromProto(model, &error_message));
   if (response->status() != MPSOLVER_MODEL_IS_VALID) {
-    LOG(WARNING) << "Loading model from protocol buffer failed, load status = "
-                 << ProtoEnumToString<MPSolverResponseStatus>(
-                        response->status())
-                 << " (" << response->status() << "); Error: " << error_message;
-
+    response->set_status_str(error_message);
+    LOG_IF(WARNING, model_request.enable_internal_solver_output())
+        << "Loading model from protocol buffer failed, load status = "
+        << ProtoEnumToString<MPSolverResponseStatus>(response->status()) << " ("
+        << response->status() << "); Error: " << error_message;
     return;
   }
   if (model_request.has_solver_time_limit_seconds()) {
diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h
index 2cd65b523c..b0f14ef522 100644
--- a/ortools/linear_solver/linear_solver.h
+++ b/ortools/linear_solver/linear_solver.h
@@ -1482,6 +1482,14 @@ class MPSolverInterface {
   // solution is optimal.
   virtual MPSolver::ResultStatus Solve(const MPSolverParameters& param) = 0;
 
+  // Directly solves a MPModelRequest, bypassing the MPSolver data structures
+  // entirely. Returns {} (eg. absl::nullopt) if the feature is not supported by
+  // the underlying solver.
+  virtual absl::optional<MPSolutionResponse> DirectlySolveProto(
+      const MPModelRequest& request) {
+    return absl::nullopt;
+  }
+
   // Writes the model using the solver internal write function.  Currently only
   // available for GurobiInterface.
   virtual void Write(const std::string& filename);
diff --git a/ortools/linear_solver/linear_solver.proto b/ortools/linear_solver/linear_solver.proto
index 4915a7c86a..5c31f106c7 100644
--- a/ortools/linear_solver/linear_solver.proto
+++ b/ortools/linear_solver/linear_solver.proto
@@ -161,6 +161,24 @@ message MPSosConstraint {
   repeated double weight = 3;
 }
 
+// Quadratic part of a model's objective. Added with other objectives (such as
+// linear), this creates the model's objective function to be optimized.
+// Note: the linear part of the objective currently needs to be specified in the
+// MPVariableProto.objective_coefficient fields. If you'd rather have a
+// dedicated linear array here, talk to or-core-team@
+message MPQuadraticObjective {
+  // Sparse representation of quadratic terms in the objective function, where
+  // term i is qvar1_index[i] * qvar2_index[i] * coefficient[i].
+  // `qvar1_index` and `qvar2_index` are variable indices w.r.t the "variable"
+  // field in MPModelProto.
+  // `qvar1_index`, `qvar2_index` and `coefficients` must have the same size.
+  // If the same unordered pair (qvar1_index, qvar2_index) appears several
+  // times, the sum of all of the associated coefficients will be applied.
+  repeated int32 qvar1_index = 1;
+  repeated int32 qvar2_index = 2;
+  repeated double coefficient = 3;  // Must be finite.
+}
+
 // This message encodes a partial (or full) assignment of the variables of a
 // MPModelProto problem. The indices in var_index should be unique and valid
 // variable indices of the associated problem.
@@ -171,12 +189,6 @@ message PartialVariableAssignment {
 
 // MPModelProto contains all the information for a Linear Programming model.
 message MPModelProto {
-  // True if the problem is a maximization problem. Minimize by default.
-  optional bool maximize = 1 [default = false];
-
-  // Offset for the objective function. Must be finite.
-  optional double objective_offset = 2 [default = 0.0];
-
   // All the variables appearing in the model.
   repeated MPVariableProto variable = 3;
 
@@ -187,6 +199,16 @@ message MPModelProto {
   // solvers support all types of general constraints.
   repeated MPGeneralConstraintProto general_constraint = 7;
 
+  // True if the problem is a maximization problem. Minimize by default.
+  optional bool maximize = 1 [default = false];
+
+  // Offset for the objective function. Must be finite.
+  optional double objective_offset = 2 [default = 0.0];
+
+  // Optionally, a quadratic objective.
+  // As of 2019/06, only SCIP and Gurobi support quadratic objectives.
+  optional MPQuadraticObjective quadratic_objective = 8;
+
   // Name of the model.
   optional string name = 5 [default = ""];
 
@@ -270,6 +292,32 @@ message MPSolverCommonParameters {
   optional OptionalBoolean scaling = 7 [default = BOOL_UNSPECIFIED];
 }
 
+// Encodes a full MPModelProto by way of referencing to a "baseline"
+// MPModelProto stored in a file, and a "delta" to apply to this model.
+message MPModelDeltaProto {
+  optional /*required*/ string baseline_model_file_path = 1;
+
+  // The variable protos listed here will override (via MergeFrom()) the ones
+  // in the baseline model: you only need to specify the fields that change.
+  // To add a new variable, add it with a new variable index (variable indices
+  // still need to span a dense integer interval).
+  // You can't "delete" a variable but you can "neutralize" it by fixing its
+  // value, setting its objective coefficient to zero, and by nullifying all
+  // the terms involving it in the constraints.
+  map<int32, MPVariableProto> variable_overrides = 2;
+
+  // Constraints can be changed (or added) in the same way as variables, see
+  // above. It's mostly like applying MergeFrom(), except that:
+  // - the "var_index" and "coefficient" fields will be overridden like a map:
+  //   if a key pre-exists, we overwrite its value, otherwise we add it.
+  // - if you set the lower bound to -inf and the upper bound to +inf, thus
+  //   effectively neutralizing the constraint, the solver will implicitly
+  //   remove all of the constraint's terms.
+  map<int32, MPConstraintProto> constraint_overrides = 3;
+
+  // NOTE(user): We may also add more deltas, eg. objective offset.
+}
+
 message MPModelRequest {
   // The model to be optimized by the server.
   optional MPModelProto model = 1;
@@ -420,6 +468,12 @@ message MPSolutionResponse {
   optional /*required*/ MPSolverResponseStatus status = 1
       [default = MPSOLVER_UNKNOWN_STATUS];
 
+  // Human-readable string giving more details about the status. For example,
+  // when the status is MPSOLVER_INVALID_MODE, this can hold a description of
+  // why the model is invalid.
+  // This isn't always filled: don't depend on its value or even its presence.
+  optional string status_str = 7;
+
   // Objective value corresponding to the "variable_value" below, taking into
   // account the source "objective_offset" and "objective_coefficient".
   // This is set iff 'status' is OPTIMAL or FEASIBLE.
diff --git a/ortools/linear_solver/model_validator.cc b/ortools/linear_solver/model_validator.cc
index 5fdb21a96b..9b0bea903a 100644
--- a/ortools/linear_solver/model_validator.cc
+++ b/ortools/linear_solver/model_validator.cc
@@ -17,7 +17,9 @@
 #include <cmath>
 #include <limits>
 
+#include "absl/container/flat_hash_set.h"
 #include "absl/strings/str_cat.h"
+#include "absl/strings/str_format.h"
 #include "ortools/base/accurate_sum.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
 #include "ortools/port/proto_utils.h"
@@ -192,6 +194,32 @@ std::string FindErrorInMPSosConstraint(const MPModelProto& model,
   return "";
 }
 
+std::string FindErrorInQuadraticObjective(const MPQuadraticObjective& qobj,
+                                          int num_vars) {
+  if (qobj.qvar1_index_size() != qobj.qvar2_index_size() ||
+      qobj.qvar1_index_size() != qobj.coefficient_size()) {
+    return "indices and coefficients must have the same size";
+  }
+
+  for (int i = 0; i < qobj.qvar1_index_size(); ++i) {
+    if (qobj.qvar1_index(i) >= num_vars || qobj.qvar1_index(i) < 0) {
+      return absl::StrCat("qvar1_index(", i, ")=", qobj.qvar1_index(i),
+                          " is invalid.", " It must be in [0, ", num_vars, ")");
+    }
+
+    if (qobj.qvar2_index(i) >= num_vars || qobj.qvar2_index(i) < 0) {
+      return absl::StrCat("qvar2_index(", i, ")=", qobj.qvar2_index(i),
+                          " is invalid.", " It must be in [0, ", num_vars, ")");
+    }
+
+    if (!std::isfinite(qobj.coefficient(i))) {
+      return absl::StrCat("coefficient(", i, ")=", (qobj.coefficient(i)),
+                          " is invalid");
+    }
+  }
+  return "";
+}
+
 std::string FindErrorInSolutionHint(
     const PartialVariableAssignment& solution_hint, int num_vars) {
   if (solution_hint.var_index_size() != solution_hint.var_value_size()) {
@@ -220,11 +248,9 @@ std::string FindErrorInSolutionHint(
 }  // namespace
 
 std::string FindErrorInMPModelProto(const MPModelProto& model) {
-  // TODO(user): enhance the error reporting: report several errors instead of
-  // stopping at the first one.
-
-  // TODO(user): clarify explicitly, at least in a comment, whether we
-  // accept models without variables and/or constraints.
+  // NOTE(user): Empty models are considered fine by this function, although
+  // it is not clear whether MPSolver::Solve() will always respond in the same
+  // way, depending on the solvers.
 
   if (!std::isfinite(model.objective_offset())) {
     return absl::StrCat("Invalid objective_offset: ",
@@ -280,6 +306,13 @@ std::string FindErrorInMPModelProto(const MPModelProto& model) {
     }
   }
 
+  // Validate objectives.
+  if (model.has_quadratic_objective()) {
+    error =
+        FindErrorInQuadraticObjective(model.quadratic_objective(), num_vars);
+    if (!error.empty()) return absl::StrCat("In quadratic_objective: ", error);
+  }
+
   // Validate the solution hint.
   error = FindErrorInSolutionHint(model.solution_hint(), num_vars);
   if (!error.empty()) {
@@ -289,6 +322,39 @@ std::string FindErrorInMPModelProto(const MPModelProto& model) {
   return std::string();
 }
 
+bool MPRequestIsEmptyOrInvalid(const MPModelRequest& request,
+                               MPSolutionResponse* response) {
+  CHECK(response != nullptr);
+
+  if (!request.has_model()) {
+    response->set_status(MPSOLVER_OPTIMAL);
+    response->set_status_str("Requests without model are considered OPTIMAL");
+    return true;
+  }
+  const MPModelProto& model = request.model();
+  if (model.variable_size() == 0 && model.constraint_size() == 0) {
+    response->set_status(MPSOLVER_OPTIMAL);
+    response->set_objective_value(request.model().objective_offset());
+    response->set_best_objective_bound(request.model().objective_offset());
+    response->set_status_str(
+        "Requests without variables and constraints are considered OPTIMAL");
+    return true;
+  }
+
+  const std::string error = FindErrorInMPModelProto(model);
+  if (!error.empty()) {
+    if (request.enable_internal_solver_output()) {
+      LOG(ERROR) << absl::StrCat("Invalid model: ", error);
+    }
+    response->set_status(error.find("Infeasible") == std::string::npos
+                             ? MPSOLVER_MODEL_INVALID
+                             : MPSOLVER_INFEASIBLE);
+    response->set_status_str(error);
+    return true;
+  }
+  return false;
+}
+
 // TODO(user): Add a general FindFeasibilityErrorInSolution() and factor out the
 // common code.
 std::string FindFeasibilityErrorInSolutionHint(const MPModelProto& model,
@@ -350,4 +416,104 @@ std::string FindFeasibilityErrorInSolutionHint(const MPModelProto& model,
   return "";
 }
 
+std::string FindErrorInMPModelDeltaProto(const MPModelDeltaProto& delta,
+                                         const MPModelProto& model) {
+  int num_vars = model.variable_size();
+  // Validate delta variables.
+  std::string error;
+  absl::flat_hash_set<int> new_var_indices;
+  int max_var_index = num_vars - 1;
+  MPVariableProto tmp_var_proto;
+  for (const auto& pair : delta.variable_overrides()) {
+    const int var_index = pair.first;
+    const MPVariableProto& var_override_proto = pair.second;
+    if (var_index < 0) {
+      error = "Invalid key";
+    } else if (var_index >= num_vars) {
+      max_var_index = std::max(max_var_index, var_index);
+      new_var_indices.insert(var_index);
+      error = FindErrorInMPVariable(var_override_proto);
+    } else {
+      tmp_var_proto = model.variable(var_index);
+      // NOTE(user): It is OK for the override proto to be empty, i.e. be a
+      // non-override.
+      tmp_var_proto.MergeFrom(var_override_proto);
+      error = FindErrorInMPVariable(tmp_var_proto);
+    }
+    if (!error.empty()) {
+      return absl::StrFormat(
+          "variable_overrides with key (eg. var index) = %d: %s", var_index,
+          error);
+    }
+  }
+  if (max_var_index != num_vars + new_var_indices.size() - 1) {
+    return absl::StrFormat(
+        "The added and existing variable indices do not form a dense integer "
+        "interval: oldmax=%d, max=%d, num added=%d",
+        num_vars - 1, max_var_index, new_var_indices.size());
+  }
+  // Now we "officially" add the new variables to "num_vars".
+  num_vars += new_var_indices.size();
+
+  // Validate delta constraints. We can avoid going over the full
+  // var_index/coefficient of the original constraint, since the overrides are
+  // self-sufficient (i.e. the override var_index/coefficients are valid iff
+  // they would be valid in a standalone, new constraint). So we use a partial
+  // proto merger to avoid those in the baseline constraint.
+  std::vector<bool> variable_appears(num_vars, false);
+  MPConstraintProto tmp_constraint_proto;
+  const int num_constraints = model.constraint_size();
+  absl::flat_hash_set<int> new_ct_indices;
+  int max_ct_index = num_constraints - 1;
+  for (const auto& pair : delta.constraint_overrides()) {
+    const int ct_index = pair.first;
+    const MPConstraintProto& constraint_override_proto = pair.second;
+    if (ct_index < 0) {
+      error = "Invalid constraint index";
+    } else if (ct_index >= num_constraints) {
+      max_ct_index = std::max(max_ct_index, ct_index);
+      new_ct_indices.insert(ct_index);
+      error =
+          FindErrorInMPConstraint(constraint_override_proto, &variable_appears);
+    } else {
+      // NOTE(user): We don't need to do the merging of var_index/coefficient:
+      // that part of the merged constraint will be valid iff the override is
+      // valid as a standalone var_index/coefficient map.
+      // So we simply validate a reduced version of the actual "merged"
+      // constraint, by removing the var_index/coefficient of the baseline.
+      // Benefit: the complexity is O(|constraint override|) even if the
+      // baseline constraint was huge.
+      tmp_constraint_proto.Clear();
+      MergeMPConstraintProtoExceptTerms(model.constraint(ct_index),
+                                        &tmp_constraint_proto);
+      tmp_constraint_proto.MergeFrom(constraint_override_proto);
+      error = FindErrorInMPConstraint(tmp_constraint_proto, &variable_appears);
+    }
+    if (!error.empty()) {
+      return absl::StrFormat(
+          "constraint_overrides with key (eg. constraint index) = %d: %s",
+          ct_index, error);
+    }
+  }
+  if (max_ct_index != num_constraints + new_ct_indices.size() - 1) {
+    return absl::StrFormat(
+        "The added and existing constraint indices do not form a dense integer "
+        "interval: oldmax=%d, max=%d, num added=%d",
+        num_constraints - 1, max_ct_index, new_ct_indices.size());
+  }
+
+  return "";
+}
+
+void MergeMPConstraintProtoExceptTerms(const MPConstraintProto& from,
+                                       MPConstraintProto* to) {
+#define COPY_FIELD_IF_PRESENT(field) \
+  if (from.has_##field()) to->set_##field(from.field())
+  COPY_FIELD_IF_PRESENT(lower_bound);
+  COPY_FIELD_IF_PRESENT(upper_bound);
+  COPY_FIELD_IF_PRESENT(name);
+  COPY_FIELD_IF_PRESENT(is_lazy);
+#undef COPY_FIELD_IF_PRESENT
+}
+
 }  // namespace operations_research
diff --git a/ortools/linear_solver/model_validator.h b/ortools/linear_solver/model_validator.h
index 8e4bf6b112..c8bc9f0ac5 100644
--- a/ortools/linear_solver/model_validator.h
+++ b/ortools/linear_solver/model_validator.h
@@ -30,6 +30,23 @@ namespace operations_research {
  */
 std::string FindErrorInMPModelProto(const MPModelProto& model);
 
+/**
+ *  Like FindErrorInMPModelProto, but for a MPModelDeltaProto applied to a given
+ * baseline model (assumed valid, eg. FindErrorInMPModelProto(model)="").
+ * Works in O(|model_delta|) + O(num_vars in model), but the latter term has a
+ * very small constant factor.
+ */
+std::string FindErrorInMPModelDeltaProto(const MPModelDeltaProto& delta,
+                                         const MPModelProto& model);
+
+/**
+ *  Updates `response` and returns true if errors, infeasibilities, or trivial
+ * optimals were found. Returns false if the model is valid and non-trivially
+ * solvable.
+ */
+bool MPRequestIsEmptyOrInvalid(const MPModelRequest& request,
+                               MPSolutionResponse* response);
+
 /**
  * Returns an empty std::string if the solution hint given in the model is a
  * feasible solution. Otherwise, returns a description of the first reason for
@@ -43,6 +60,15 @@ std::string FindErrorInMPModelProto(const MPModelProto& model);
 std::string FindFeasibilityErrorInSolutionHint(const MPModelProto& model,
                                                double tolerance);
 
+// PUBLIC ONLY FOR TESTING.
+// Partially merges a MPConstraintProto onto another, skipping only the
+// repeated fields "var_index" and "coefficients". This is used within
+// FindErrorInMPModelDeltaProto.
+// See the unit test MergeMPConstraintProtoExceptTermsTest that explains why we
+// need this.
+void MergeMPConstraintProtoExceptTerms(const MPConstraintProto& from,
+                                       MPConstraintProto* to);
+
 }  // namespace operations_research
 
 #endif  // OR_TOOLS_LINEAR_SOLVER_MODEL_VALIDATOR_H_
diff --git a/ortools/linear_solver/samples/integer_programming_example.py b/ortools/linear_solver/samples/integer_programming_example.py
index bc8d4f68c8..18354bf20f 100644
--- a/ortools/linear_solver/samples/integer_programming_example.py
+++ b/ortools/linear_solver/samples/integer_programming_example.py
@@ -15,7 +15,6 @@
 from __future__ import print_function
 # [START import]
 from ortools.linear_solver import pywraplp
-
 # [END import]
 
 
diff --git a/ortools/linear_solver/samples/simple_lp_program.py b/ortools/linear_solver/samples/simple_lp_program.py
index 68912fe76d..0f4b4561f1 100644
--- a/ortools/linear_solver/samples/simple_lp_program.py
+++ b/ortools/linear_solver/samples/simple_lp_program.py
@@ -15,7 +15,6 @@
 # [START import]
 from __future__ import print_function
 from ortools.linear_solver import pywraplp
-
 # [END import]
 
 
diff --git a/ortools/linear_solver/samples/simple_mip_program.py b/ortools/linear_solver/samples/simple_mip_program.py
index d6598ace74..167b5584db 100644
--- a/ortools/linear_solver/samples/simple_mip_program.py
+++ b/ortools/linear_solver/samples/simple_mip_program.py
@@ -15,7 +15,6 @@
 # [START import]
 from __future__ import print_function
 from ortools.linear_solver import pywraplp
-
 # [END import]
 
 
diff --git a/ortools/linear_solver/scip_helper_macros.h b/ortools/linear_solver/scip_helper_macros.h
new file mode 100644
index 0000000000..00fadb2dfd
--- /dev/null
+++ b/ortools/linear_solver/scip_helper_macros.h
@@ -0,0 +1,46 @@
+// Copyright 2010-2018 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_LINEAR_SOLVER_SCIP_HELPER_MACROS_H_
+#define OR_TOOLS_LINEAR_SOLVER_SCIP_HELPER_MACROS_H_
+
+#include "absl/strings/str_format.h"
+#include "ortools/base/canonical_errors.h"
+#include "ortools/base/status.h"
+
+namespace operations_research {
+namespace internal {
+// Our own version of SCIP_CALL to do error management.
+// NOTE(user): There are so many SCIP error codes, in so many different
+// situations. We don't try to match them perfectly to google3 error codes.
+// Instead, we use the most likely/generic code "invalid argument" and surface
+// the internal SCIP error code to the user.
+inline util::Status ScipCodeToUtilStatus(/*SCIP_Retcode*/ int retcode,
+                                         const char* source_file,
+                                         int source_line,
+                                         const char* scip_statement) {
+  if (retcode == /*SCIP_OKAY*/ 1) return util::OkStatus();
+  return util::InvalidArgumentError(
+      absl::StrFormat("SCIP error code %d (file '%s', line %d) on '%s'",
+                      retcode, source_file, source_line, scip_statement));
+}
+}  // namespace internal
+
+#define SCIP_TO_STATUS(x) \
+  internal::ScipCodeToUtilStatus(x, __FILE__, __LINE__, #x)
+
+#define RETURN_IF_SCIP_ERROR(x) RETURN_IF_ERROR(SCIP_TO_STATUS(x));
+
+}  // namespace operations_research
+
+#endif  // OR_TOOLS_LINEAR_SOLVER_SCIP_HELPER_MACROS_H_
diff --git a/ortools/linear_solver/scip_interface.cc b/ortools/linear_solver/scip_interface.cc
index 9b21500d6e..c3e4646fa7 100644
--- a/ortools/linear_solver/scip_interface.cc
+++ b/ortools/linear_solver/scip_interface.cc
@@ -21,6 +21,7 @@
 #include <string>
 #include <vector>
 
+#include "absl/types/optional.h"
 #include "ortools/base/canonical_errors.h"
 #include "ortools/base/commandlineflags.h"
 #include "ortools/base/hash.h"
@@ -30,6 +31,10 @@
 #include "ortools/base/status_macros.h"
 #include "ortools/base/timer.h"
 #include "ortools/linear_solver/linear_solver.h"
+#include "ortools/linear_solver/linear_solver.pb.h"
+#include "ortools/linear_solver/scip_helper_macros.h"
+#include "ortools/linear_solver/scip_proto_solver.h"
+#include "scip/cons_indicator.h"
 #include "scip/scip.h"
 #include "scip/scipdefplugins.h"
 
@@ -46,6 +51,8 @@ class SCIPInterface : public MPSolverInterface {
 
   void SetOptimizationDirection(bool maximize) override;
   MPSolver::ResultStatus Solve(const MPSolverParameters& param) override;
+  absl::optional<MPSolutionResponse> DirectlySolveProto(
+      const MPModelRequest& request) override;
   void Reset() override;
 
   void SetVariableBounds(int var_index, double lb, double ub) override;
@@ -136,28 +143,11 @@ class SCIPInterface : public MPSolverInterface {
   util::Status status_;
 
   SCIP* scip_;
-  SCIP_VAR* objective_offset_variable_;
   std::vector<SCIP_VAR*> scip_variables_;
   std::vector<SCIP_CONS*> scip_constraints_;
   bool branching_priority_reset_ = false;
 };
 
-// Our own version of SCIP_CALL to do error management.
-// NOTE(user): There are so many SCIP error codes, in so many different
-// situations.. We don't try to match them perfectly to google3 error codes.
-// Instead, we use the most likely/generic code "invalid argument" and surface
-// the internal SCIP error code to the user.
-#define TO_STATUS(x) ScipReturnCodeToUtilStatus(x, __FILE__, __LINE__, #x)
-util::Status ScipReturnCodeToUtilStatus(SCIP_Retcode retcode,
-                                        const char* source_file,
-                                        int source_line,
-                                        const char* scip_statement) {
-  if (retcode == SCIP_OKAY) return util::OkStatus();
-  return util::InvalidArgumentError(
-      absl::StrFormat("SCIP error code %d (file '%s', line %d) on '%s'",
-                      retcode, source_file, source_line, scip_statement));
-}
-
 SCIPInterface::SCIPInterface(MPSolver* solver)
     : MPSolverInterface(solver), scip_(nullptr) {
   status_ = CreateSCIP();
@@ -172,34 +162,23 @@ void SCIPInterface::Reset() {
 }
 
 util::Status SCIPInterface::CreateSCIP() {
-  RETURN_IF_ERROR(TO_STATUS(SCIPcreate(&scip_)));
-  RETURN_IF_ERROR(TO_STATUS(SCIPincludeDefaultPlugins(scip_)));
+  RETURN_IF_SCIP_ERROR(SCIPcreate(&scip_));
+  RETURN_IF_SCIP_ERROR(SCIPincludeDefaultPlugins(scip_));
   // Set the emphasis to enum SCIP_PARAMEMPHASIS_FEASIBILITY. Do not print
   // the new parameter (quiet = true).
   if (FLAGS_scip_feasibility_emphasis) {
-    RETURN_IF_ERROR(
-        TO_STATUS(SCIPsetEmphasis(scip_, SCIP_PARAMEMPHASIS_FEASIBILITY,
-                                  /*quiet=*/true)));
+    RETURN_IF_SCIP_ERROR(SCIPsetEmphasis(scip_, SCIP_PARAMEMPHASIS_FEASIBILITY,
+                                         /*quiet=*/true));
   }
   // Default clock type. We use wall clock time because getting CPU user seconds
   // involves calling times() which is very expensive.
-  RETURN_IF_ERROR(TO_STATUS(
-      SCIPsetIntParam(scip_, "timing/clocktype", SCIP_CLOCKTYPE_WALL)));
-  RETURN_IF_ERROR(
-      TO_STATUS(SCIPcreateProb(scip_, solver_->name_.c_str(), nullptr, nullptr,
-                               nullptr, nullptr, nullptr, nullptr, nullptr)));
-  RETURN_IF_ERROR(TO_STATUS(SCIPsetObjsense(
-      scip_, maximize_ ? SCIP_OBJSENSE_MAXIMIZE : SCIP_OBJSENSE_MINIMIZE)));
-  // SCIPaddObjoffset cannot be used at the problem building stage. So we handle
-  // the objective offset by creating a dummy variable.
-  objective_offset_variable_ = nullptr;
-  // The true objective coefficient will be set in ExtractObjective.
-  double dummy_obj_coef = 0.0;
-  RETURN_IF_ERROR(TO_STATUS(
-      SCIPcreateVar(scip_, &objective_offset_variable_, "dummy", 1.0, 1.0,
-                    dummy_obj_coef, SCIP_VARTYPE_CONTINUOUS, true, false,
-                    nullptr, nullptr, nullptr, nullptr, nullptr)));
-  RETURN_IF_ERROR(TO_STATUS(SCIPaddVar(scip_, objective_offset_variable_)));
+  RETURN_IF_SCIP_ERROR(
+      SCIPsetIntParam(scip_, "timing/clocktype", SCIP_CLOCKTYPE_WALL));
+  RETURN_IF_SCIP_ERROR(SCIPcreateProb(scip_, solver_->name_.c_str(), nullptr,
+                                      nullptr, nullptr, nullptr, nullptr,
+                                      nullptr, nullptr));
+  RETURN_IF_SCIP_ERROR(SCIPsetObjsense(
+      scip_, maximize_ ? SCIP_OBJSENSE_MAXIMIZE : SCIP_OBJSENSE_MINIMIZE));
   return util::OkStatus();
 }
 
@@ -209,9 +188,6 @@ void SCIPInterface::DeleteSCIP() {
   // errors. The current code isn't perfect, since some CHECKs() remain, but
   // hopefully they'll never be triggered in practice.
   CHECK(scip_ != nullptr);
-  if (objective_offset_variable_ != nullptr) {
-    CHECK_EQ(SCIPreleaseVar(scip_, &objective_offset_variable_), SCIP_OKAY);
-  }
   for (int i = 0; i < scip_variables_.size(); ++i) {
     CHECK_EQ(SCIPreleaseVar(scip_, &scip_variables_[i]), SCIP_OKAY);
   }
@@ -232,18 +208,18 @@ void SCIPInterface::DeleteSCIP() {
     }                                                                    \
   } while (false)
 
-#define RETURN_IF_SCIP_ERROR(x) \
-  do {                          \
-    status_ = TO_STATUS(x);     \
-    if (!status_.ok()) return;  \
+#define RETURN_AND_STORE_IF_SCIP_ERROR(x) \
+  do {                                    \
+    status_ = SCIP_TO_STATUS(x);          \
+    if (!status_.ok()) return;            \
   } while (false)
 
 // Not cached.
 void SCIPInterface::SetOptimizationDirection(bool maximize) {
   RETURN_IF_ALREADY_IN_ERROR_STATE;
   InvalidateSolutionSynchronization();
-  RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
-  RETURN_IF_SCIP_ERROR(SCIPsetObjsense(
+  RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+  RETURN_AND_STORE_IF_SCIP_ERROR(SCIPsetObjsense(
       scip_, maximize ? SCIP_OBJSENSE_MAXIMIZE : SCIP_OBJSENSE_MINIMIZE));
 }
 
@@ -253,9 +229,11 @@ void SCIPInterface::SetVariableBounds(int var_index, double lb, double ub) {
   if (variable_is_extracted(var_index)) {
     // Not cached if the variable has been extracted.
     DCHECK_LT(var_index, last_variable_index_);
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
-    RETURN_IF_SCIP_ERROR(SCIPchgVarLb(scip_, scip_variables_[var_index], lb));
-    RETURN_IF_SCIP_ERROR(SCIPchgVarUb(scip_, scip_variables_[var_index], ub));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(
+        SCIPchgVarLb(scip_, scip_variables_[var_index], lb));
+    RETURN_AND_STORE_IF_SCIP_ERROR(
+        SCIPchgVarUb(scip_, scip_variables_[var_index], ub));
   } else {
     sync_status_ = MUST_RELOAD;
   }
@@ -266,14 +244,14 @@ void SCIPInterface::SetVariableInteger(int var_index, bool integer) {
   InvalidateSolutionSynchronization();
   if (variable_is_extracted(var_index)) {
     // Not cached if the variable has been extracted.
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
 #if (SCIP_VERSION >= 210)
     SCIP_Bool infeasible = false;
-    RETURN_IF_SCIP_ERROR(SCIPchgVarType(
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPchgVarType(
         scip_, scip_variables_[var_index],
         integer ? SCIP_VARTYPE_INTEGER : SCIP_VARTYPE_CONTINUOUS, &infeasible));
 #else
-    RETURN_IF_SCIP_ERROR(SCIPchgVarType(
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPchgVarType(
         scip_, scip_variables_[var_index],
         integer ? SCIP_VARTYPE_INTEGER : SCIP_VARTYPE_CONTINUOUS));
 #endif  // SCIP_VERSION >= 210
@@ -288,9 +266,11 @@ void SCIPInterface::SetConstraintBounds(int index, double lb, double ub) {
   if (constraint_is_extracted(index)) {
     // Not cached if the row has been extracted.
     DCHECK_LT(index, last_constraint_index_);
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
-    RETURN_IF_SCIP_ERROR(SCIPchgLhsLinear(scip_, scip_constraints_[index], lb));
-    RETURN_IF_SCIP_ERROR(SCIPchgRhsLinear(scip_, scip_constraints_[index], ub));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(
+        SCIPchgLhsLinear(scip_, scip_constraints_[index], lb));
+    RETURN_AND_STORE_IF_SCIP_ERROR(
+        SCIPchgRhsLinear(scip_, scip_constraints_[index], ub));
   } else {
     sync_status_ = MUST_RELOAD;
   }
@@ -309,8 +289,8 @@ void SCIPInterface::SetCoefficient(MPConstraint* constraint,
     DCHECK_LT(variable->index(), last_variable_index_);
     // SCIP does not allow to set a coefficient directly, so we add the
     // difference between the new and the old value instead.
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
-    RETURN_IF_SCIP_ERROR(SCIPaddCoefLinear(
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddCoefLinear(
         scip_, scip_constraints_[constraint->index()],
         scip_variables_[variable->index()], new_value - old_value));
   } else {
@@ -331,9 +311,9 @@ void SCIPInterface::ClearConstraint(MPConstraint* constraint) {
     const int var_index = entry.first->index();
     const double old_coef_value = entry.second;
     DCHECK(variable_is_extracted(var_index));
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
     // Set coefficient to zero by substracting the old coefficient value.
-    RETURN_IF_SCIP_ERROR(
+    RETURN_AND_STORE_IF_SCIP_ERROR(
         SCIPaddCoefLinear(scip_, scip_constraints_[constraint_index],
                           scip_variables_[var_index], -old_coef_value));
   }
@@ -353,8 +333,10 @@ void SCIPInterface::SetObjectiveOffset(double value) {
 // Clear objective of all its terms.
 void SCIPInterface::ClearObjective() {
   RETURN_IF_ALREADY_IN_ERROR_STATE;
+  sync_status_ = MUST_RELOAD;
+
   InvalidateSolutionSynchronization();
-  RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+  RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
   // Clear linear terms
   for (const auto& entry : solver_->objective_->coefficients_) {
     const int var_index = entry.first->index();
@@ -362,12 +344,15 @@ void SCIPInterface::ClearObjective() {
     if (!variable_is_extracted(var_index)) {
       DCHECK_NE(MODEL_SYNCHRONIZED, sync_status_);
     } else {
-      RETURN_IF_SCIP_ERROR(
+      RETURN_AND_STORE_IF_SCIP_ERROR(
           SCIPchgVarObj(scip_, scip_variables_[var_index], 0.0));
     }
   }
-  // Constant term: change objective offset variable.
-  RETURN_IF_SCIP_ERROR(SCIPchgVarObj(scip_, objective_offset_variable_, 0.0));
+  // Note: we don't clear the objective offset here because it's not necessary
+  // (it's always reset anyway in ExtractObjective) and we sometimes run into
+  // crashes when clearing the whole model (see
+  // http://test/OCL:253365573:BASE:253566457:1560777456754:e181f4ab).
+  // It's not worth to spend time investigating this issue.
 }
 
 void SCIPInterface::BranchingPriorityChangedForVariable(int var_index) {
@@ -397,7 +382,7 @@ void SCIPInterface::ExtractNewVariables() {
   RETURN_IF_ALREADY_IN_ERROR_STATE;
   int total_num_vars = solver_->variables_.size();
   if (total_num_vars > last_variable_index_) {
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
     // Define new variables
     for (int j = last_variable_index_; j < total_num_vars; ++j) {
       MPVariable* const var = solver_->variables_[j];
@@ -406,17 +391,17 @@ void SCIPInterface::ExtractNewVariables() {
       SCIP_VAR* scip_var = nullptr;
       // The true objective coefficient will be set later in ExtractObjective.
       double tmp_obj_coef = 0.0;
-      RETURN_IF_SCIP_ERROR(SCIPcreateVar(
+      RETURN_AND_STORE_IF_SCIP_ERROR(SCIPcreateVar(
           scip_, &scip_var, var->name().c_str(), var->lb(), var->ub(),
           tmp_obj_coef,
           var->integer() ? SCIP_VARTYPE_INTEGER : SCIP_VARTYPE_CONTINUOUS, true,
           false, nullptr, nullptr, nullptr, nullptr, nullptr));
-      RETURN_IF_SCIP_ERROR(SCIPaddVar(scip_, scip_var));
+      RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddVar(scip_, scip_var));
       scip_variables_.push_back(scip_var);
       const int branching_priority = var->branching_priority();
       if (branching_priority != 0) {
         const int index = var->index();
-        RETURN_IF_SCIP_ERROR(SCIPchgVarBranchPriority(
+        RETURN_AND_STORE_IF_SCIP_ERROR(SCIPchgVarBranchPriority(
             scip_, scip_variables_[index], branching_priority));
       }
     }
@@ -429,9 +414,9 @@ void SCIPInterface::ExtractNewVariables() {
         if (var_index >= last_variable_index_) {
           // The variable is new, so we know the previous coefficient
           // value was 0 and we can directly add the coefficient.
-          RETURN_IF_SCIP_ERROR(SCIPaddCoefLinear(scip_, scip_constraints_[i],
-                                                 scip_variables_[var_index],
-                                                 entry.second));
+          RETURN_AND_STORE_IF_SCIP_ERROR(
+              SCIPaddCoefLinear(scip_, scip_constraints_[i],
+                                scip_variables_[var_index], entry.second));
         }
       }
     }
@@ -442,7 +427,7 @@ void SCIPInterface::ExtractNewConstraints() {
   RETURN_IF_ALREADY_IN_ERROR_STATE;
   int total_num_rows = solver_->constraints_.size();
   if (last_constraint_index_ < total_num_rows) {
-    RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+    RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
     // Find the length of the longest row.
     int max_row_length = 0;
     for (int i = last_constraint_index_; i < total_num_rows; ++i) {
@@ -475,12 +460,12 @@ void SCIPInterface::ExtractNewConstraints() {
         DCHECK(variable_is_extracted(ind_index));
         SCIP_VAR* ind_var = scip_variables_[ind_index];
         if (ct->indicator_value() == 0) {
-          RETURN_IF_SCIP_ERROR(
+          RETURN_AND_STORE_IF_SCIP_ERROR(
               SCIPgetNegatedVar(scip_, scip_variables_[ind_index], &ind_var));
         }
 
         if (ct->ub() < std::numeric_limits<double>::infinity()) {
-          RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator(
+          RETURN_AND_STORE_IF_SCIP_ERROR(SCIPcreateConsIndicator(
               scip_, &scip_constraint, ct->name().c_str(), ind_var, size,
               vars.get(), coeffs.get(), ct->ub(),
               /*initial=*/!is_lazy,
@@ -492,14 +477,14 @@ void SCIPInterface::ExtractNewConstraints() {
               /*dynamic=*/false,
               /*removable=*/is_lazy,
               /*stickingatnode=*/false));
-          RETURN_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
+          RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
           scip_constraints_.push_back(scip_constraint);
         }
         if (ct->lb() > -std::numeric_limits<double>::infinity()) {
           for (int i = 0; i < size; ++i) {
             coeffs[i] *= -1;
           }
-          RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator(
+          RETURN_AND_STORE_IF_SCIP_ERROR(SCIPcreateConsIndicator(
               scip_, &scip_constraint, ct->name().c_str(), ind_var, size,
               vars.get(), coeffs.get(), -ct->lb(),
               /*initial=*/!is_lazy,
@@ -511,14 +496,14 @@ void SCIPInterface::ExtractNewConstraints() {
               /*dynamic=*/false,
               /*removable=*/is_lazy,
               /*stickingatnode=*/false));
-          RETURN_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
+          RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
           scip_constraints_.push_back(scip_constraint);
         }
       } else {
         // See
         // http://scip.zib.de/doc/html/cons__linear_8h.php#aa7aed137a4130b35b168812414413481
         // for an explanation of the parameters.
-        RETURN_IF_SCIP_ERROR(SCIPcreateConsLinear(
+        RETURN_AND_STORE_IF_SCIP_ERROR(SCIPcreateConsLinear(
             scip_, &scip_constraint, ct->name().c_str(), size, vars.get(),
             coeffs.get(), ct->lb(), ct->ub(),
             /*initial=*/!is_lazy,
@@ -531,7 +516,7 @@ void SCIPInterface::ExtractNewConstraints() {
             /*dynamic=*/false,
             /*removable=*/is_lazy,
             /*stickingatnode=*/false));
-        RETURN_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
+        RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddCons(scip_, scip_constraint));
         scip_constraints_.push_back(scip_constraint);
       }
     }
@@ -540,19 +525,19 @@ void SCIPInterface::ExtractNewConstraints() {
 
 void SCIPInterface::ExtractObjective() {
   RETURN_IF_ALREADY_IN_ERROR_STATE;
-  RETURN_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
+  RETURN_AND_STORE_IF_SCIP_ERROR(SCIPfreeTransform(scip_));
   // Linear objective: set objective coefficients for all variables (some might
   // have been modified).
   for (const auto& entry : solver_->objective_->coefficients_) {
     const int var_index = entry.first->index();
     const double obj_coef = entry.second;
-    RETURN_IF_SCIP_ERROR(
+    RETURN_AND_STORE_IF_SCIP_ERROR(
         SCIPchgVarObj(scip_, scip_variables_[var_index], obj_coef));
   }
 
-  // Constant term: change objective offset variable.
-  RETURN_IF_SCIP_ERROR(SCIPchgVarObj(scip_, objective_offset_variable_,
-                                     solver_->Objective().offset()));
+  // Constant term: change objective offset.
+  RETURN_AND_STORE_IF_SCIP_ERROR(SCIPaddOrigObjoffset(
+      scip_, solver_->Objective().offset() - SCIPgetOrigObjoffset(scip_)));
 }
 
 #define RETURN_ABNORMAL_IF_BAD_STATUS             \
@@ -567,7 +552,7 @@ void SCIPInterface::ExtractObjective() {
 #define RETURN_ABNORMAL_IF_SCIP_ERROR(x) \
   do {                                   \
     RETURN_ABNORMAL_IF_BAD_STATUS;       \
-    status_ = TO_STATUS(x);              \
+    status_ = SCIP_TO_STATUS(x);         \
     RETURN_ABNORMAL_IF_BAD_STATUS;       \
   } while (false);
 
@@ -638,11 +623,6 @@ MPSolver::ResultStatus SCIPInterface::Solve(const MPSolverParameters& param) {
       RETURN_ABNORMAL_IF_SCIP_ERROR(SCIPcreateSol(scip_, &solution, nullptr));
     }
 
-    // The variable representing the objective offset should always be one!!
-    // See CreateSCIP().
-    RETURN_ABNORMAL_IF_SCIP_ERROR(
-        SCIPsetSolVal(scip_, solution, objective_offset_variable_, 1.0));
-
     // Fill the other variables from the given solution hint.
     for (const std::pair<const MPVariable*, double>& p :
          solver_->solution_hint_) {
@@ -740,6 +720,23 @@ MPSolver::ResultStatus SCIPInterface::Solve(const MPSolverParameters& param) {
   return result_status_;
 }
 
+absl::optional<MPSolutionResponse> SCIPInterface::DirectlySolveProto(
+    const MPModelRequest& request) {
+  const auto status_or = ScipSolveProto(request);
+  if (status_or.ok()) return status_or.ValueOrDie();
+  // Special case: if something is not implemented yet, fall back to solving
+  // through MPSolver.
+  if (util::IsUnimplemented(status_or.status())) return absl::nullopt;
+
+  if (request.enable_internal_solver_output()) {
+    LOG(INFO) << "Invalid SCIP status: " << status_or.status();
+  }
+  MPSolutionResponse response;
+  response.set_status(MPSOLVER_NOT_SOLVED);
+  response.set_status_str(status_or.status().ToString());
+  return response;
+}
+
 int64 SCIPInterface::iterations() const {
   // NOTE(user): As of 2018-12 it doesn't run in the stubby server, and is
   // a specialized call, so it's ok to crash if the status is broken.
@@ -785,7 +782,8 @@ void SCIPInterface::SetRelativeMipGap(double value) {
   //   won't crash even if we're in an error state. I did *not* verify this).
   // - if that call yielded an error *and* we weren't already in an error state,
   //   set the state to that error we just got.
-  const auto status = TO_STATUS(SCIPsetRealParam(scip_, "limits/gap", value));
+  const auto status =
+      SCIP_TO_STATUS(SCIPsetRealParam(scip_, "limits/gap", value));
   if (status_.ok()) status_ = status;
 }
 
@@ -799,18 +797,18 @@ void SCIPInterface::SetPrimalTolerance(double value) {
   if (value < current_lpfeastol) {
     // See the NOTE on SetRelativeMipGap().
     const auto status =
-        TO_STATUS(SCIPsetRealParam(scip_, "numerics/lpfeastol", value));
+        SCIP_TO_STATUS(SCIPsetRealParam(scip_, "numerics/lpfeastol", value));
     if (status_.ok()) status_ = status;
   }
   // See the NOTE on SetRelativeMipGap().
   const auto status =
-      TO_STATUS(SCIPsetRealParam(scip_, "numerics/feastol", value));
+      SCIP_TO_STATUS(SCIPsetRealParam(scip_, "numerics/feastol", value));
   if (status_.ok()) status_ = status;
 }
 
 void SCIPInterface::SetDualTolerance(double value) {
   const auto status =
-      TO_STATUS(SCIPsetRealParam(scip_, "numerics/dualfeastol", value));
+      SCIP_TO_STATUS(SCIPsetRealParam(scip_, "numerics/dualfeastol", value));
   if (status_.ok()) status_ = status;
 }
 
@@ -819,13 +817,13 @@ void SCIPInterface::SetPresolveMode(int value) {
   switch (value) {
     case MPSolverParameters::PRESOLVE_OFF: {
       const auto status =
-          TO_STATUS(SCIPsetIntParam(scip_, "presolving/maxrounds", 0));
+          SCIP_TO_STATUS(SCIPsetIntParam(scip_, "presolving/maxrounds", 0));
       if (status_.ok()) status_ = status;
       return;
     }
     case MPSolverParameters::PRESOLVE_ON: {
       const auto status =
-          TO_STATUS(SCIPsetIntParam(scip_, "presolving/maxrounds", -1));
+          SCIP_TO_STATUS(SCIPsetIntParam(scip_, "presolving/maxrounds", -1));
       if (status_.ok()) status_ = status;
       return;
     }
@@ -848,20 +846,20 @@ void SCIPInterface::SetLpAlgorithm(int value) {
   switch (value) {
     case MPSolverParameters::DUAL: {
       const auto status =
-          TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'd'));
+          SCIP_TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'd'));
       if (status_.ok()) status_ = status;
       return;
     }
     case MPSolverParameters::PRIMAL: {
       const auto status =
-          TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'p'));
+          SCIP_TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'p'));
       if (status_.ok()) status_ = status;
       return;
     }
     case MPSolverParameters::BARRIER: {
       // Barrier with crossover.
       const auto status =
-          TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'p'));
+          SCIP_TO_STATUS(SCIPsetCharParam(scip_, "lp/initalgorithm", 'p'));
       if (status_.ok()) status_ = status;
       return;
     }
@@ -917,8 +915,7 @@ MPSolverInterface* BuildSCIPInterface(MPSolver* const solver) {
 }  // namespace operations_research
 #endif  //  #if defined(USE_SCIP)
 
-#undef TO_STATUS
-#undef RETURN_IF_SCIP_ERROR
+#undef RETURN_AND_STORE_IF_SCIP_ERROR
 #undef RETURN_IF_ALREADY_IN_ERROR_STATE
 #undef RETURN_ABNORMAL_IF_BAD_STATUS
 #undef RETURN_ABNORMAL_IF_SCIP_ERROR
diff --git a/ortools/linear_solver/scip_proto_solver.h b/ortools/linear_solver/scip_proto_solver.h
new file mode 100644
index 0000000000..065ed95e3e
--- /dev/null
+++ b/ortools/linear_solver/scip_proto_solver.h
@@ -0,0 +1,27 @@
+// Copyright 2010-2018 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_LINEAR_SOLVER_SCIP_PROTO_SOLVER_H_
+#define OR_TOOLS_LINEAR_SOLVER_SCIP_PROTO_SOLVER_H_
+
+#include "ortools/base/statusor.h"
+#include "ortools/linear_solver/linear_solver.pb.h"
+
+namespace operations_research {
+
+util::StatusOr<MPSolutionResponse> ScipSolveProto(
+    const MPModelRequest& request);
+
+}  // namespace operations_research
+
+#endif  // OR_TOOLS_LINEAR_SOLVER_SCIP_PROTO_SOLVER_H_