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"
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
55ABSL_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
59namespace operations_research {
60
61namespace {
62// This function will create a new constraint if the indicator constraint has
63// both a lower bound and an upper bound.
64absl::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
138absl::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
198absl::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
259absl::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
313absl::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
337absl::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).
366absl::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
456absl::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
501absl::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.
533std::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
657absl::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