1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 //     http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13 
14 // Gurobi backend to MPSolver.
15 //
16 // Implementation Notes:
17 //
18 // Incrementalism (last updated June 29, 2020): For solving both LPs and MIPs,
19 // Gurobi attempts to reuse information from previous solves, potentially
20 // giving a faster solve time. MPSolver supports this for the following problem
21 // modification types:
22 //   * Adding a variable,
23 //   * Adding a linear constraint,
24 //   * Updating a variable bound,
25 //   * Updating an objective coefficient or the objective offset (note that in
26 //     Gurobi 7.5 LP solver, there is a bug if you update only the objective
27 //     offset and nothing else).
28 //   * Updating a coefficient in the constraint matrix.
29 //   * Updating the type of variable (integer, continuous)
30 //   * Changing the optimization direction.
31 // Updates of the following types will force a resolve from scratch:
32 //   * Updating the upper or lower bounds of a linear constraint. Note that in
33 //     MPSolver's model, this includes updating the sense (le, ge, eq, range) of
34 //     a linear constraint.
35 //   * Clearing a constraint
36 // Any model containing indicator constraints is considered "non-incremental"
37 // and will always solve from scratch.
38 //
39 // The above limitations are largely due MPSolver and this file, not Gurobi.
40 //
41 // Warning(rander): the interactions between callbacks and incrementalism are
42 // poorly tested, proceed with caution.
43 //
44 
45 #include <cmath>
46 #include <cstddef>
47 #include <cstdint>
48 #include <limits>
49 #include <memory>
50 #include <stdexcept>
51 #include <string>
52 #include <utility>
53 #include <vector>
54 
55 #include "absl/base/attributes.h"
56 #include "absl/status/status.h"
57 #include "absl/strings/match.h"
58 #include "absl/strings/str_format.h"
59 #include "ortools/base/commandlineflags.h"
60 #include "ortools/base/integral_types.h"
61 #include "ortools/base/logging.h"
62 #include "ortools/base/map_util.h"
63 #include "ortools/base/timer.h"
64 #include "ortools/gurobi/environment.h"
65 #include "ortools/linear_solver/gurobi_proto_solver.h"
66 #include "ortools/linear_solver/linear_solver.h"
67 #include "ortools/linear_solver/linear_solver_callback.h"
68 #include "ortools/util/time_limit.h"
69 
70 ABSL_FLAG(int, num_gurobi_threads, 4,
71           "Number of threads available for Gurobi.");
72 
73 namespace operations_research {
74 
75 class GurobiInterface : public MPSolverInterface {
76  public:
77   // Constructor that takes a name for the underlying GRB solver.
78   explicit GurobiInterface(MPSolver* const solver, bool mip);
79   ~GurobiInterface() override;
80 
81   // Sets the optimization direction (min/max).
82   void SetOptimizationDirection(bool maximize) override;
83 
84   // ----- Solve -----
85   // Solves the problem using the parameter values specified.
86   MPSolver::ResultStatus Solve(const MPSolverParameters& param) override;
87   absl::optional<MPSolutionResponse> DirectlySolveProto(
88       const MPModelRequest& request, std::atomic<bool>* interrupt) override;
89   // Writes the model.
90   void Write(const std::string& filename) override;
91 
92   // ----- Model modifications and extraction -----
93   // Resets extracted model
94   void Reset() override;
95 
96   // Modifies bounds.
97   void SetVariableBounds(int var_index, double lb, double ub) override;
98   void SetVariableInteger(int var_index, bool integer) override;
99   void SetConstraintBounds(int row_index, double lb, double ub) override;
100 
101   // Adds Constraint incrementally.
102   void AddRowConstraint(MPConstraint* const ct) override;
103   bool AddIndicatorConstraint(MPConstraint* const ct) override;
104   // Adds variable incrementally.
105   void AddVariable(MPVariable* const var) override;
106   // Changes a coefficient in a constraint.
107   void SetCoefficient(MPConstraint* const constraint,
108                       const MPVariable* const variable, double new_value,
109                       double old_value) override;
110   // Clears a constraint from all its terms.
111   void ClearConstraint(MPConstraint* const constraint) override;
112   // Changes a coefficient in the linear objective
113   void SetObjectiveCoefficient(const MPVariable* const variable,
114                                double coefficient) override;
115   // Changes the constant term in the linear objective.
116   void SetObjectiveOffset(double value) override;
117   // Clears the objective from all its terms.
118   void ClearObjective() override;
119   void BranchingPriorityChangedForVariable(int var_index) override;
120 
121   // ------ Query statistics on the solution and the solve ------
122   // Number of simplex or interior-point iterations
123   int64_t iterations() const override;
124   // Number of branch-and-bound nodes. Only available for discrete problems.
125   int64_t nodes() const override;
126 
127   // Returns the basis status of a row.
128   MPSolver::BasisStatus row_status(int constraint_index) const override;
129   // Returns the basis status of a column.
130   MPSolver::BasisStatus column_status(int variable_index) const override;
131 
132   // ----- Misc -----
133   // Queries problem type.
IsContinuous() const134   bool IsContinuous() const override { return IsLP(); }
IsLP() const135   bool IsLP() const override { return !mip_; }
IsMIP() const136   bool IsMIP() const override { return mip_; }
137 
138   void ExtractNewVariables() override;
139   void ExtractNewConstraints() override;
140   void ExtractObjective() override;
141 
SolverVersion() const142   std::string SolverVersion() const override {
143     int major, minor, technical;
144     GRBversion(&major, &minor, &technical);
145     return absl::StrFormat("Gurobi library version %d.%d.%d\n", major, minor,
146                            technical);
147   }
148 
InterruptSolve()149   bool InterruptSolve() override {
150     const absl::MutexLock lock(&hold_interruptions_mutex_);
151     if (model_ != nullptr) GRBterminate(model_);
152     return true;
153   }
154 
underlying_solver()155   void* underlying_solver() override { return reinterpret_cast<void*>(model_); }
156 
ComputeExactConditionNumber() const157   double ComputeExactConditionNumber() const override {
158     if (!IsContinuous()) {
159       LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
160                   << " GUROBI_MIXED_INTEGER_PROGRAMMING";
161       return 0.0;
162     }
163 
164     // TODO(user): Not yet working.
165     LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
166                 << " GUROBI_LINEAR_PROGRAMMING";
167     return 0.0;
168 
169     // double cond = 0.0;
170     // const int status = GRBgetdblattr(model_, GRB_DBL_ATTR_KAPPA, &cond);
171     // if (0 == status) {
172     //   return cond;
173     // } else {
174     //   LOG(DFATAL) << "Condition number only available for "
175     //               << "continuous problems";
176     //   return 0.0;
177     // }
178   }
179 
180   // Iterates through the solutions in Gurobi's solution pool.
181   bool NextSolution() override;
182 
183   void SetCallback(MPCallback* mp_callback) override;
SupportsCallbacks() const184   bool SupportsCallbacks() const override { return true; }
185 
186  private:
187   // Sets all parameters in the underlying solver.
188   void SetParameters(const MPSolverParameters& param) override;
189   // Sets solver-specific parameters (avoiding using files). The previous
190   // implementations supported multi-line strings of the form:
191   // parameter_i value_i\n
192   // We extend support for strings of the form:
193   // parameter1=value1,....,parametern=valuen
194   // or for strings of the form:
195   // parameter1 value1, ... ,parametern valuen
196   // which are easier to set in the command line.
197   // This implementations relies on SetSolverSpecificParameters, which has the
198   // extra benefit of unifying the way we handle specific parameters for both
199   // proto-based solves and for MPModel solves.
200   bool SetSolverSpecificParametersAsString(
201       const std::string& parameters) override;
202   // Sets each parameter in the underlying solver.
203   void SetRelativeMipGap(double value) override;
204   void SetPrimalTolerance(double value) override;
205   void SetDualTolerance(double value) override;
206   void SetPresolveMode(int value) override;
207   void SetScalingMode(int value) override;
208   void SetLpAlgorithm(int value) override;
209 
210   MPSolver::BasisStatus TransformGRBVarBasisStatus(
211       int gurobi_basis_status) const;
212   MPSolver::BasisStatus TransformGRBConstraintBasisStatus(
213       int gurobi_basis_status, int constraint_index) const;
214 
215   // See the implementation note at the top of file on incrementalism.
216   bool ModelIsNonincremental() const;
217 
218   void SetIntAttr(const char* name, int value);
219   int GetIntAttr(const char* name) const;
220   void SetDoubleAttr(const char* name, double value);
221   double GetDoubleAttr(const char* name) const;
222   void SetIntAttrElement(const char* name, int index, int value);
223   int GetIntAttrElement(const char* name, int index) const;
224   void SetDoubleAttrElement(const char* name, int index, double value);
225   double GetDoubleAttrElement(const char* name, int index) const;
226   std::vector<double> GetDoubleAttrArray(const char* name, int elements);
227   void SetCharAttrElement(const char* name, int index, char value);
228   char GetCharAttrElement(const char* name, int index) const;
229 
230   void CheckedGurobiCall(int err) const;
231 
232   int SolutionCount() const;
233 
234   GRBmodel* model_;
235   GRBenv* env_;
236   bool mip_;
237   int current_solution_index_;
238   MPCallback* callback_ = nullptr;
239   bool update_branching_priorities_ = false;
240   // Has length equal to the number of MPVariables in
241   // MPSolverInterface::solver_. Values are the index of the corresponding
242   // Gurobi variable. Note that Gurobi may have additional auxiliary variables
243   // not represented by MPVariables, such as those created by two-sided range
244   // constraints.
245   std::vector<int> mp_var_to_gurobi_var_;
246   // Has length equal to the number of MPConstraints in
247   // MPSolverInterface::solver_. Values are the index of the corresponding
248   // linear (or range) constraint in Gurobi, or -1 if no such constraint exists
249   // (e.g. for indicator constraints).
250   std::vector<int> mp_cons_to_gurobi_linear_cons_;
251   // Should match the Gurobi model after it is updated.
252   int num_gurobi_vars_ = 0;
253   // Should match the Gurobi model after it is updated.
254   // NOTE(user): indicator constraints are not counted below.
255   int num_gurobi_linear_cons_ = 0;
256   // See the implementation note at the top of file on incrementalism.
257   bool had_nonincremental_change_ = false;
258 
259   // Mutex is held to prevent InterruptSolve() to call GRBterminate() when
260   // model_ is not completely built. It also prevents model_ to be changed
261   // during the execution of GRBterminate().
262   mutable absl::Mutex hold_interruptions_mutex_;
263 };
264 
265 namespace {
266 
CheckedGurobiCall(int err,GRBenv * const env)267 void CheckedGurobiCall(int err, GRBenv* const env) {
268   CHECK_EQ(0, err) << "Fatal error with code " << err << ", due to "
269                    << GRBgeterrormsg(env);
270 }
271 
272 // For interacting directly with the Gurobi C API for callbacks.
273 struct GurobiInternalCallbackContext {
274   GRBmodel* model;
275   void* gurobi_internal_callback_data;
276   int where;
277 };
278 
279 class GurobiMPCallbackContext : public MPCallbackContext {
280  public:
281   GurobiMPCallbackContext(GRBenv* env,
282                           const std::vector<int>* mp_var_to_gurobi_var,
283                           int num_gurobi_vars, bool might_add_cuts,
284                           bool might_add_lazy_constraints);
285 
286   // Implementation of the interface.
287   MPCallbackEvent Event() override;
288   bool CanQueryVariableValues() override;
289   double VariableValue(const MPVariable* variable) override;
290   void AddCut(const LinearRange& cutting_plane) override;
291   void AddLazyConstraint(const LinearRange& lazy_constraint) override;
292   double SuggestSolution(
293       const absl::flat_hash_map<const MPVariable*, double>& solution) override;
294   int64_t NumExploredNodes() override;
295 
296   // Call this method to update the internal state of the callback context
297   // before passing it to MPCallback::RunCallback().
298   void UpdateFromGurobiState(
299       const GurobiInternalCallbackContext& gurobi_internal_context);
300 
301  private:
302   // Wraps GRBcbget(), used to query the state of the solver.  See
303   // http://www.gurobi.com/documentation/8.0/refman/callback_codes.html#sec:CallbackCodes
304   // for callback_code values.
305   template <typename T>
306   T GurobiCallbackGet(
307       const GurobiInternalCallbackContext& gurobi_internal_context,
308       int callback_code);
309   void CheckedGurobiCall(int gurobi_error_code) const;
310 
311   template <typename GRBConstraintFunction>
312   void AddGeneratedConstraint(const LinearRange& linear_range,
313                               GRBConstraintFunction grb_constraint_function);
314 
315   GRBenv* const env_;
316   const std::vector<int>* const mp_var_to_gurobi_var_;
317   const int num_gurobi_vars_;
318 
319   const bool might_add_cuts_;
320   const bool might_add_lazy_constraints_;
321 
322   // Stateful, updated before each call to the callback.
323   GurobiInternalCallbackContext current_gurobi_internal_callback_context_;
324   bool variable_values_extracted_ = false;
325   std::vector<double> gurobi_variable_values_;
326 };
327 
CheckedGurobiCall(int gurobi_error_code) const328 void GurobiMPCallbackContext::CheckedGurobiCall(int gurobi_error_code) const {
329   ::operations_research::CheckedGurobiCall(gurobi_error_code, env_);
330 }
331 
GurobiMPCallbackContext(GRBenv * env,const std::vector<int> * mp_var_to_gurobi_var,int num_gurobi_vars,bool might_add_cuts,bool might_add_lazy_constraints)332 GurobiMPCallbackContext::GurobiMPCallbackContext(
333     GRBenv* env, const std::vector<int>* mp_var_to_gurobi_var,
334     int num_gurobi_vars, bool might_add_cuts, bool might_add_lazy_constraints)
335     : env_(ABSL_DIE_IF_NULL(env)),
336       mp_var_to_gurobi_var_(ABSL_DIE_IF_NULL(mp_var_to_gurobi_var)),
337       num_gurobi_vars_(num_gurobi_vars),
338       might_add_cuts_(might_add_cuts),
339       might_add_lazy_constraints_(might_add_lazy_constraints) {}
340 
UpdateFromGurobiState(const GurobiInternalCallbackContext & gurobi_internal_context)341 void GurobiMPCallbackContext::UpdateFromGurobiState(
342     const GurobiInternalCallbackContext& gurobi_internal_context) {
343   current_gurobi_internal_callback_context_ = gurobi_internal_context;
344   variable_values_extracted_ = false;
345 }
346 
NumExploredNodes()347 int64_t GurobiMPCallbackContext::NumExploredNodes() {
348   switch (Event()) {
349     case MPCallbackEvent::kMipNode:
350       return static_cast<int64_t>(GurobiCallbackGet<double>(
351           current_gurobi_internal_callback_context_, GRB_CB_MIPNODE_NODCNT));
352     case MPCallbackEvent::kMipSolution:
353       return static_cast<int64_t>(GurobiCallbackGet<double>(
354           current_gurobi_internal_callback_context_, GRB_CB_MIPSOL_NODCNT));
355     default:
356       LOG(FATAL) << "Node count is supported only for callback events MIP_NODE "
357                     "and MIP_SOL, but was requested at: "
358                  << ToString(Event());
359   }
360 }
361 
362 template <typename T>
GurobiCallbackGet(const GurobiInternalCallbackContext & gurobi_internal_context,const int callback_code)363 T GurobiMPCallbackContext::GurobiCallbackGet(
364     const GurobiInternalCallbackContext& gurobi_internal_context,
365     const int callback_code) {
366   T result = 0;
367   CheckedGurobiCall(
368       GRBcbget(gurobi_internal_context.gurobi_internal_callback_data,
369                gurobi_internal_context.where, callback_code,
370                static_cast<void*>(&result)));
371   return result;
372 }
373 
Event()374 MPCallbackEvent GurobiMPCallbackContext::Event() {
375   switch (current_gurobi_internal_callback_context_.where) {
376     case GRB_CB_POLLING:
377       return MPCallbackEvent::kPolling;
378     case GRB_CB_PRESOLVE:
379       return MPCallbackEvent::kPresolve;
380     case GRB_CB_SIMPLEX:
381       return MPCallbackEvent::kSimplex;
382     case GRB_CB_MIP:
383       return MPCallbackEvent::kMip;
384     case GRB_CB_MIPSOL:
385       return MPCallbackEvent::kMipSolution;
386     case GRB_CB_MIPNODE:
387       return MPCallbackEvent::kMipNode;
388     case GRB_CB_MESSAGE:
389       return MPCallbackEvent::kMessage;
390     case GRB_CB_BARRIER:
391       return MPCallbackEvent::kBarrier;
392       // TODO(b/112427356): in Gurobi 8.0, there is a new callback location.
393       // case GRB_CB_MULTIOBJ:
394       //   return MPCallbackEvent::kMultiObj;
395     default:
396       LOG_FIRST_N(ERROR, 1) << "Gurobi callback at unknown where="
397                             << current_gurobi_internal_callback_context_.where;
398       return MPCallbackEvent::kUnknown;
399   }
400 }
401 
CanQueryVariableValues()402 bool GurobiMPCallbackContext::CanQueryVariableValues() {
403   const MPCallbackEvent where = Event();
404   if (where == MPCallbackEvent::kMipSolution) {
405     return true;
406   }
407   if (where == MPCallbackEvent::kMipNode) {
408     const int gurobi_node_status = GurobiCallbackGet<int>(
409         current_gurobi_internal_callback_context_, GRB_CB_MIPNODE_STATUS);
410     return gurobi_node_status == GRB_OPTIMAL;
411   }
412   return false;
413 }
414 
VariableValue(const MPVariable * variable)415 double GurobiMPCallbackContext::VariableValue(const MPVariable* variable) {
416   CHECK(variable != nullptr);
417   if (!variable_values_extracted_) {
418     const MPCallbackEvent where = Event();
419     CHECK(where == MPCallbackEvent::kMipSolution ||
420           where == MPCallbackEvent::kMipNode)
421         << "You can only call VariableValue at "
422         << ToString(MPCallbackEvent::kMipSolution) << " or "
423         << ToString(MPCallbackEvent::kMipNode)
424         << " but called from: " << ToString(where);
425     const int gurobi_get_var_param = where == MPCallbackEvent::kMipNode
426                                          ? GRB_CB_MIPNODE_REL
427                                          : GRB_CB_MIPSOL_SOL;
428 
429     gurobi_variable_values_.resize(num_gurobi_vars_);
430     CheckedGurobiCall(GRBcbget(
431         current_gurobi_internal_callback_context_.gurobi_internal_callback_data,
432         current_gurobi_internal_callback_context_.where, gurobi_get_var_param,
433         static_cast<void*>(gurobi_variable_values_.data())));
434     variable_values_extracted_ = true;
435   }
436   return gurobi_variable_values_[mp_var_to_gurobi_var_->at(variable->index())];
437 }
438 
439 template <typename GRBConstraintFunction>
AddGeneratedConstraint(const LinearRange & linear_range,GRBConstraintFunction grb_constraint_function)440 void GurobiMPCallbackContext::AddGeneratedConstraint(
441     const LinearRange& linear_range,
442     GRBConstraintFunction grb_constraint_function) {
443   std::vector<int> variable_indices;
444   std::vector<double> variable_coefficients;
445   const int num_terms = linear_range.linear_expr().terms().size();
446   variable_indices.reserve(num_terms);
447   variable_coefficients.reserve(num_terms);
448   for (const auto& var_coef_pair : linear_range.linear_expr().terms()) {
449     variable_indices.push_back(
450         mp_var_to_gurobi_var_->at(var_coef_pair.first->index()));
451     variable_coefficients.push_back(var_coef_pair.second);
452   }
453   if (std::isfinite(linear_range.upper_bound())) {
454     CheckedGurobiCall(grb_constraint_function(
455         current_gurobi_internal_callback_context_.gurobi_internal_callback_data,
456         variable_indices.size(), variable_indices.data(),
457         variable_coefficients.data(), GRB_LESS_EQUAL,
458         linear_range.upper_bound()));
459   }
460   if (std::isfinite(linear_range.lower_bound())) {
461     CheckedGurobiCall(grb_constraint_function(
462         current_gurobi_internal_callback_context_.gurobi_internal_callback_data,
463         variable_indices.size(), variable_indices.data(),
464         variable_coefficients.data(), GRB_GREATER_EQUAL,
465         linear_range.lower_bound()));
466   }
467 }
468 
AddCut(const LinearRange & cutting_plane)469 void GurobiMPCallbackContext::AddCut(const LinearRange& cutting_plane) {
470   CHECK(might_add_cuts_);
471   const MPCallbackEvent where = Event();
472   CHECK(where == MPCallbackEvent::kMipNode)
473       << "Cuts can only be added at MIP_NODE, tried to add cut at: "
474       << ToString(where);
475   AddGeneratedConstraint(cutting_plane, GRBcbcut);
476 }
477 
AddLazyConstraint(const LinearRange & lazy_constraint)478 void GurobiMPCallbackContext::AddLazyConstraint(
479     const LinearRange& lazy_constraint) {
480   CHECK(might_add_lazy_constraints_);
481   const MPCallbackEvent where = Event();
482   CHECK(where == MPCallbackEvent::kMipNode ||
483         where == MPCallbackEvent::kMipSolution)
484       << "Lazy constraints can only be added at MIP_NODE or MIP_SOL, tried to "
485          "add lazy constraint at: "
486       << ToString(where);
487   AddGeneratedConstraint(lazy_constraint, GRBcblazy);
488 }
489 
SuggestSolution(const absl::flat_hash_map<const MPVariable *,double> & solution)490 double GurobiMPCallbackContext::SuggestSolution(
491     const absl::flat_hash_map<const MPVariable*, double>& solution) {
492   const MPCallbackEvent where = Event();
493   CHECK(where == MPCallbackEvent::kMipNode)
494       << "Feasible solutions can only be added at MIP_NODE, tried to add "
495          "solution at: "
496       << ToString(where);
497 
498   std::vector<double> full_solution(num_gurobi_vars_, GRB_UNDEFINED);
499   for (const auto& variable_value : solution) {
500     const MPVariable* var = variable_value.first;
501     full_solution[mp_var_to_gurobi_var_->at(var->index())] =
502         variable_value.second;
503   }
504 
505   double objval;
506   CheckedGurobiCall(GRBcbsolution(
507       current_gurobi_internal_callback_context_.gurobi_internal_callback_data,
508       full_solution.data(), &objval));
509 
510   return objval;
511 }
512 
513 struct MPCallbackWithGurobiContext {
514   GurobiMPCallbackContext* context;
515   MPCallback* callback;
516 };
517 
518 // NOTE(user): This function must have this exact API, because we are passing
519 // it to Gurobi as a callback.
CallbackImpl(GRBmodel * model,void * gurobi_internal_callback_data,int where,void * raw_model_and_callback)520 int GUROBI_STDCALL CallbackImpl(GRBmodel* model,
521                                 void* gurobi_internal_callback_data, int where,
522                                 void* raw_model_and_callback) {
523   MPCallbackWithGurobiContext* const callback_with_context =
524       static_cast<MPCallbackWithGurobiContext*>(raw_model_and_callback);
525   CHECK(callback_with_context != nullptr);
526   CHECK(callback_with_context->context != nullptr);
527   CHECK(callback_with_context->callback != nullptr);
528   GurobiInternalCallbackContext gurobi_internal_context{
529       model, gurobi_internal_callback_data, where};
530   callback_with_context->context->UpdateFromGurobiState(
531       gurobi_internal_context);
532   callback_with_context->callback->RunCallback(callback_with_context->context);
533   return 0;
534 }
535 
536 }  // namespace
537 
CheckedGurobiCall(int err) const538 void GurobiInterface::CheckedGurobiCall(int err) const {
539   ::operations_research::CheckedGurobiCall(err, env_);
540 }
541 
SetIntAttr(const char * name,int value)542 void GurobiInterface::SetIntAttr(const char* name, int value) {
543   CheckedGurobiCall(GRBsetintattr(model_, name, value));
544 }
545 
GetIntAttr(const char * name) const546 int GurobiInterface::GetIntAttr(const char* name) const {
547   int value;
548   CheckedGurobiCall(GRBgetintattr(model_, name, &value));
549   return value;
550 }
551 
SetDoubleAttr(const char * name,double value)552 void GurobiInterface::SetDoubleAttr(const char* name, double value) {
553   CheckedGurobiCall(GRBsetdblattr(model_, name, value));
554 }
555 
GetDoubleAttr(const char * name) const556 double GurobiInterface::GetDoubleAttr(const char* name) const {
557   double value;
558   CheckedGurobiCall(GRBgetdblattr(model_, name, &value));
559   return value;
560 }
561 
SetIntAttrElement(const char * name,int index,int value)562 void GurobiInterface::SetIntAttrElement(const char* name, int index,
563                                         int value) {
564   CheckedGurobiCall(GRBsetintattrelement(model_, name, index, value));
565 }
566 
GetIntAttrElement(const char * name,int index) const567 int GurobiInterface::GetIntAttrElement(const char* name, int index) const {
568   int value;
569   CheckedGurobiCall(GRBgetintattrelement(model_, name, index, &value));
570   return value;
571 }
572 
SetDoubleAttrElement(const char * name,int index,double value)573 void GurobiInterface::SetDoubleAttrElement(const char* name, int index,
574                                            double value) {
575   CheckedGurobiCall(GRBsetdblattrelement(model_, name, index, value));
576 }
GetDoubleAttrElement(const char * name,int index) const577 double GurobiInterface::GetDoubleAttrElement(const char* name,
578                                              int index) const {
579   double value;
580   CheckedGurobiCall(GRBgetdblattrelement(model_, name, index, &value));
581   return value;
582 }
583 
GetDoubleAttrArray(const char * name,int elements)584 std::vector<double> GurobiInterface::GetDoubleAttrArray(const char* name,
585                                                         int elements) {
586   std::vector<double> results(elements);
587   CheckedGurobiCall(
588       GRBgetdblattrarray(model_, name, 0, elements, results.data()));
589   return results;
590 }
591 
SetCharAttrElement(const char * name,int index,char value)592 void GurobiInterface::SetCharAttrElement(const char* name, int index,
593                                          char value) {
594   CheckedGurobiCall(GRBsetcharattrelement(model_, name, index, value));
595 }
GetCharAttrElement(const char * name,int index) const596 char GurobiInterface::GetCharAttrElement(const char* name, int index) const {
597   char value;
598   CheckedGurobiCall(GRBgetcharattrelement(model_, name, index, &value));
599   return value;
600 }
601 
602 // Creates a LP/MIP instance with the specified name and minimization objective.
GurobiInterface(MPSolver * const solver,bool mip)603 GurobiInterface::GurobiInterface(MPSolver* const solver, bool mip)
604     : MPSolverInterface(solver),
605       model_(nullptr),
606       env_(nullptr),
607       mip_(mip),
608       current_solution_index_(0) {
609   env_ = GetGurobiEnv().value();
610   CheckedGurobiCall(GRBnewmodel(env_, &model_, solver_->name_.c_str(),
611                                 0,          // numvars
612                                 nullptr,    // obj
613                                 nullptr,    // lb
614                                 nullptr,    // ub
615                                 nullptr,    // vtype
616                                 nullptr));  // varnanes
617   SetIntAttr(GRB_INT_ATTR_MODELSENSE, maximize_ ? GRB_MAXIMIZE : GRB_MINIMIZE);
618   CheckedGurobiCall(GRBsetintparam(env_, GRB_INT_PAR_THREADS,
619                                    absl::GetFlag(FLAGS_num_gurobi_threads)));
620 }
621 
~GurobiInterface()622 GurobiInterface::~GurobiInterface() {
623   CheckedGurobiCall(GRBfreemodel(model_));
624   GRBfreeenv(env_);
625 }
626 
627 // ------ Model modifications and extraction -----
628 
Reset()629 void GurobiInterface::Reset() {
630   // We hold calls to GRBterminate() until the new model_ is ready.
631   const absl::MutexLock lock(&hold_interruptions_mutex_);
632 
633   GRBmodel* old_model = model_;
634   CheckedGurobiCall(GRBnewmodel(env_, &model_, solver_->name_.c_str(),
635                                 0,          // numvars
636                                 nullptr,    // obj
637                                 nullptr,    // lb
638                                 nullptr,    // ub
639                                 nullptr,    // vtype
640                                 nullptr));  // varnames
641 
642   // Copy all existing parameters from the previous model to the new one. This
643   // ensures that if a user calls multiple times
644   // SetSolverSpecificParametersAsString() and then Reset() is called, we still
645   // take into account all parameters.
646   //
647   // The current code only reapplies the parameters stored in
648   // solver_specific_parameter_string_ at the start of the solve; other
649   // parameters set by previous calls are only kept in the Gurobi model.
650   CheckedGurobiCall(GRBcopyparams(GRBgetenv(model_), GRBgetenv(old_model)));
651 
652   CheckedGurobiCall(GRBfreemodel(old_model));
653   old_model = nullptr;
654 
655   ResetExtractionInformation();
656   mp_var_to_gurobi_var_.clear();
657   mp_cons_to_gurobi_linear_cons_.clear();
658   num_gurobi_vars_ = 0;
659   num_gurobi_linear_cons_ = 0;
660   had_nonincremental_change_ = false;
661 }
662 
SetOptimizationDirection(bool maximize)663 void GurobiInterface::SetOptimizationDirection(bool maximize) {
664   InvalidateSolutionSynchronization();
665   SetIntAttr(GRB_INT_ATTR_MODELSENSE, maximize_ ? GRB_MAXIMIZE : GRB_MINIMIZE);
666 }
667 
SetVariableBounds(int var_index,double lb,double ub)668 void GurobiInterface::SetVariableBounds(int var_index, double lb, double ub) {
669   InvalidateSolutionSynchronization();
670   if (!had_nonincremental_change_ && variable_is_extracted(var_index)) {
671     SetDoubleAttrElement(GRB_DBL_ATTR_LB, mp_var_to_gurobi_var_.at(var_index),
672                          lb);
673     SetDoubleAttrElement(GRB_DBL_ATTR_UB, mp_var_to_gurobi_var_.at(var_index),
674                          ub);
675   } else {
676     sync_status_ = MUST_RELOAD;
677   }
678 }
679 
SetVariableInteger(int index,bool integer)680 void GurobiInterface::SetVariableInteger(int index, bool integer) {
681   InvalidateSolutionSynchronization();
682   if (!had_nonincremental_change_ && variable_is_extracted(index)) {
683     char type_var;
684     if (integer) {
685       type_var = GRB_INTEGER;
686     } else {
687       type_var = GRB_CONTINUOUS;
688     }
689     SetCharAttrElement(GRB_CHAR_ATTR_VTYPE, mp_var_to_gurobi_var_.at(index),
690                        type_var);
691   } else {
692     sync_status_ = MUST_RELOAD;
693   }
694 }
695 
SetConstraintBounds(int index,double lb,double ub)696 void GurobiInterface::SetConstraintBounds(int index, double lb, double ub) {
697   sync_status_ = MUST_RELOAD;
698   if (constraint_is_extracted(index)) {
699     had_nonincremental_change_ = true;
700   }
701   // TODO(user): this is nontrivial to make incremental:
702   //   1. Make sure it is a linear constraint (not an indicator or indicator
703   //      range constraint).
704   //   2. Check if the sense of the constraint changes. If it was previously a
705   //      range constraint, we can do nothing, and if it becomes a range
706   //      constraint, we can do nothing. We could support range constraints if
707   //      we tracked the auxiliary variable that is added with range
708   //      constraints.
709 }
710 
AddRowConstraint(MPConstraint * const ct)711 void GurobiInterface::AddRowConstraint(MPConstraint* const ct) {
712   sync_status_ = MUST_RELOAD;
713 }
714 
AddIndicatorConstraint(MPConstraint * const ct)715 bool GurobiInterface::AddIndicatorConstraint(MPConstraint* const ct) {
716   had_nonincremental_change_ = true;
717   sync_status_ = MUST_RELOAD;
718   return !IsContinuous();
719 }
720 
AddVariable(MPVariable * const var)721 void GurobiInterface::AddVariable(MPVariable* const var) {
722   sync_status_ = MUST_RELOAD;
723 }
724 
SetCoefficient(MPConstraint * const constraint,const MPVariable * const variable,double new_value,double old_value)725 void GurobiInterface::SetCoefficient(MPConstraint* const constraint,
726                                      const MPVariable* const variable,
727                                      double new_value, double old_value) {
728   InvalidateSolutionSynchronization();
729   if (!had_nonincremental_change_ && variable_is_extracted(variable->index()) &&
730       constraint_is_extracted(constraint->index())) {
731     // Cannot be const, GRBchgcoeffs needs non-const pointer.
732     int grb_var = mp_var_to_gurobi_var_.at(variable->index());
733     int grb_cons = mp_cons_to_gurobi_linear_cons_.at(constraint->index());
734     if (grb_cons < 0) {
735       had_nonincremental_change_ = true;
736       sync_status_ = MUST_RELOAD;
737     } else {
738       // TODO(user): investigate if this has bad performance.
739       CheckedGurobiCall(
740           GRBchgcoeffs(model_, 1, &grb_cons, &grb_var, &new_value));
741     }
742   } else {
743     sync_status_ = MUST_RELOAD;
744   }
745 }
746 
ClearConstraint(MPConstraint * const constraint)747 void GurobiInterface::ClearConstraint(MPConstraint* const constraint) {
748   had_nonincremental_change_ = true;
749   sync_status_ = MUST_RELOAD;
750   // TODO(user): this is difficult to make incremental, like
751   //  SetConstraintBounds(), because of the auxiliary Gurobi variables that
752   //  range constraints introduce.
753 }
754 
SetObjectiveCoefficient(const MPVariable * const variable,double coefficient)755 void GurobiInterface::SetObjectiveCoefficient(const MPVariable* const variable,
756                                               double coefficient) {
757   InvalidateSolutionSynchronization();
758   if (!had_nonincremental_change_ && variable_is_extracted(variable->index())) {
759     SetDoubleAttrElement(GRB_DBL_ATTR_OBJ,
760                          mp_var_to_gurobi_var_.at(variable->index()),
761                          coefficient);
762   } else {
763     sync_status_ = MUST_RELOAD;
764   }
765 }
766 
SetObjectiveOffset(double value)767 void GurobiInterface::SetObjectiveOffset(double value) {
768   InvalidateSolutionSynchronization();
769   if (!had_nonincremental_change_) {
770     SetDoubleAttr(GRB_DBL_ATTR_OBJCON, value);
771   } else {
772     sync_status_ = MUST_RELOAD;
773   }
774 }
775 
ClearObjective()776 void GurobiInterface::ClearObjective() {
777   InvalidateSolutionSynchronization();
778   if (!had_nonincremental_change_) {
779     SetObjectiveOffset(0.0);
780     for (const auto& entry : solver_->objective_->coefficients_) {
781       SetObjectiveCoefficient(entry.first, 0.0);
782     }
783   } else {
784     sync_status_ = MUST_RELOAD;
785   }
786 }
787 
BranchingPriorityChangedForVariable(int var_index)788 void GurobiInterface::BranchingPriorityChangedForVariable(int var_index) {
789   update_branching_priorities_ = true;
790 }
791 
792 // ------ Query statistics on the solution and the solve ------
793 
iterations() const794 int64_t GurobiInterface::iterations() const {
795   double iter;
796   if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations;
797   CheckedGurobiCall(GRBgetdblattr(model_, GRB_DBL_ATTR_ITERCOUNT, &iter));
798   return static_cast<int64_t>(iter);
799 }
800 
nodes() const801 int64_t GurobiInterface::nodes() const {
802   if (mip_) {
803     if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes;
804     return static_cast<int64_t>(GetDoubleAttr(GRB_DBL_ATTR_NODECOUNT));
805   } else {
806     LOG(DFATAL) << "Number of nodes only available for discrete problems.";
807     return kUnknownNumberOfNodes;
808   }
809 }
810 
TransformGRBVarBasisStatus(int gurobi_basis_status) const811 MPSolver::BasisStatus GurobiInterface::TransformGRBVarBasisStatus(
812     int gurobi_basis_status) const {
813   switch (gurobi_basis_status) {
814     case GRB_BASIC:
815       return MPSolver::BASIC;
816     case GRB_NONBASIC_LOWER:
817       return MPSolver::AT_LOWER_BOUND;
818     case GRB_NONBASIC_UPPER:
819       return MPSolver::AT_UPPER_BOUND;
820     case GRB_SUPERBASIC:
821       return MPSolver::FREE;
822     default:
823       LOG(DFATAL) << "Unknown GRB basis status.";
824       return MPSolver::FREE;
825   }
826 }
827 
TransformGRBConstraintBasisStatus(int gurobi_basis_status,int constraint_index) const828 MPSolver::BasisStatus GurobiInterface::TransformGRBConstraintBasisStatus(
829     int gurobi_basis_status, int constraint_index) const {
830   const int grb_index = mp_cons_to_gurobi_linear_cons_.at(constraint_index);
831   if (grb_index < 0) {
832     LOG(DFATAL) << "Basis status not available for nonlinear constraints.";
833     return MPSolver::FREE;
834   }
835   switch (gurobi_basis_status) {
836     case GRB_BASIC:
837       return MPSolver::BASIC;
838     default: {
839       // Non basic.
840       double tolerance = 0.0;
841       CheckedGurobiCall(GRBgetdblparam(GRBgetenv(model_),
842                                        GRB_DBL_PAR_FEASIBILITYTOL, &tolerance));
843       const double slack = GetDoubleAttrElement(GRB_DBL_ATTR_SLACK, grb_index);
844       const char sense = GetCharAttrElement(GRB_CHAR_ATTR_SENSE, grb_index);
845       VLOG(4) << "constraint " << constraint_index << " , slack = " << slack
846               << " , sense = " << sense;
847       if (fabs(slack) <= tolerance) {
848         switch (sense) {
849           case GRB_EQUAL:
850           case GRB_LESS_EQUAL:
851             return MPSolver::AT_UPPER_BOUND;
852           case GRB_GREATER_EQUAL:
853             return MPSolver::AT_LOWER_BOUND;
854           default:
855             return MPSolver::FREE;
856         }
857       } else {
858         return MPSolver::FREE;
859       }
860     }
861   }
862 }
863 
864 // Returns the basis status of a row.
row_status(int constraint_index) const865 MPSolver::BasisStatus GurobiInterface::row_status(int constraint_index) const {
866   const int optim_status = GetIntAttr(GRB_INT_ATTR_STATUS);
867   if (optim_status != GRB_OPTIMAL && optim_status != GRB_SUBOPTIMAL) {
868     LOG(DFATAL) << "Basis status only available after a solution has "
869                 << "been found.";
870     return MPSolver::FREE;
871   }
872   if (mip_) {
873     LOG(DFATAL) << "Basis status only available for continuous problems.";
874     return MPSolver::FREE;
875   }
876   const int grb_index = mp_cons_to_gurobi_linear_cons_.at(constraint_index);
877   if (grb_index < 0) {
878     LOG(DFATAL) << "Basis status not available for nonlinear constraints.";
879     return MPSolver::FREE;
880   }
881   const int gurobi_basis_status =
882       GetIntAttrElement(GRB_INT_ATTR_CBASIS, grb_index);
883   return TransformGRBConstraintBasisStatus(gurobi_basis_status,
884                                            constraint_index);
885 }
886 
887 // Returns the basis status of a column.
column_status(int variable_index) const888 MPSolver::BasisStatus GurobiInterface::column_status(int variable_index) const {
889   const int optim_status = GetIntAttr(GRB_INT_ATTR_STATUS);
890   if (optim_status != GRB_OPTIMAL && optim_status != GRB_SUBOPTIMAL) {
891     LOG(DFATAL) << "Basis status only available after a solution has "
892                 << "been found.";
893     return MPSolver::FREE;
894   }
895   if (mip_) {
896     LOG(DFATAL) << "Basis status only available for continuous problems.";
897     return MPSolver::FREE;
898   }
899   const int grb_index = mp_var_to_gurobi_var_.at(variable_index);
900   const int gurobi_basis_status =
901       GetIntAttrElement(GRB_INT_ATTR_VBASIS, grb_index);
902   return TransformGRBVarBasisStatus(gurobi_basis_status);
903 }
904 
905 // Extracts new variables.
ExtractNewVariables()906 void GurobiInterface::ExtractNewVariables() {
907   const int total_num_vars = solver_->variables_.size();
908   if (total_num_vars > last_variable_index_) {
909     // Define new variables.
910     for (int j = last_variable_index_; j < total_num_vars; ++j) {
911       const MPVariable* const var = solver_->variables_.at(j);
912       set_variable_as_extracted(var->index(), true);
913       CheckedGurobiCall(GRBaddvar(
914           model_, 0,  // numnz
915           nullptr,    // vind
916           nullptr,    // vval
917           solver_->objective_->GetCoefficient(var), var->lb(), var->ub(),
918           var->integer() && mip_ ? GRB_INTEGER : GRB_CONTINUOUS,
919           var->name().empty() ? nullptr : var->name().c_str()));
920       mp_var_to_gurobi_var_.push_back(num_gurobi_vars_++);
921     }
922     CheckedGurobiCall(GRBupdatemodel(model_));
923     // Add new variables to existing constraints.
924     std::vector<int> grb_cons_ind;
925     std::vector<int> grb_var_ind;
926     std::vector<double> coef;
927     for (int i = 0; i < last_constraint_index_; ++i) {
928       // If there was a nonincremental change/the model is not incremental (e.g.
929       // there is an indicator constraint), we should never enter this loop, as
930       // last_variable_index_ will be reset to zero before ExtractNewVariables()
931       // is called.
932       MPConstraint* const ct = solver_->constraints_[i];
933       const int grb_ct_idx = mp_cons_to_gurobi_linear_cons_.at(ct->index());
934       DCHECK_GE(grb_ct_idx, 0);
935       DCHECK(ct->indicator_variable() == nullptr);
936       for (const auto& entry : ct->coefficients_) {
937         const int var_index = entry.first->index();
938         DCHECK(variable_is_extracted(var_index));
939 
940         if (var_index >= last_variable_index_) {
941           grb_cons_ind.push_back(grb_ct_idx);
942           grb_var_ind.push_back(mp_var_to_gurobi_var_.at(var_index));
943           coef.push_back(entry.second);
944         }
945       }
946     }
947     if (!grb_cons_ind.empty()) {
948       CheckedGurobiCall(GRBchgcoeffs(model_, grb_cons_ind.size(),
949                                      grb_cons_ind.data(), grb_var_ind.data(),
950                                      coef.data()));
951     }
952   }
953   CheckedGurobiCall(GRBupdatemodel(model_));
954   DCHECK_EQ(GetIntAttr(GRB_INT_ATTR_NUMVARS), num_gurobi_vars_);
955 }
956 
ExtractNewConstraints()957 void GurobiInterface::ExtractNewConstraints() {
958   int total_num_rows = solver_->constraints_.size();
959   if (last_constraint_index_ < total_num_rows) {
960     // Add each new constraint.
961     for (int row = last_constraint_index_; row < total_num_rows; ++row) {
962       MPConstraint* const ct = solver_->constraints_[row];
963       set_constraint_as_extracted(row, true);
964       const int size = ct->coefficients_.size();
965       std::vector<int> grb_vars;
966       std::vector<double> coefs;
967       grb_vars.reserve(size);
968       coefs.reserve(size);
969       for (const auto& entry : ct->coefficients_) {
970         const int var_index = entry.first->index();
971         CHECK(variable_is_extracted(var_index));
972         grb_vars.push_back(mp_var_to_gurobi_var_.at(var_index));
973         coefs.push_back(entry.second);
974       }
975       char* const name =
976           ct->name().empty() ? nullptr : const_cast<char*>(ct->name().c_str());
977       if (ct->indicator_variable() != nullptr) {
978         const int grb_ind_var =
979             mp_var_to_gurobi_var_.at(ct->indicator_variable()->index());
980         if (ct->lb() > -std::numeric_limits<double>::infinity()) {
981           CheckedGurobiCall(GRBaddgenconstrIndicator(
982               model_, name, grb_ind_var, ct->indicator_value(), size,
983               grb_vars.data(), coefs.data(),
984               ct->ub() == ct->lb() ? GRB_EQUAL : GRB_GREATER_EQUAL, ct->lb()));
985         }
986         if (ct->ub() < std::numeric_limits<double>::infinity() &&
987             ct->lb() != ct->ub()) {
988           CheckedGurobiCall(GRBaddgenconstrIndicator(
989               model_, name, grb_ind_var, ct->indicator_value(), size,
990               grb_vars.data(), coefs.data(), GRB_LESS_EQUAL, ct->ub()));
991         }
992         mp_cons_to_gurobi_linear_cons_.push_back(-1);
993       } else {
994         // Using GRBaddrangeconstr for constraints that don't require it adds
995         // a slack which is not always removed by presolve.
996         if (ct->lb() == ct->ub()) {
997           CheckedGurobiCall(GRBaddconstr(model_, size, grb_vars.data(),
998                                          coefs.data(), GRB_EQUAL, ct->lb(),
999                                          name));
1000         } else if (ct->lb() == -std::numeric_limits<double>::infinity()) {
1001           CheckedGurobiCall(GRBaddconstr(model_, size, grb_vars.data(),
1002                                          coefs.data(), GRB_LESS_EQUAL, ct->ub(),
1003                                          name));
1004         } else if (ct->ub() == std::numeric_limits<double>::infinity()) {
1005           CheckedGurobiCall(GRBaddconstr(model_, size, grb_vars.data(),
1006                                          coefs.data(), GRB_GREATER_EQUAL,
1007                                          ct->lb(), name));
1008         } else {
1009           CheckedGurobiCall(GRBaddrangeconstr(model_, size, grb_vars.data(),
1010                                               coefs.data(), ct->lb(), ct->ub(),
1011                                               name));
1012           // NOTE(user): range constraints implicitly add an extra variable
1013           // to the model.
1014           num_gurobi_vars_++;
1015         }
1016         mp_cons_to_gurobi_linear_cons_.push_back(num_gurobi_linear_cons_++);
1017       }
1018     }
1019   }
1020   CheckedGurobiCall(GRBupdatemodel(model_));
1021   DCHECK_EQ(GetIntAttr(GRB_INT_ATTR_NUMCONSTRS), num_gurobi_linear_cons_);
1022 }
1023 
ExtractObjective()1024 void GurobiInterface::ExtractObjective() {
1025   SetIntAttr(GRB_INT_ATTR_MODELSENSE, maximize_ ? GRB_MAXIMIZE : GRB_MINIMIZE);
1026   SetDoubleAttr(GRB_DBL_ATTR_OBJCON, solver_->Objective().offset());
1027 }
1028 
1029 // ------ Parameters  -----
1030 
SetParameters(const MPSolverParameters & param)1031 void GurobiInterface::SetParameters(const MPSolverParameters& param) {
1032   SetCommonParameters(param);
1033   if (mip_) {
1034     SetMIPParameters(param);
1035   }
1036 }
1037 
SetSolverSpecificParametersAsString(const std::string & parameters)1038 bool GurobiInterface::SetSolverSpecificParametersAsString(
1039     const std::string& parameters) {
1040   return SetSolverSpecificParameters(parameters, GRBgetenv(model_)).ok();
1041 }
1042 
SetRelativeMipGap(double value)1043 void GurobiInterface::SetRelativeMipGap(double value) {
1044   if (mip_) {
1045     CheckedGurobiCall(
1046         GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_MIPGAP, value));
1047   } else {
1048     LOG(WARNING) << "The relative MIP gap is only available "
1049                  << "for discrete problems.";
1050   }
1051 }
1052 
1053 // Gurobi has two different types of primal tolerance (feasibility tolerance):
1054 // constraint and integrality. We need to set them both.
1055 // See:
1056 // http://www.gurobi.com/documentation/6.0/refman/feasibilitytol.html
1057 // and
1058 // http://www.gurobi.com/documentation/6.0/refman/intfeastol.html
SetPrimalTolerance(double value)1059 void GurobiInterface::SetPrimalTolerance(double value) {
1060   CheckedGurobiCall(
1061       GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_FEASIBILITYTOL, value));
1062   CheckedGurobiCall(
1063       GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_INTFEASTOL, value));
1064 }
1065 
1066 // As opposed to primal (feasibility) tolerance, the dual (optimality) tolerance
1067 // applies only to the reduced costs in the improving direction.
1068 // See:
1069 // http://www.gurobi.com/documentation/6.0/refman/optimalitytol.html
SetDualTolerance(double value)1070 void GurobiInterface::SetDualTolerance(double value) {
1071   CheckedGurobiCall(
1072       GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_OPTIMALITYTOL, value));
1073 }
1074 
SetPresolveMode(int value)1075 void GurobiInterface::SetPresolveMode(int value) {
1076   switch (value) {
1077     case MPSolverParameters::PRESOLVE_OFF: {
1078       CheckedGurobiCall(
1079           GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_PRESOLVE, false));
1080       break;
1081     }
1082     case MPSolverParameters::PRESOLVE_ON: {
1083       CheckedGurobiCall(
1084           GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_PRESOLVE, true));
1085       break;
1086     }
1087     default: {
1088       SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value);
1089     }
1090   }
1091 }
1092 
1093 // Sets the scaling mode.
SetScalingMode(int value)1094 void GurobiInterface::SetScalingMode(int value) {
1095   switch (value) {
1096     case MPSolverParameters::SCALING_OFF:
1097       CheckedGurobiCall(
1098           GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_SCALEFLAG, false));
1099       break;
1100     case MPSolverParameters::SCALING_ON:
1101       CheckedGurobiCall(
1102           GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_SCALEFLAG, true));
1103       CheckedGurobiCall(
1104           GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_OBJSCALE, 0.0));
1105       break;
1106     default:
1107       // Leave the parameters untouched.
1108       break;
1109   }
1110 }
1111 
1112 // Sets the LP algorithm : primal, dual or barrier. Note that GRB
1113 // offers automatic selection
SetLpAlgorithm(int value)1114 void GurobiInterface::SetLpAlgorithm(int value) {
1115   switch (value) {
1116     case MPSolverParameters::DUAL:
1117       CheckedGurobiCall(GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_METHOD,
1118                                        GRB_METHOD_DUAL));
1119       break;
1120     case MPSolverParameters::PRIMAL:
1121       CheckedGurobiCall(GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_METHOD,
1122                                        GRB_METHOD_PRIMAL));
1123       break;
1124     case MPSolverParameters::BARRIER:
1125       CheckedGurobiCall(GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_METHOD,
1126                                        GRB_METHOD_BARRIER));
1127       break;
1128     default:
1129       SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM,
1130                                         value);
1131   }
1132 }
1133 
SolutionCount() const1134 int GurobiInterface::SolutionCount() const {
1135   return GetIntAttr(GRB_INT_ATTR_SOLCOUNT);
1136 }
1137 
ModelIsNonincremental() const1138 bool GurobiInterface::ModelIsNonincremental() const {
1139   for (const MPConstraint* c : solver_->constraints()) {
1140     if (c->indicator_variable() != nullptr) {
1141       return true;
1142     }
1143   }
1144   return false;
1145 }
1146 
Solve(const MPSolverParameters & param)1147 MPSolver::ResultStatus GurobiInterface::Solve(const MPSolverParameters& param) {
1148   WallTimer timer;
1149   timer.Start();
1150 
1151   if (param.GetIntegerParam(MPSolverParameters::INCREMENTALITY) ==
1152           MPSolverParameters::INCREMENTALITY_OFF ||
1153       ModelIsNonincremental() || had_nonincremental_change_) {
1154     Reset();
1155   }
1156 
1157   // Set log level.
1158   CheckedGurobiCall(
1159       GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_OUTPUTFLAG, !quiet_));
1160 
1161   ExtractModel();
1162   // Sync solver.
1163   CheckedGurobiCall(GRBupdatemodel(model_));
1164   VLOG(1) << absl::StrFormat("Model built in %s.",
1165                              absl::FormatDuration(timer.GetDuration()));
1166 
1167   // Set solution hints if any.
1168   for (const std::pair<const MPVariable*, double>& p :
1169        solver_->solution_hint_) {
1170     SetDoubleAttrElement(GRB_DBL_ATTR_START,
1171                          mp_var_to_gurobi_var_.at(p.first->index()), p.second);
1172   }
1173 
1174   // Pass branching priority annotations if at least one has been updated.
1175   if (update_branching_priorities_) {
1176     for (const MPVariable* var : solver_->variables_) {
1177       SetIntAttrElement(GRB_INT_ATTR_BRANCHPRIORITY,
1178                         mp_var_to_gurobi_var_.at(var->index()),
1179                         var->branching_priority());
1180     }
1181     update_branching_priorities_ = false;
1182   }
1183 
1184   // Time limit.
1185   if (solver_->time_limit() != 0) {
1186     VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
1187     CheckedGurobiCall(GRBsetdblparam(GRBgetenv(model_), GRB_DBL_PAR_TIMELIMIT,
1188                                      solver_->time_limit_in_secs()));
1189   }
1190 
1191   // We first set our internal MPSolverParameters from 'param' and then set
1192   // any user-specified internal solver parameters via
1193   // solver_specific_parameter_string_.
1194   // Default MPSolverParameters can override custom parameters (for example for
1195   // presolving) and therefore we apply MPSolverParameters first.
1196   SetParameters(param);
1197   solver_->SetSolverSpecificParametersAsString(
1198       solver_->solver_specific_parameter_string_);
1199 
1200   std::unique_ptr<GurobiMPCallbackContext> gurobi_context;
1201   MPCallbackWithGurobiContext mp_callback_with_context;
1202   int gurobi_precrush = 0;
1203   int gurobi_lazy_constraint = 0;
1204   if (callback_ == nullptr) {
1205     CheckedGurobiCall(GRBsetcallbackfunc(model_, nullptr, nullptr));
1206   } else {
1207     gurobi_context = absl::make_unique<GurobiMPCallbackContext>(
1208         env_, &mp_var_to_gurobi_var_, num_gurobi_vars_,
1209         callback_->might_add_cuts(), callback_->might_add_lazy_constraints());
1210     mp_callback_with_context.context = gurobi_context.get();
1211     mp_callback_with_context.callback = callback_;
1212     CheckedGurobiCall(GRBsetcallbackfunc(
1213         model_, CallbackImpl, static_cast<void*>(&mp_callback_with_context)));
1214     gurobi_precrush = callback_->might_add_cuts();
1215     gurobi_lazy_constraint = callback_->might_add_lazy_constraints();
1216   }
1217   CheckedGurobiCall(
1218       GRBsetintparam(GRBgetenv(model_), GRB_INT_PAR_PRECRUSH, gurobi_precrush));
1219   CheckedGurobiCall(GRBsetintparam(
1220       GRBgetenv(model_), GRB_INT_PAR_LAZYCONSTRAINTS, gurobi_lazy_constraint));
1221 
1222   // Solve
1223   timer.Restart();
1224   const int status = GRBoptimize(model_);
1225 
1226   if (status) {
1227     VLOG(1) << "Failed to optimize MIP." << GRBgeterrormsg(env_);
1228   } else {
1229     VLOG(1) << absl::StrFormat("Solved in %s.",
1230                                absl::FormatDuration(timer.GetDuration()));
1231   }
1232 
1233   // Get the status.
1234   const int optimization_status = GetIntAttr(GRB_INT_ATTR_STATUS);
1235   VLOG(1) << absl::StrFormat("Solution status %d.\n", optimization_status);
1236   const int solution_count = SolutionCount();
1237 
1238   switch (optimization_status) {
1239     case GRB_OPTIMAL:
1240       result_status_ = MPSolver::OPTIMAL;
1241       break;
1242     case GRB_INFEASIBLE:
1243       result_status_ = MPSolver::INFEASIBLE;
1244       break;
1245     case GRB_UNBOUNDED:
1246       result_status_ = MPSolver::UNBOUNDED;
1247       break;
1248     case GRB_INF_OR_UNBD:
1249       // TODO(user): We could introduce our own "infeasible or
1250       // unbounded" status.
1251       result_status_ = MPSolver::INFEASIBLE;
1252       break;
1253     default: {
1254       if (solution_count > 0) {
1255         result_status_ = MPSolver::FEASIBLE;
1256       } else {
1257         result_status_ = MPSolver::NOT_SOLVED;
1258       }
1259       break;
1260     }
1261   }
1262 
1263   if (IsMIP() && (result_status_ != MPSolver::UNBOUNDED &&
1264                   result_status_ != MPSolver::INFEASIBLE)) {
1265     const int error =
1266         GRBgetdblattr(model_, GRB_DBL_ATTR_OBJBOUND, &best_objective_bound_);
1267     LOG_IF(WARNING, error != 0)
1268         << "Best objective bound is not available, error=" << error
1269         << ", message=" << GRBgeterrormsg(env_);
1270     VLOG(1) << "best bound = " << best_objective_bound_;
1271   }
1272 
1273   if (solution_count > 0 && (result_status_ == MPSolver::FEASIBLE ||
1274                              result_status_ == MPSolver::OPTIMAL)) {
1275     current_solution_index_ = 0;
1276     // Get the results.
1277     objective_value_ = GetDoubleAttr(GRB_DBL_ATTR_OBJVAL);
1278     VLOG(1) << "objective = " << objective_value_;
1279 
1280     {
1281       const std::vector<double> grb_variable_values =
1282           GetDoubleAttrArray(GRB_DBL_ATTR_X, num_gurobi_vars_);
1283       for (int i = 0; i < solver_->variables_.size(); ++i) {
1284         MPVariable* const var = solver_->variables_[i];
1285         const double val = grb_variable_values.at(mp_var_to_gurobi_var_.at(i));
1286         var->set_solution_value(val);
1287         VLOG(3) << var->name() << ", value = " << val;
1288       }
1289     }
1290     if (!mip_) {
1291       {
1292         const std::vector<double> grb_reduced_costs =
1293             GetDoubleAttrArray(GRB_DBL_ATTR_RC, num_gurobi_vars_);
1294         for (int i = 0; i < solver_->variables_.size(); ++i) {
1295           MPVariable* const var = solver_->variables_[i];
1296           const double rc = grb_reduced_costs.at(mp_var_to_gurobi_var_.at(i));
1297           var->set_reduced_cost(rc);
1298           VLOG(4) << var->name() << ", reduced cost = " << rc;
1299         }
1300       }
1301 
1302       {
1303         std::vector<double> grb_dual_values =
1304             GetDoubleAttrArray(GRB_DBL_ATTR_PI, num_gurobi_linear_cons_);
1305         for (int i = 0; i < solver_->constraints_.size(); ++i) {
1306           MPConstraint* const ct = solver_->constraints_[i];
1307           const double dual_value =
1308               grb_dual_values.at(mp_cons_to_gurobi_linear_cons_.at(i));
1309           ct->set_dual_value(dual_value);
1310           VLOG(4) << "row " << ct->index() << ", dual value = " << dual_value;
1311         }
1312       }
1313     }
1314   }
1315 
1316   sync_status_ = SOLUTION_SYNCHRONIZED;
1317   GRBresetparams(GRBgetenv(model_));
1318   return result_status_;
1319 }
1320 
DirectlySolveProto(const MPModelRequest & request,std::atomic<bool> * interrupt)1321 absl::optional<MPSolutionResponse> GurobiInterface::DirectlySolveProto(
1322     const MPModelRequest& request, std::atomic<bool>* interrupt) {
1323   // Interruption via atomic<bool> is not directly supported by Gurobi.
1324   if (interrupt != nullptr) return absl::nullopt;
1325 
1326   // Here we reuse the Gurobi environment to support single-use license that
1327   // forbids creating a second environment if one already exists.
1328   const auto status_or = GurobiSolveProto(request, env_);
1329   if (status_or.ok()) return status_or.value();
1330   // Special case: if something is not implemented yet, fall back to solving
1331   // through MPSolver.
1332   if (absl::IsUnimplemented(status_or.status())) return absl::nullopt;
1333 
1334   if (request.enable_internal_solver_output()) {
1335     LOG(INFO) << "Invalid Gurobi status: " << status_or.status();
1336   }
1337   MPSolutionResponse response;
1338   response.set_status(MPSOLVER_NOT_SOLVED);
1339   response.set_status_str(status_or.status().ToString());
1340   return response;
1341 }
1342 
NextSolution()1343 bool GurobiInterface::NextSolution() {
1344   // Next solution only supported for MIP
1345   if (!mip_) return false;
1346 
1347   // Make sure we have successfully solved the problem and not modified it.
1348   if (!CheckSolutionIsSynchronizedAndExists()) {
1349     return false;
1350   }
1351   // Check if we are out of solutions.
1352   if (current_solution_index_ + 1 >= SolutionCount()) {
1353     return false;
1354   }
1355   current_solution_index_++;
1356 
1357   CheckedGurobiCall(GRBsetintparam(
1358       GRBgetenv(model_), GRB_INT_PAR_SOLUTIONNUMBER, current_solution_index_));
1359 
1360   objective_value_ = GetDoubleAttr(GRB_DBL_ATTR_POOLOBJVAL);
1361   const std::vector<double> grb_variable_values =
1362       GetDoubleAttrArray(GRB_DBL_ATTR_XN, num_gurobi_vars_);
1363 
1364   for (int i = 0; i < solver_->variables_.size(); ++i) {
1365     MPVariable* const var = solver_->variables_[i];
1366     var->set_solution_value(
1367         grb_variable_values.at(mp_var_to_gurobi_var_.at(i)));
1368   }
1369   // TODO(user): This reset may not be necessary, investigate.
1370   GRBresetparams(GRBgetenv(model_));
1371   return true;
1372 }
1373 
Write(const std::string & filename)1374 void GurobiInterface::Write(const std::string& filename) {
1375   if (sync_status_ == MUST_RELOAD) {
1376     Reset();
1377   }
1378   ExtractModel();
1379   // Sync solver.
1380   CheckedGurobiCall(GRBupdatemodel(model_));
1381   VLOG(1) << "Writing Gurobi model file \"" << filename << "\".";
1382   const int status = GRBwrite(model_, filename.c_str());
1383   if (status) {
1384     LOG(WARNING) << "Failed to write MIP." << GRBgeterrormsg(env_);
1385   }
1386 }
1387 
BuildGurobiInterface(bool mip,MPSolver * const solver)1388 MPSolverInterface* BuildGurobiInterface(bool mip, MPSolver* const solver) {
1389   return new GurobiInterface(solver, mip);
1390 }
1391 
SetCallback(MPCallback * mp_callback)1392 void GurobiInterface::SetCallback(MPCallback* mp_callback) {
1393   callback_ = mp_callback;
1394 }
1395 
1396 }  // namespace operations_research
1397