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