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