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