OR-Tools  8.2
scip_proto_solver.cc
Go to the documentation of this file.
1 // Copyright 2010-2018 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 #if defined(USE_SCIP)
15 
17 
18 #include <cmath>
19 #include <limits>
20 #include <memory>
21 #include <numeric>
22 #include <set>
23 #include <string>
24 #include <vector>
25 
26 #include "absl/status/status.h"
27 #include "absl/status/statusor.h"
28 #include "absl/strings/ascii.h"
29 #include "absl/strings/numbers.h"
30 #include "absl/strings/str_cat.h"
31 #include "absl/strings/str_format.h"
32 #include "absl/strings/str_split.h"
33 #include "ortools/base/cleanup.h"
37 #include "ortools/linear_solver/linear_solver.pb.h"
41 #include "scip/cons_disjunction.h"
42 #include "scip/cons_linear.h"
43 #include "scip/pub_var.h"
44 #include "scip/scip.h"
45 #include "scip/scip_param.h"
46 #include "scip/scip_prob.h"
47 #include "scip/scip_var.h"
48 #include "scip/scipdefplugins.h"
49 #include "scip/set.h"
50 #include "scip/struct_paramset.h"
51 #include "scip/type_cons.h"
52 #include "scip/type_paramset.h"
53 #include "scip/type_var.h"
54 
55 ABSL_FLAG(std::string, scip_proto_solver_output_cip_file, "",
56  "If given, saves the generated CIP file here. Useful for "
57  "reporting bugs to SCIP.");
58 
59 namespace operations_research {
60 
61 namespace {
62 // This function will create a new constraint if the indicator constraint has
63 // both a lower bound and an upper bound.
64 absl::Status AddIndicatorConstraint(const MPGeneralConstraintProto& gen_cst,
65  SCIP* scip, SCIP_CONS** scip_cst,
66  std::vector<SCIP_VAR*>* scip_variables,
67  std::vector<SCIP_CONS*>* scip_constraints,
68  std::vector<SCIP_VAR*>* tmp_variables,
69  std::vector<double>* tmp_coefficients) {
70  CHECK(scip != nullptr);
71  CHECK(scip_cst != nullptr);
72  CHECK(scip_variables != nullptr);
73  CHECK(scip_constraints != nullptr);
74  CHECK(tmp_variables != nullptr);
75  CHECK(tmp_coefficients != nullptr);
76  CHECK(gen_cst.has_indicator_constraint());
77  constexpr double kInfinity = std::numeric_limits<double>::infinity();
78 
79  const auto& ind = gen_cst.indicator_constraint();
80  if (!ind.has_constraint()) return absl::OkStatus();
81 
82  const MPConstraintProto& constraint = ind.constraint();
83  const int size = constraint.var_index_size();
84  tmp_variables->resize(size, nullptr);
85  tmp_coefficients->resize(size, 0);
86  for (int i = 0; i < size; ++i) {
87  (*tmp_variables)[i] = (*scip_variables)[constraint.var_index(i)];
88  (*tmp_coefficients)[i] = constraint.coefficient(i);
89  }
90 
91  SCIP_VAR* ind_var = (*scip_variables)[ind.var_index()];
92  if (ind.var_value() == 0) {
94  SCIPgetNegatedVar(scip, (*scip_variables)[ind.var_index()], &ind_var));
95  }
96 
97  if (ind.constraint().upper_bound() < kInfinity) {
98  RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator(
99  scip, scip_cst, gen_cst.name().c_str(), ind_var, size,
100  tmp_variables->data(), tmp_coefficients->data(),
101  ind.constraint().upper_bound(),
102  /*initial=*/!ind.constraint().is_lazy(),
103  /*separate=*/true,
104  /*enforce=*/true,
105  /*check=*/true,
106  /*propagate=*/true,
107  /*local=*/false,
108  /*dynamic=*/false,
109  /*removable=*/ind.constraint().is_lazy(),
110  /*stickingatnode=*/false));
111  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
112  scip_constraints->push_back(nullptr);
113  scip_cst = &scip_constraints->back();
114  }
115  if (ind.constraint().lower_bound() > -kInfinity) {
116  for (int i = 0; i < size; ++i) {
117  (*tmp_coefficients)[i] *= -1;
118  }
119  RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator(
120  scip, scip_cst, gen_cst.name().c_str(), ind_var, size,
121  tmp_variables->data(), tmp_coefficients->data(),
122  -ind.constraint().lower_bound(),
123  /*initial=*/!ind.constraint().is_lazy(),
124  /*separate=*/true,
125  /*enforce=*/true,
126  /*check=*/true,
127  /*propagate=*/true,
128  /*local=*/false,
129  /*dynamic=*/false,
130  /*removable=*/ind.constraint().is_lazy(),
131  /*stickingatnode=*/false));
132  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
133  }
134 
135  return absl::OkStatus();
136 }
137 
138 absl::Status AddSosConstraint(const MPGeneralConstraintProto& gen_cst,
139  const std::vector<SCIP_VAR*>& scip_variables,
140  SCIP* scip, SCIP_CONS** scip_cst,
141  std::vector<SCIP_VAR*>* tmp_variables,
142  std::vector<double>* tmp_weights) {
143  CHECK(scip != nullptr);
144  CHECK(scip_cst != nullptr);
145  CHECK(tmp_variables != nullptr);
146  CHECK(tmp_weights != nullptr);
147 
148  CHECK(gen_cst.has_sos_constraint());
149  const MPSosConstraint& sos_cst = gen_cst.sos_constraint();
150 
151  // SOS constraints of type N indicate at most N variables are non-zero.
152  // Constraints with N variables or less are valid, but useless. They also
153  // crash SCIP, so we skip them.
154  if (sos_cst.var_index_size() <= 1) return absl::OkStatus();
155  if (sos_cst.type() == MPSosConstraint::SOS2 &&
156  sos_cst.var_index_size() <= 2) {
157  return absl::OkStatus();
158  }
159 
160  tmp_variables->resize(sos_cst.var_index_size(), nullptr);
161  for (int v = 0; v < sos_cst.var_index_size(); ++v) {
162  (*tmp_variables)[v] = scip_variables[sos_cst.var_index(v)];
163  }
164  tmp_weights->resize(sos_cst.var_index_size(), 0);
165  if (sos_cst.weight_size() == sos_cst.var_index_size()) {
166  for (int w = 0; w < sos_cst.weight_size(); ++w) {
167  (*tmp_weights)[w] = sos_cst.weight(w);
168  }
169  } else {
170  // In theory, SCIP should accept empty weight arrays and use natural
171  // ordering, but in practice, this crashes their code.
172  std::iota(tmp_weights->begin(), tmp_weights->end(), 1);
173  }
174  switch (sos_cst.type()) {
175  case MPSosConstraint::SOS1_DEFAULT:
177  SCIPcreateConsBasicSOS1(scip,
178  /*cons=*/scip_cst,
179  /*name=*/gen_cst.name().c_str(),
180  /*nvars=*/sos_cst.var_index_size(),
181  /*vars=*/tmp_variables->data(),
182  /*weights=*/tmp_weights->data()));
183  break;
184  case MPSosConstraint::SOS2:
186  SCIPcreateConsBasicSOS2(scip,
187  /*cons=*/scip_cst,
188  /*name=*/gen_cst.name().c_str(),
189  /*nvars=*/sos_cst.var_index_size(),
190  /*vars=*/tmp_variables->data(),
191  /*weights=*/tmp_weights->data()));
192  break;
193  }
194  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
195  return absl::OkStatus();
196 }
197 
198 absl::Status AddQuadraticConstraint(
199  const MPGeneralConstraintProto& gen_cst,
200  const std::vector<SCIP_VAR*>& scip_variables, SCIP* scip,
201  SCIP_CONS** scip_cst, std::vector<SCIP_VAR*>* tmp_variables,
202  std::vector<double>* tmp_coefficients,
203  std::vector<SCIP_VAR*>* tmp_qvariables1,
204  std::vector<SCIP_VAR*>* tmp_qvariables2,
205  std::vector<double>* tmp_qcoefficients) {
206  CHECK(scip != nullptr);
207  CHECK(scip_cst != nullptr);
208  CHECK(tmp_variables != nullptr);
209  CHECK(tmp_coefficients != nullptr);
210  CHECK(tmp_qvariables1 != nullptr);
211  CHECK(tmp_qvariables2 != nullptr);
212  CHECK(tmp_qcoefficients != nullptr);
213 
214  CHECK(gen_cst.has_quadratic_constraint());
215  const MPQuadraticConstraint& quad_cst = gen_cst.quadratic_constraint();
216 
217  // Process linear part of the constraint.
218  const int lsize = quad_cst.var_index_size();
219  CHECK_EQ(quad_cst.coefficient_size(), lsize);
220  tmp_variables->resize(lsize, nullptr);
221  tmp_coefficients->resize(lsize, 0.0);
222  for (int i = 0; i < lsize; ++i) {
223  (*tmp_variables)[i] = scip_variables[quad_cst.var_index(i)];
224  (*tmp_coefficients)[i] = quad_cst.coefficient(i);
225  }
226 
227  // Process quadratic part of the constraint.
228  const int qsize = quad_cst.qvar1_index_size();
229  CHECK_EQ(quad_cst.qvar2_index_size(), qsize);
230  CHECK_EQ(quad_cst.qcoefficient_size(), qsize);
231  tmp_qvariables1->resize(qsize, nullptr);
232  tmp_qvariables2->resize(qsize, nullptr);
233  tmp_qcoefficients->resize(qsize, 0.0);
234  for (int i = 0; i < qsize; ++i) {
235  (*tmp_qvariables1)[i] = scip_variables[quad_cst.qvar1_index(i)];
236  (*tmp_qvariables2)[i] = scip_variables[quad_cst.qvar2_index(i)];
237  (*tmp_qcoefficients)[i] = quad_cst.qcoefficient(i);
238  }
239 
241  SCIPcreateConsBasicQuadratic(scip,
242  /*cons=*/scip_cst,
243  /*name=*/gen_cst.name().c_str(),
244  /*nlinvars=*/lsize,
245  /*linvars=*/tmp_variables->data(),
246  /*lincoefs=*/tmp_coefficients->data(),
247  /*nquadterms=*/qsize,
248  /*quadvars1=*/tmp_qvariables1->data(),
249  /*quadvars2=*/tmp_qvariables2->data(),
250  /*quadcoefs=*/tmp_qcoefficients->data(),
251  /*lhs=*/quad_cst.lower_bound(),
252  /*rhs=*/quad_cst.upper_bound()));
253  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
254  return absl::OkStatus();
255 }
256 
257 // Models the constraint y = |x| as y >= 0 plus one disjunction constraint:
258 // y = x OR y = -x
259 absl::Status AddAbsConstraint(const MPGeneralConstraintProto& gen_cst,
260  const std::vector<SCIP_VAR*>& scip_variables,
261  SCIP* scip, SCIP_CONS** scip_cst) {
262  CHECK(scip != nullptr);
263  CHECK(scip_cst != nullptr);
264  CHECK(gen_cst.has_abs_constraint());
265  const auto& abs = gen_cst.abs_constraint();
266  SCIP_VAR* scip_var = scip_variables[abs.var_index()];
267  SCIP_VAR* scip_resultant_var = scip_variables[abs.resultant_var_index()];
268 
269  // Set the resultant variable's lower bound to zero if it's negative.
270  if (SCIPvarGetLbLocal(scip_resultant_var) < 0.0) {
271  RETURN_IF_SCIP_ERROR(SCIPchgVarLb(scip, scip_resultant_var, 0.0));
272  }
273 
274  std::vector<SCIP_VAR*> vars;
275  std::vector<double> vals;
276  std::vector<SCIP_CONS*> cons;
277  auto add_abs_constraint =
278  [&](const std::string& name_prefix) -> absl::Status {
279  SCIP_CONS* scip_cons = nullptr;
280  CHECK(vars.size() == vals.size());
281  const std::string name =
282  gen_cst.has_name() ? absl::StrCat(gen_cst.name(), name_prefix) : "";
283  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicLinear(
284  scip, /*cons=*/&scip_cons,
285  /*name=*/name.c_str(), /*nvars=*/vars.size(), /*vars=*/vars.data(),
286  /*vals=*/vals.data(), /*lhs=*/0.0, /*rhs=*/0.0));
287  // Note that the constraints are, by design, not added into the model using
288  // SCIPaddCons.
289  cons.push_back(scip_cons);
290  return absl::OkStatus();
291  };
292 
293  // Create an intermediary constraint such that y = -x
294  vars = {scip_resultant_var, scip_var};
295  vals = {1, 1};
296  RETURN_IF_ERROR(add_abs_constraint("_neg"));
297 
298  // Create an intermediary constraint such that y = x
299  vals = {1, -1};
300  RETURN_IF_ERROR(add_abs_constraint("_pos"));
301 
302  // Activate at least one of the two above constraints.
303  const std::string name =
304  gen_cst.has_name() ? absl::StrCat(gen_cst.name(), "_disj") : "";
305  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicDisjunction(
306  scip, /*cons=*/scip_cst, /*name=*/name.c_str(),
307  /*nconss=*/cons.size(), /*conss=*/cons.data(), /*relaxcons=*/nullptr));
308  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
309 
310  return absl::OkStatus();
311 }
312 
313 absl::Status AddAndConstraint(const MPGeneralConstraintProto& gen_cst,
314  const std::vector<SCIP_VAR*>& scip_variables,
315  SCIP* scip, SCIP_CONS** scip_cst,
316  std::vector<SCIP_VAR*>* tmp_variables) {
317  CHECK(scip != nullptr);
318  CHECK(scip_cst != nullptr);
319  CHECK(tmp_variables != nullptr);
320  CHECK(gen_cst.has_and_constraint());
321  const auto& andcst = gen_cst.and_constraint();
322 
323  tmp_variables->resize(andcst.var_index_size(), nullptr);
324  for (int i = 0; i < andcst.var_index_size(); ++i) {
325  (*tmp_variables)[i] = scip_variables[andcst.var_index(i)];
326  }
327  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicAnd(
328  scip, /*cons=*/scip_cst,
329  /*name=*/gen_cst.name().c_str(),
330  /*resvar=*/scip_variables[andcst.resultant_var_index()],
331  /*nvars=*/andcst.var_index_size(),
332  /*vars=*/tmp_variables->data()));
333  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
334  return absl::OkStatus();
335 }
336 
337 absl::Status AddOrConstraint(const MPGeneralConstraintProto& gen_cst,
338  const std::vector<SCIP_VAR*>& scip_variables,
339  SCIP* scip, SCIP_CONS** scip_cst,
340  std::vector<SCIP_VAR*>* tmp_variables) {
341  CHECK(scip != nullptr);
342  CHECK(scip_cst != nullptr);
343  CHECK(tmp_variables != nullptr);
344  CHECK(gen_cst.has_or_constraint());
345  const auto& orcst = gen_cst.or_constraint();
346 
347  tmp_variables->resize(orcst.var_index_size(), nullptr);
348  for (int i = 0; i < orcst.var_index_size(); ++i) {
349  (*tmp_variables)[i] = scip_variables[orcst.var_index(i)];
350  }
351  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicOr(
352  scip, /*cons=*/scip_cst,
353  /*name=*/gen_cst.name().c_str(),
354  /*resvar=*/scip_variables[orcst.resultant_var_index()],
355  /*nvars=*/orcst.var_index_size(),
356  /*vars=*/tmp_variables->data()));
357  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
358  return absl::OkStatus();
359 }
360 
361 // Models the constraint y = min(x1, x2, ... xn, c) with c being a constant with
362 // - n + 1 constraints to ensure y <= min(x1, x2, ... xn, c)
363 // - one disjunction constraint among all of the possible y = x1, y = x2, ...
364 // y = xn, y = c constraints
365 // Does the equivalent thing for max (with y >= max(...) instead).
366 absl::Status AddMinMaxConstraint(const MPGeneralConstraintProto& gen_cst,
367  const std::vector<SCIP_VAR*>& scip_variables,
368  SCIP* scip, SCIP_CONS** scip_cst,
369  std::vector<SCIP_CONS*>* scip_constraints,
370  std::vector<SCIP_VAR*>* tmp_variables) {
371  CHECK(scip != nullptr);
372  CHECK(scip_cst != nullptr);
373  CHECK(tmp_variables != nullptr);
374  CHECK(gen_cst.has_min_constraint() || gen_cst.has_max_constraint());
375  const auto& minmax = gen_cst.has_min_constraint() ? gen_cst.min_constraint()
376  : gen_cst.max_constraint();
377  const std::set<int> unique_var_indices(minmax.var_index().begin(),
378  minmax.var_index().end());
379  SCIP_VAR* scip_resultant_var = scip_variables[minmax.resultant_var_index()];
380 
381  std::vector<SCIP_VAR*> vars;
382  std::vector<double> vals;
383  std::vector<SCIP_CONS*> cons;
384  auto add_lin_constraint = [&](const std::string& name_prefix,
385  double lower_bound = 0.0,
386  double upper_bound = 0.0) -> absl::Status {
387  SCIP_CONS* scip_cons = nullptr;
388  CHECK(vars.size() == vals.size());
389  const std::string name =
390  gen_cst.has_name() ? absl::StrCat(gen_cst.name(), name_prefix) : "";
391  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicLinear(
392  scip, /*cons=*/&scip_cons,
393  /*name=*/name.c_str(), /*nvars=*/vars.size(), /*vars=*/vars.data(),
394  /*vals=*/vals.data(), /*lhs=*/lower_bound, /*rhs=*/upper_bound));
395  // Note that the constraints are, by design, not added into the model using
396  // SCIPaddCons.
397  cons.push_back(scip_cons);
398  return absl::OkStatus();
399  };
400 
401  // Create intermediary constraints such that y = xi
402  for (const int var_index : unique_var_indices) {
403  vars = {scip_resultant_var, scip_variables[var_index]};
404  vals = {1, -1};
405  RETURN_IF_ERROR(add_lin_constraint(absl::StrCat("_", var_index)));
406  }
407 
408  // Create an intermediary constraint such that y = c
409  if (minmax.has_constant()) {
410  vars = {scip_resultant_var};
411  vals = {1};
413  add_lin_constraint("_constant", minmax.constant(), minmax.constant()));
414  }
415 
416  // Activate at least one of the above constraints.
417  const std::string name =
418  gen_cst.has_name() ? absl::StrCat(gen_cst.name(), "_disj") : "";
419  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicDisjunction(
420  scip, /*cons=*/scip_cst, /*name=*/name.c_str(),
421  /*nconss=*/cons.size(), /*conss=*/cons.data(), /*relaxcons=*/nullptr));
422  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst));
423 
424  // Add all of the inequality constraints.
425  constexpr double kInfinity = std::numeric_limits<double>::infinity();
426  cons.clear();
427  for (const int var_index : unique_var_indices) {
428  vars = {scip_resultant_var, scip_variables[var_index]};
429  vals = {1, -1};
430  if (gen_cst.has_min_constraint()) {
431  RETURN_IF_ERROR(add_lin_constraint(absl::StrCat("_ineq_", var_index),
432  -kInfinity, 0.0));
433  } else {
434  RETURN_IF_ERROR(add_lin_constraint(absl::StrCat("_ineq_", var_index), 0.0,
435  kInfinity));
436  }
437  }
438  if (minmax.has_constant()) {
439  vars = {scip_resultant_var};
440  vals = {1};
441  if (gen_cst.has_min_constraint()) {
442  RETURN_IF_ERROR(add_lin_constraint(absl::StrCat("_ineq_constant"),
443  -kInfinity, minmax.constant()));
444  } else {
445  RETURN_IF_ERROR(add_lin_constraint(absl::StrCat("_ineq_constant"),
446  minmax.constant(), kInfinity));
447  }
448  }
449  for (SCIP_CONS* scip_cons : cons) {
450  scip_constraints->push_back(scip_cons);
451  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, scip_cons));
452  }
453  return absl::OkStatus();
454 }
455 
456 absl::Status AddQuadraticObjective(const MPQuadraticObjective& quadobj,
457  SCIP* scip,
458  std::vector<SCIP_VAR*>* scip_variables,
459  std::vector<SCIP_CONS*>* scip_constraints) {
460  CHECK(scip != nullptr);
461  CHECK(scip_variables != nullptr);
462  CHECK(scip_constraints != nullptr);
463 
464  constexpr double kInfinity = std::numeric_limits<double>::infinity();
465 
466  const int size = quadobj.coefficient_size();
467  if (size == 0) return absl::OkStatus();
468 
469  // SCIP supports quadratic objectives by adding a quadratic constraint. We
470  // need to create an extra variable to hold this quadratic objective.
471  scip_variables->push_back(nullptr);
472  RETURN_IF_SCIP_ERROR(SCIPcreateVarBasic(scip, /*var=*/&scip_variables->back(),
473  /*name=*/"quadobj",
474  /*lb=*/-kInfinity, /*ub=*/kInfinity,
475  /*obj=*/1,
476  /*vartype=*/SCIP_VARTYPE_CONTINUOUS));
477  RETURN_IF_SCIP_ERROR(SCIPaddVar(scip, scip_variables->back()));
478 
479  scip_constraints->push_back(nullptr);
480  SCIP_VAR* linvars[1] = {scip_variables->back()};
481  double lincoefs[1] = {-1};
482  std::vector<SCIP_VAR*> quadvars1(size, nullptr);
483  std::vector<SCIP_VAR*> quadvars2(size, nullptr);
484  std::vector<double> quadcoefs(size, 0);
485  for (int i = 0; i < size; ++i) {
486  quadvars1[i] = scip_variables->at(quadobj.qvar1_index(i));
487  quadvars2[i] = scip_variables->at(quadobj.qvar2_index(i));
488  quadcoefs[i] = quadobj.coefficient(i);
489  }
490  RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicQuadratic(
491  scip, /*cons=*/&scip_constraints->back(), /*name=*/"quadobj",
492  /*nlinvars=*/1, /*linvars=*/linvars, /*lincoefs=*/lincoefs,
493  /*nquadterms=*/size, /*quadvars1=*/quadvars1.data(),
494  /*quadvars2=*/quadvars2.data(), /*quadcoefs=*/quadcoefs.data(),
495  /*lhs=*/0, /*rhs=*/0));
496  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, scip_constraints->back()));
497 
498  return absl::OkStatus();
499 }
500 
501 absl::Status AddSolutionHint(const MPModelProto& model, SCIP* scip,
502  const std::vector<SCIP_VAR*>& scip_variables) {
503  CHECK(scip != nullptr);
504  if (!model.has_solution_hint()) return absl::OkStatus();
505 
506  const PartialVariableAssignment& solution_hint = model.solution_hint();
507  SCIP_SOL* solution;
508  bool is_solution_partial =
509  solution_hint.var_index_size() != model.variable_size();
510  if (is_solution_partial) {
512  SCIPcreatePartialSol(scip, /*sol=*/&solution, /*heur=*/nullptr));
513  } else {
515  SCIPcreateSol(scip, /*sol=*/&solution, /*heur=*/nullptr));
516  }
517 
518  for (int i = 0; i < solution_hint.var_index_size(); ++i) {
519  RETURN_IF_SCIP_ERROR(SCIPsetSolVal(
520  scip, solution, scip_variables[solution_hint.var_index(i)],
521  solution_hint.var_value(i)));
522  }
523 
524  SCIP_Bool is_stored;
525  RETURN_IF_SCIP_ERROR(SCIPaddSolFree(scip, &solution, &is_stored));
526 
527  return absl::OkStatus();
528 }
529 } // namespace
530 
531 // Returns "" iff the model seems valid for SCIP, else returns a human-readable
532 // error message. Assumes that FindErrorInMPModelProto(model) found no error.
533 std::string FindErrorInMPModelForScip(const MPModelProto& model, SCIP* scip) {
534  CHECK(scip != nullptr);
535  const double infinity = SCIPinfinity(scip);
536 
537  for (int v = 0; v < model.variable_size(); ++v) {
538  const MPVariableProto& variable = model.variable(v);
539  if (variable.lower_bound() >= infinity) {
540  return absl::StrFormat(
541  "Variable %i's lower bound is considered +infinity", v);
542  }
543  if (variable.upper_bound() <= -infinity) {
544  return absl::StrFormat(
545  "Variable %i's upper bound is considered -infinity", v);
546  }
547  const double coeff = variable.objective_coefficient();
548  if (coeff >= infinity || coeff <= -infinity) {
549  return absl::StrFormat(
550  "Variable %i's objective coefficient is considered infinite", v);
551  }
552  }
553 
554  for (int c = 0; c < model.constraint_size(); ++c) {
555  const MPConstraintProto& cst = model.constraint(c);
556  if (cst.lower_bound() >= infinity) {
557  return absl::StrFormat(
558  "Constraint %d's lower_bound is considered +infinity", c);
559  }
560  if (cst.upper_bound() <= -infinity) {
561  return absl::StrFormat(
562  "Constraint %d's upper_bound is considered -infinity", c);
563  }
564  for (int i = 0; i < cst.coefficient_size(); ++i) {
565  if (std::abs(cst.coefficient(i)) >= infinity) {
566  return absl::StrFormat(
567  "Constraint %d's coefficient #%d is considered infinite", c, i);
568  }
569  }
570  }
571 
572  for (int c = 0; c < model.general_constraint_size(); ++c) {
573  const MPGeneralConstraintProto& cst = model.general_constraint(c);
574  switch (cst.general_constraint_case()) {
575  case MPGeneralConstraintProto::kQuadraticConstraint:
576  if (cst.quadratic_constraint().lower_bound() >= infinity) {
577  return absl::StrFormat(
578  "Quadratic constraint %d's lower_bound is considered +infinity",
579  c);
580  }
581  if (cst.quadratic_constraint().upper_bound() <= -infinity) {
582  return absl::StrFormat(
583  "Quadratic constraint %d's upper_bound is considered -infinity",
584  c);
585  }
586  for (int i = 0; i < cst.quadratic_constraint().coefficient_size();
587  ++i) {
588  const double coefficient = cst.quadratic_constraint().coefficient(i);
589  if (coefficient >= infinity || coefficient <= -infinity) {
590  return absl::StrFormat(
591  "Quadratic constraint %d's linear coefficient #%d considered "
592  "infinite",
593  c, i);
594  }
595  }
596  for (int i = 0; i < cst.quadratic_constraint().qcoefficient_size();
597  ++i) {
598  const double qcoefficient =
599  cst.quadratic_constraint().qcoefficient(i);
600  if (qcoefficient >= infinity || qcoefficient <= -infinity) {
601  return absl::StrFormat(
602  "Quadratic constraint %d's quadratic coefficient #%d "
603  "considered infinite",
604  c, i);
605  }
606  }
607  break;
608  case MPGeneralConstraintProto::kMinConstraint:
609  if (cst.min_constraint().constant() >= infinity ||
610  cst.min_constraint().constant() <= -infinity) {
611  return absl::StrFormat(
612  "Min constraint %d's coefficient constant considered infinite",
613  c);
614  }
615  break;
616  case MPGeneralConstraintProto::kMaxConstraint:
617  if (cst.max_constraint().constant() >= infinity ||
618  cst.max_constraint().constant() <= -infinity) {
619  return absl::StrFormat(
620  "Max constraint %d's coefficient constant considered infinite",
621  c);
622  }
623  break;
624  default:
625  continue;
626  }
627  }
628 
629  const MPQuadraticObjective& quad_obj = model.quadratic_objective();
630  for (int i = 0; i < quad_obj.coefficient_size(); ++i) {
631  if (std::abs(quad_obj.coefficient(i)) >= infinity) {
632  return absl::StrFormat(
633  "Quadratic objective term #%d's coefficient is considered infinite",
634  i);
635  }
636  }
637 
638  if (model.has_solution_hint()) {
639  for (int i = 0; i < model.solution_hint().var_value_size(); ++i) {
640  const double value = model.solution_hint().var_value(i);
641  if (value >= infinity || value <= -infinity) {
642  return absl::StrFormat(
643  "Variable %i's solution hint is considered infinite",
644  model.solution_hint().var_index(i));
645  }
646  }
647  }
648 
649  if (model.objective_offset() >= infinity ||
650  model.objective_offset() <= -infinity) {
651  return "Model's objective offset is considered infinite.";
652  }
653 
654  return "";
655 }
656 
657 absl::StatusOr<MPSolutionResponse> ScipSolveProto(
658  const MPModelRequest& request) {
659  MPSolutionResponse response;
660  const absl::optional<LazyMutableCopy<MPModelProto>> optional_model =
662  if (!optional_model) return response;
663  const MPModelProto& model = optional_model->get();
664  SCIP* scip = nullptr;
665  std::vector<SCIP_VAR*> scip_variables(model.variable_size(), nullptr);
666  std::vector<SCIP_CONS*> scip_constraints(
667  model.constraint_size() + model.general_constraint_size(), nullptr);
668 
669  auto delete_scip_objects = [&]() -> absl::Status {
670  // Release all created pointers.
671  if (scip == nullptr) return absl::OkStatus();
672  for (SCIP_VAR* variable : scip_variables) {
673  if (variable != nullptr) {
674  RETURN_IF_SCIP_ERROR(SCIPreleaseVar(scip, &variable));
675  }
676  }
677  for (SCIP_CONS* constraint : scip_constraints) {
678  if (constraint != nullptr) {
679  RETURN_IF_SCIP_ERROR(SCIPreleaseCons(scip, &constraint));
680  }
681  }
682  RETURN_IF_SCIP_ERROR(SCIPfree(&scip));
683  return absl::OkStatus();
684  };
685 
686  auto scip_deleter = absl::MakeCleanup([delete_scip_objects]() {
687  const absl::Status deleter_status = delete_scip_objects();
688  LOG_IF(DFATAL, !deleter_status.ok()) << deleter_status;
689  });
690 
691  RETURN_IF_SCIP_ERROR(SCIPcreate(&scip));
692  RETURN_IF_SCIP_ERROR(SCIPincludeDefaultPlugins(scip));
693  const std::string scip_model_invalid_error =
695  if (!scip_model_invalid_error.empty()) {
696  response.set_status(MPSOLVER_MODEL_INVALID);
697  response.set_status_str(scip_model_invalid_error);
698  return response;
699  }
700 
701  const auto parameters_status = LegacyScipSetSolverSpecificParameters(
702  request.solver_specific_parameters(), scip);
703  if (!parameters_status.ok()) {
704  response.set_status(MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS);
705  response.set_status_str(
706  std::string(parameters_status.message())); // NOLINT
707  return response;
708  }
709  // Default clock type. We use wall clock time because getting CPU user seconds
710  // involves calling times() which is very expensive.
711  // NOTE(user): Also, time limit based on CPU user seconds is *NOT* thread
712  // safe. We observed that different instances of SCIP running concurrently
713  // in different threads consume the time limit *together*. E.g., 2 threads
714  // running SCIP with time limit 10s each will both terminate after ~5s.
716  SCIPsetIntParam(scip, "timing/clocktype", SCIP_CLOCKTYPE_WALL));
717  if (request.solver_time_limit_seconds() > 0 &&
718  request.solver_time_limit_seconds() < 1e20) {
719  RETURN_IF_SCIP_ERROR(SCIPsetRealParam(scip, "limits/time",
720  request.solver_time_limit_seconds()));
721  }
722  SCIPsetMessagehdlrQuiet(scip, !request.enable_internal_solver_output());
723 
724  RETURN_IF_SCIP_ERROR(SCIPcreateProbBasic(scip, model.name().c_str()));
725  if (model.maximize()) {
726  RETURN_IF_SCIP_ERROR(SCIPsetObjsense(scip, SCIP_OBJSENSE_MAXIMIZE));
727  }
728 
729  for (int v = 0; v < model.variable_size(); ++v) {
730  const MPVariableProto& variable = model.variable(v);
731  RETURN_IF_SCIP_ERROR(SCIPcreateVarBasic(
732  scip, /*var=*/&scip_variables[v], /*name=*/variable.name().c_str(),
733  /*lb=*/variable.lower_bound(), /*ub=*/variable.upper_bound(),
734  /*obj=*/variable.objective_coefficient(),
735  /*vartype=*/variable.is_integer() ? SCIP_VARTYPE_INTEGER
736  : SCIP_VARTYPE_CONTINUOUS));
737  RETURN_IF_SCIP_ERROR(SCIPaddVar(scip, scip_variables[v]));
738  }
739 
740  {
741  std::vector<SCIP_VAR*> ct_variables;
742  std::vector<double> ct_coefficients;
743  for (int c = 0; c < model.constraint_size(); ++c) {
744  const MPConstraintProto& constraint = model.constraint(c);
745  const int size = constraint.var_index_size();
746  ct_variables.resize(size, nullptr);
747  ct_coefficients.resize(size, 0);
748  for (int i = 0; i < size; ++i) {
749  ct_variables[i] = scip_variables[constraint.var_index(i)];
750  ct_coefficients[i] = constraint.coefficient(i);
751  }
752  RETURN_IF_SCIP_ERROR(SCIPcreateConsLinear(
753  scip, /*cons=*/&scip_constraints[c],
754  /*name=*/constraint.name().c_str(),
755  /*nvars=*/constraint.var_index_size(), /*vars=*/ct_variables.data(),
756  /*vals=*/ct_coefficients.data(),
757  /*lhs=*/constraint.lower_bound(), /*rhs=*/constraint.upper_bound(),
758  /*initial=*/!constraint.is_lazy(),
759  /*separate=*/true,
760  /*enforce=*/true,
761  /*check=*/true,
762  /*propagate=*/true,
763  /*local=*/false,
764  /*modifiable=*/false,
765  /*dynamic=*/false,
766  /*removable=*/constraint.is_lazy(),
767  /*stickingatnode=*/false));
768  RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, scip_constraints[c]));
769  }
770 
771  // These extra arrays are used by quadratic constraints.
772  std::vector<SCIP_VAR*> ct_qvariables1;
773  std::vector<SCIP_VAR*> ct_qvariables2;
774  std::vector<double> ct_qcoefficients;
775  const int lincst_size = model.constraint_size();
776  for (int c = 0; c < model.general_constraint_size(); ++c) {
777  const MPGeneralConstraintProto& gen_cst = model.general_constraint(c);
778  switch (gen_cst.general_constraint_case()) {
779  case MPGeneralConstraintProto::kIndicatorConstraint: {
780  RETURN_IF_ERROR(AddIndicatorConstraint(
781  gen_cst, scip, &scip_constraints[lincst_size + c],
782  &scip_variables, &scip_constraints, &ct_variables,
783  &ct_coefficients));
784  break;
785  }
786  case MPGeneralConstraintProto::kSosConstraint: {
787  RETURN_IF_ERROR(AddSosConstraint(gen_cst, scip_variables, scip,
788  &scip_constraints[lincst_size + c],
789  &ct_variables, &ct_coefficients));
790  break;
791  }
792  case MPGeneralConstraintProto::kQuadraticConstraint: {
793  RETURN_IF_ERROR(AddQuadraticConstraint(
794  gen_cst, scip_variables, scip, &scip_constraints[lincst_size + c],
795  &ct_variables, &ct_coefficients, &ct_qvariables1, &ct_qvariables2,
796  &ct_qcoefficients));
797  break;
798  }
799  case MPGeneralConstraintProto::kAbsConstraint: {
800  RETURN_IF_ERROR(AddAbsConstraint(gen_cst, scip_variables, scip,
801  &scip_constraints[lincst_size + c]));
802  break;
803  }
804  case MPGeneralConstraintProto::kAndConstraint: {
805  RETURN_IF_ERROR(AddAndConstraint(gen_cst, scip_variables, scip,
806  &scip_constraints[lincst_size + c],
807  &ct_variables));
808  break;
809  }
810  case MPGeneralConstraintProto::kOrConstraint: {
811  RETURN_IF_ERROR(AddOrConstraint(gen_cst, scip_variables, scip,
812  &scip_constraints[lincst_size + c],
813  &ct_variables));
814  break;
815  }
816  case MPGeneralConstraintProto::kMinConstraint:
817  case MPGeneralConstraintProto::kMaxConstraint: {
818  RETURN_IF_ERROR(AddMinMaxConstraint(
819  gen_cst, scip_variables, scip, &scip_constraints[lincst_size + c],
820  &scip_constraints, &ct_variables));
821  break;
822  }
823  default:
824  return absl::UnimplementedError(
825  absl::StrFormat("General constraints of type %i not supported.",
826  gen_cst.general_constraint_case()));
827  }
828  }
829  }
830 
831  if (model.has_quadratic_objective()) {
832  RETURN_IF_ERROR(AddQuadraticObjective(model.quadratic_objective(), scip,
833  &scip_variables, &scip_constraints));
834  }
835  RETURN_IF_SCIP_ERROR(SCIPaddOrigObjoffset(scip, model.objective_offset()));
836  RETURN_IF_ERROR(AddSolutionHint(model, scip, scip_variables));
837 
838  if (!absl::GetFlag(FLAGS_scip_proto_solver_output_cip_file).empty()) {
839  SCIPwriteOrigProblem(
840  scip, absl::GetFlag(FLAGS_scip_proto_solver_output_cip_file).c_str(),
841  nullptr, true);
842  }
843  RETURN_IF_SCIP_ERROR(SCIPsolve(scip));
844 
845  SCIP_SOL* const solution = SCIPgetBestSol(scip);
846  if (solution != nullptr) {
847  response.set_objective_value(SCIPgetSolOrigObj(scip, solution));
848  response.set_best_objective_bound(SCIPgetDualbound(scip));
849  for (int v = 0; v < model.variable_size(); ++v) {
850  double value = SCIPgetSolVal(scip, solution, scip_variables[v]);
851  if (model.variable(v).is_integer()) value = std::round(value);
852  response.add_variable_value(value);
853  }
854  }
855 
856  const SCIP_STATUS scip_status = SCIPgetStatus(scip);
857  switch (scip_status) {
858  case SCIP_STATUS_OPTIMAL:
859  response.set_status(MPSOLVER_OPTIMAL);
860  break;
861  case SCIP_STATUS_GAPLIMIT:
862  // To be consistent with the other solvers.
863  response.set_status(MPSOLVER_OPTIMAL);
864  break;
865  case SCIP_STATUS_INFORUNBD:
866  // NOTE(user): After looking at the SCIP code on 2019-06-14, it seems
867  // that this will mostly happen for INFEASIBLE problems in practice.
868  // Since most (all?) users shouldn't have their application behave very
869  // differently upon INFEASIBLE or UNBOUNDED, the potential error that we
870  // are making here seems reasonable (and not worth a LOG, unless in
871  // debug mode).
872  DLOG(INFO) << "SCIP solve returned SCIP_STATUS_INFORUNBD, which we treat "
873  "as INFEASIBLE even though it may mean UNBOUNDED.";
874  response.set_status_str(
875  "The model may actually be unbounded: SCIP returned "
876  "SCIP_STATUS_INFORUNBD");
877  ABSL_FALLTHROUGH_INTENDED;
878  case SCIP_STATUS_INFEASIBLE:
879  response.set_status(MPSOLVER_INFEASIBLE);
880  break;
881  case SCIP_STATUS_UNBOUNDED:
882  response.set_status(MPSOLVER_UNBOUNDED);
883  break;
884  default:
885  if (solution != nullptr) {
886  response.set_status(MPSOLVER_FEASIBLE);
887  } else {
888  response.set_status(MPSOLVER_NOT_SOLVED);
889  response.set_status_str(absl::StrFormat("SCIP status code %d",
890  static_cast<int>(scip_status)));
891  }
892  break;
893  }
894 
895  VLOG(1) << "ScipSolveProto() status="
896  << MPSolverResponseStatus_Name(response.status()) << ".";
897  return response;
898 }
899 
900 } // namespace operations_research
901 
902 #endif // #if defined(USE_SCIP)
#define LOG_IF(severity, condition)
Definition: base/logging.h:479
#define DLOG(severity)
Definition: base/logging.h:875
#define CHECK(condition)
Definition: base/logging.h:495
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define VLOG(verboselevel)
Definition: base/logging.h:978
SharedResponseManager * response
const std::string name
int64 value
GRBmodel * model
const int INFO
Definition: log_severity.h:31
absl::Cleanup< absl::decay_t< Callback > > MakeCleanup(Callback &&callback)
Definition: cleanup.h:120
const double kInfinity
Definition: lp_types.h:83
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::string FindErrorInMPModelForScip(const MPModelProto &model, SCIP *scip)
absl::StatusOr< MPSolutionResponse > ScipSolveProto(const MPModelRequest &request)
absl::Status LegacyScipSetSolverSpecificParameters(const std::string &parameters, SCIP *scip)
absl::optional< LazyMutableCopy< MPModelProto > > ExtractValidMPModelOrPopulateResponseStatus(const MPModelRequest &request, MPSolutionResponse *response)
If the model is valid and non-empty, returns it (possibly after extracting the model_delta).
int64 coefficient
#define RETURN_IF_SCIP_ERROR(x)
ABSL_FLAG(std::string, scip_proto_solver_output_cip_file, "", "If given, saves the generated CIP file here. Useful for " "reporting bugs to SCIP.")
#define RETURN_IF_ERROR(expr)
Definition: status_macros.h:27