Google OR-Tools v9.15
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
gurobi_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
15
16#include <algorithm>
17#include <cmath>
18#include <cstddef>
19#include <cstdint>
20#include <limits>
21#include <memory>
22#include <optional>
23#include <string>
24#include <string_view>
25#include <utility>
26#include <vector>
27
28#include "absl/algorithm/container.h"
29#include "absl/base/nullability.h"
30#include "absl/container/flat_hash_map.h"
31#include "absl/container/flat_hash_set.h"
32#include "absl/log/check.h"
33#include "absl/log/log.h"
34#include "absl/memory/memory.h"
35#include "absl/status/status.h"
36#include "absl/status/statusor.h"
37#include "absl/strings/escaping.h"
38#include "absl/strings/str_cat.h"
39#include "absl/strings/str_join.h"
40#include "absl/strings/string_view.h"
41#include "absl/time/clock.h"
42#include "absl/time/time.h"
43#include "absl/types/span.h"
44#include "google/protobuf/repeated_ptr_field.h"
73
74namespace operations_research {
75namespace math_opt {
76namespace {
77
78constexpr SupportedProblemStructures kGurobiSupportedStructures = {
79 .integer_variables = SupportType::kSupported,
80 .multi_objectives = SupportType::kSupported,
81 .quadratic_objectives = SupportType::kSupported,
82 .quadratic_constraints = SupportType::kSupported,
83 .second_order_cone_constraints = SupportType::kSupported,
84 .sos1_constraints = SupportType::kSupported,
85 .sos2_constraints = SupportType::kSupported,
86 .indicator_constraints = SupportType::kSupported};
87
88absl::StatusOr<std::unique_ptr<Gurobi>> GurobiFromInitArgs(
89 const SolverInterface::InitArgs& init_args) {
90 if (kAnyXsanEnabled) {
91 return absl::FailedPreconditionError(
92 "the Gurobi library is not compatible with any sanitizer (MSAN, ASAN "
93 "or TSAN)");
94 }
95
96 // We don't test or return an error for incorrect non streamable arguments
97 // type since it is already tested by the Solver class.
98 const NonStreamableGurobiInitArguments* const non_streamable_args =
99 init_args.non_streamable != nullptr
100 ? init_args.non_streamable->ToNonStreamableGurobiInitArguments()
101 : nullptr;
102 std::unique_ptr<Gurobi> gurobi;
103 if (non_streamable_args != nullptr &&
104 non_streamable_args->primary_env != nullptr) {
105 return Gurobi::NewWithSharedPrimaryEnv(non_streamable_args->primary_env);
106 } else if (init_args.streamable.has_gurobi() &&
107 init_args.streamable.gurobi().has_isv_key()) {
109 GRBenvUniquePtr env,
110 NewPrimaryEnvironment(init_args.streamable.gurobi().isv_key()));
111 return Gurobi::New(std::move(env));
112 } else {
113 return Gurobi::New();
114 }
115}
116
117inline BasisStatusProto ConvertVariableStatus(const int status) {
118 switch (status) {
119 case GRB_BASIC:
120 return BASIS_STATUS_BASIC;
125 case GRB_SUPERBASIC:
126 return BASIS_STATUS_FREE;
127 default:
129 }
130}
131
132inline int GrbVariableStatus(const BasisStatusProto status) {
133 switch (status) {
135 return GRB_BASIC;
138 return GRB_NONBASIC_LOWER;
140 return GRB_NONBASIC_UPPER;
142 return GRB_SUPERBASIC;
144 default:
145 LOG(FATAL) << "Unexpected invalid initial_basis.";
146 return 0;
147 }
148}
149
150// is_mip indicates if the problem has integer variables or constraints that
151// would cause Gurobi to treat the problem as a MIP, e.g. SOS, indicator.
152absl::StatusOr<GurobiParametersProto> MergeParameters(
153 const SolveParametersProto& solve_parameters,
154 const ModelSolveParametersProto& model_parameters, const bool is_mip,
155 const bool is_multi_objective_mode) {
156 GurobiParametersProto merged_parameters;
157
158 {
159 GurobiParametersProto::Parameter* const parameter =
160 merged_parameters.add_parameters();
161 parameter->set_name(GRB_INT_PAR_LOGTOCONSOLE);
162 parameter->set_value(solve_parameters.enable_output() ? "1" : "0");
163 }
164
165 {
166 absl::Duration time_limit = absl::InfiniteDuration();
167 if (solve_parameters.has_time_limit()) {
169 solve_parameters.time_limit()));
170 }
171 // If we do not have a multi-objective model, but the user has specified a
172 // time limit on the primary objective, this is functionally a second way to
173 // specify the global time limit. So, we should take the min of the two.
174 if (!is_multi_objective_mode &&
175 model_parameters.primary_objective_parameters().has_time_limit()) {
177 const absl::Duration primary_objective_time_limit,
179 model_parameters.primary_objective_parameters().time_limit()));
180 time_limit = std::min(time_limit, primary_objective_time_limit);
181 }
182 if (time_limit < absl::InfiniteDuration()) {
183 GurobiParametersProto::Parameter* const parameter =
184 merged_parameters.add_parameters();
185 parameter->set_name(GRB_DBL_PAR_TIMELIMIT);
186 parameter->set_value(absl::StrCat(absl::ToDoubleSeconds(time_limit)));
187 }
188 }
189
190 if (solve_parameters.has_node_limit()) {
191 GurobiParametersProto::Parameter* const parameter =
192 merged_parameters.add_parameters();
193 parameter->set_name(GRB_DBL_PAR_NODELIMIT);
194 parameter->set_value(absl::StrCat(solve_parameters.node_limit()));
195 }
196
197 if (solve_parameters.has_threads()) {
198 const int threads = solve_parameters.threads();
199 GurobiParametersProto::Parameter* const parameter =
200 merged_parameters.add_parameters();
201 parameter->set_name(GRB_INT_PAR_THREADS);
202 parameter->set_value(absl::StrCat(threads));
203 }
204
205 if (solve_parameters.has_absolute_gap_tolerance()) {
206 const double absolute_gap_tolerance =
207 solve_parameters.absolute_gap_tolerance();
208 GurobiParametersProto::Parameter* const parameter =
209 merged_parameters.add_parameters();
210 parameter->set_name(GRB_DBL_PAR_MIPGAPABS);
211 parameter->set_value(absl::StrCat(absolute_gap_tolerance));
212 }
213
214 if (solve_parameters.has_relative_gap_tolerance()) {
215 const double relative_gap_tolerance =
216 solve_parameters.relative_gap_tolerance();
217 GurobiParametersProto::Parameter* const parameter =
218 merged_parameters.add_parameters();
219 parameter->set_name(GRB_DBL_PAR_MIPGAP);
220 parameter->set_value(absl::StrCat(relative_gap_tolerance));
221 }
222
223 if (solve_parameters.has_cutoff_limit()) {
224 GurobiParametersProto::Parameter* const parameter =
225 merged_parameters.add_parameters();
226 parameter->set_name(GRB_DBL_PAR_CUTOFF);
227 parameter->set_value(absl::StrCat(solve_parameters.cutoff_limit()));
228 }
229
230 if (solve_parameters.has_objective_limit()) {
231 if (!is_mip) {
232 return absl::InvalidArgumentError(
233 "objective_limit is only supported for Gurobi on MIP models");
234 }
235 GurobiParametersProto::Parameter* const parameter =
236 merged_parameters.add_parameters();
237 parameter->set_name(GRB_DBL_PAR_BESTOBJSTOP);
238 parameter->set_value(absl::StrCat(solve_parameters.objective_limit()));
239 }
240
241 if (solve_parameters.has_best_bound_limit()) {
242 if (!is_mip) {
243 return absl::InvalidArgumentError(
244 "best_bound_limit is only supported for Gurobi on MIP models");
245 }
246 GurobiParametersProto::Parameter* const parameter =
247 merged_parameters.add_parameters();
248 parameter->set_name(GRB_DBL_PAR_BESTBDSTOP);
249 parameter->set_value(absl::StrCat(solve_parameters.best_bound_limit()));
250 }
251
252 if (solve_parameters.has_solution_limit()) {
253 GurobiParametersProto::Parameter* const parameter =
254 merged_parameters.add_parameters();
255 parameter->set_name(GRB_INT_PAR_SOLUTIONLIMIT);
256 parameter->set_value(absl::StrCat(solve_parameters.solution_limit()));
257 }
258
259 if (solve_parameters.has_random_seed()) {
260 const int random_seed =
261 std::min(GRB_MAXINT, std::max(solve_parameters.random_seed(), 0));
262 GurobiParametersProto::Parameter* const parameter =
263 merged_parameters.add_parameters();
264 parameter->set_name(GRB_INT_PAR_SEED);
265 parameter->set_value(absl::StrCat(random_seed));
266 }
267
268 if (solve_parameters.has_solution_pool_size()) {
269 GurobiParametersProto::Parameter* const solution_pool_size =
270 merged_parameters.add_parameters();
271 solution_pool_size->set_name(GRB_INT_PAR_POOLSOLUTIONS);
272 solution_pool_size->set_value(
273 absl::StrCat(solve_parameters.solution_pool_size()));
274 }
275
276 if (solve_parameters.lp_algorithm() != LP_ALGORITHM_UNSPECIFIED) {
277 GurobiParametersProto::Parameter* const parameter =
278 merged_parameters.add_parameters();
279 parameter->set_name(GRB_INT_PAR_METHOD);
280 switch (solve_parameters.lp_algorithm()) {
282 parameter->set_value(absl::StrCat(GRB_METHOD_PRIMAL));
283 break;
285 parameter->set_value(absl::StrCat(GRB_METHOD_DUAL));
286 break;
288 parameter->set_value(absl::StrCat(GRB_METHOD_BARRIER));
289 break;
291 return absl::InvalidArgumentError(
292 "lp_algorithm FIRST_ORDER is not supported for gurobi");
293 default:
294 LOG(FATAL) << "LPAlgorithm: "
295 << ProtoEnumToString(solve_parameters.lp_algorithm())
296 << " unknown, error setting Gurobi parameters";
297 }
298 }
299
300 if (solve_parameters.scaling() != EMPHASIS_UNSPECIFIED) {
301 GurobiParametersProto::Parameter* const parameter =
302 merged_parameters.add_parameters();
303 parameter->set_name(GRB_INT_PAR_SCALEFLAG);
304 switch (solve_parameters.scaling()) {
305 case EMPHASIS_OFF:
306 parameter->set_value(absl::StrCat(0));
307 break;
308 case EMPHASIS_LOW:
309 case EMPHASIS_MEDIUM:
310 parameter->set_value(absl::StrCat(1));
311 break;
312 case EMPHASIS_HIGH:
313 parameter->set_value(absl::StrCat(2));
314 break;
316 parameter->set_value(absl::StrCat(3));
317 break;
318 default:
319 LOG(FATAL) << "Scaling emphasis: "
320 << ProtoEnumToString(solve_parameters.scaling())
321 << " unknown, error setting Gurobi parameters";
322 }
323 }
324
325 if (solve_parameters.cuts() != EMPHASIS_UNSPECIFIED) {
326 GurobiParametersProto::Parameter* const parameter =
327 merged_parameters.add_parameters();
328 parameter->set_name(GRB_INT_PAR_CUTS);
329 switch (solve_parameters.cuts()) {
330 case EMPHASIS_OFF:
331 parameter->set_value(absl::StrCat(0));
332 break;
333 case EMPHASIS_LOW:
334 case EMPHASIS_MEDIUM:
335 parameter->set_value(absl::StrCat(1));
336 break;
337 case EMPHASIS_HIGH:
338 parameter->set_value(absl::StrCat(2));
339 break;
341 parameter->set_value(absl::StrCat(3));
342 break;
343 default:
344 LOG(FATAL) << "Cuts emphasis: "
345 << ProtoEnumToString(solve_parameters.cuts())
346 << " unknown, error setting Gurobi parameters";
347 }
348 }
349
350 if (solve_parameters.heuristics() != EMPHASIS_UNSPECIFIED) {
351 GurobiParametersProto::Parameter* const parameter =
352 merged_parameters.add_parameters();
353 parameter->set_name(GRB_DBL_PAR_HEURISTICS);
354 switch (solve_parameters.heuristics()) {
355 case EMPHASIS_OFF:
356 parameter->set_value(absl::StrCat(0.0));
357 break;
358 case EMPHASIS_LOW:
359 parameter->set_value(absl::StrCat(0.025));
360 break;
361 case EMPHASIS_MEDIUM:
362 // As of Gurobi 9.1 this is the default value.
363 // https://www.gurobi.com/documentation/9.1/refman/heuristics.html
364 parameter->set_value(absl::StrCat(0.05));
365 break;
366 case EMPHASIS_HIGH:
367 parameter->set_value(absl::StrCat(0.1));
368 break;
370 parameter->set_value(absl::StrCat(0.2));
371 break;
372 default:
373 LOG(FATAL) << "Heuristics emphasis: "
374 << ProtoEnumToString(solve_parameters.heuristics())
375 << " unknown, error setting Gurobi parameters";
376 }
377 }
378
379 if (solve_parameters.presolve() != EMPHASIS_UNSPECIFIED) {
380 GurobiParametersProto::Parameter* const parameter =
381 merged_parameters.add_parameters();
382 parameter->set_name(GRB_INT_PAR_PRESOLVE);
383 switch (solve_parameters.presolve()) {
384 case EMPHASIS_OFF:
385 parameter->set_value(absl::StrCat(0));
386 break;
387 case EMPHASIS_LOW:
388 case EMPHASIS_MEDIUM:
389 parameter->set_value(absl::StrCat(1));
390 break;
391 case EMPHASIS_HIGH:
393 parameter->set_value(absl::StrCat(2));
394 break;
395 default:
396 LOG(FATAL) << "Presolve emphasis: "
397 << ProtoEnumToString(solve_parameters.presolve())
398 << " unknown, error setting Gurobi parameters";
399 }
400 }
401
402 if (solve_parameters.has_iteration_limit()) {
403 GurobiParametersProto::Parameter* const iterationlimit =
404 merged_parameters.add_parameters();
405 iterationlimit->set_name(GRB_DBL_PAR_ITERATIONLIMIT);
406 iterationlimit->set_value(absl::StrCat(solve_parameters.iteration_limit()));
407 GurobiParametersProto::Parameter* const bariterlimit =
408 merged_parameters.add_parameters();
409 bariterlimit->set_name(GRB_INT_PAR_BARITERLIMIT);
410 double val = std::min<double>(std::numeric_limits<int32_t>::max(),
411 solve_parameters.iteration_limit());
412 bariterlimit->set_value(absl::StrCat(val));
413 }
414
415 for (const GurobiParametersProto::Parameter& parameter :
416 solve_parameters.gurobi().parameters()) {
417 *merged_parameters.add_parameters() = parameter;
418 }
419
420 return merged_parameters;
421}
422
423absl::StatusOr<int64_t> SafeInt64FromDouble(const double d) {
424 const int64_t result = static_cast<int64_t>(d);
425 if (static_cast<double>(result) != d) {
426 return absl::InternalError(
427 absl::StrCat("Expected double ", d, " to contain an int64_t."));
428 }
429 return result;
430}
431
432const absl::flat_hash_set<CallbackEventProto>& SupportedMIPEvents() {
433 static const auto* const kEvents =
434 new absl::flat_hash_set<CallbackEventProto>({
440 // CALLBACK_EVENT_BARRIER is not supported when solving MIPs; it turns
441 // out that Gurobi uses a barrier algorithm to solve the root node
442 // relaxation (from the traces) but does not call the associated
443 // callback.
444 });
445 return *kEvents;
446}
447
448const absl::flat_hash_set<CallbackEventProto>& SupportedLPEvents() {
449 static const auto* const kEvents =
450 new absl::flat_hash_set<CallbackEventProto>({
454 });
455 return *kEvents;
456}
457
458// Gurobi names (model, variables and constraints) must be no longer than 255
459// characters; or Gurobi fails with an error.
460constexpr std::size_t kMaxNameSize = 255;
461
462// Returns a string of at most kMaxNameSize max size.
463std::string TruncateName(const std::string_view original_name) {
464 return std::string(
465 original_name.substr(0, std::min(kMaxNameSize, original_name.size())));
466}
467
468// Truncate the names of variables and constraints.
469std::vector<std::string> TruncateNames(
470 const google::protobuf::RepeatedPtrField<std::string>& original_names) {
471 std::vector<std::string> result;
472 result.reserve(original_names.size());
473 for (const std::string& original_name : original_names) {
474 result.push_back(TruncateName(original_name));
475 }
476 return result;
477}
478
479absl::Status SafeGurobiDouble(const double d) {
480 if (std::isfinite(d) && std::abs(d) >= GRB_INFINITY) {
482 << "finite value: " << d << " will be treated as infinite by Gurobi";
483 }
484 return absl::OkStatus();
485}
486
487std::string EscapedNameForLogging(const absl::string_view name) {
488 return absl::StrCat("\"", absl::Utf8SafeCEscape(name), "\"");
489}
490
491constexpr int kDeletedIndex = -1;
492constexpr int kUnsetIndex = -2;
493// Returns a vector of length `size_before_delete` that logically provides a
494// mapping from the starting contiguous range [0, ..., size_before_delete) to
495// a potentially smaller range [0, ..., num_remaining_elems) after deleting
496// each element in `deletes` and shifting the remaining elements such that they
497// are contiguous starting at 0. The elements in the output point to the new
498// shifted index, or `kDeletedIndex` if the starting index was deleted.
499std::vector<int> IndexUpdateMap(const int size_before_delete,
500 absl::Span<const int> deletes) {
501 std::vector<int> result(size_before_delete, kUnsetIndex);
502 for (const int del : deletes) {
503 result[del] = kDeletedIndex;
504 }
505 int next_free = 0;
506 for (int& r : result) {
507 if (r != kDeletedIndex) {
508 r = next_free;
509 ++next_free;
510 }
511 CHECK_GT(r, kUnsetIndex);
512 }
513 return result;
514}
515
516absl::StatusOr<bool> GetIntAttrElementAsBool(Gurobi& gurobi,
517 const char* const name,
518 const int element) {
519 ASSIGN_OR_RETURN(const int value, gurobi.GetIntAttrElement(name, element));
520 const bool cast_value(value);
521 if (static_cast<int>(cast_value) != value) {
523 << "Gurobi unexpectedly returned non-boolean value for " << name
524 << ": " << value;
525 }
526 return cast_value;
527}
528
529class OrAccumulator {
530 public:
531 OrAccumulator() = default;
532 // Propagates any error in `update`. Otherwise, ORs the internal state with
533 // the value in `update`.
534 absl::Status Or(absl::StatusOr<bool> update) {
535 RETURN_IF_ERROR(update.status());
536 value_ |= *update;
537 return absl::OkStatus();
538 }
539 bool value() const { return value_; }
540
541 private:
542 bool value_ = false;
543};
544
545// Propagate any error in `maybe_value`. Otherwise, if `maybe_value` is set,
546// add (`id`, `maybe_value`) as an entry to `target`.
547absl::Status AddMapEntryIfPresent(
548 const int64_t map_id,
549 absl::StatusOr<std::optional<ModelSubsetProto::Bounds>> maybe_value,
550 google::protobuf::Map<int64_t, ModelSubsetProto::Bounds>& target) {
551 RETURN_IF_ERROR(maybe_value.status());
552 if (maybe_value->has_value()) {
553 target[map_id] = **std::move(maybe_value);
554 }
555 return absl::OkStatus();
556}
557
558// Propagate any error in `should_append`. Otherwise, if `should_append` is
559// true, append `id` to `target`.
560absl::Status AppendEntryIfTrue(
561 const int64_t id, absl::StatusOr<bool> should_append,
562 google::protobuf::RepeatedField<int64_t>& target) {
563 RETURN_IF_ERROR(should_append.status());
564 if (*should_append) {
565 target.Add(id);
566 }
567 return absl::OkStatus();
568}
569
570} // namespace
571
572GurobiSolver::GurobiSolver(std::unique_ptr<Gurobi> g_gurobi)
573 : gurobi_(std::move(g_gurobi)) {}
574
575GurobiSolver::GurobiModelElements
576GurobiSolver::LinearConstraintData::DependentElements() const {
577 GurobiModelElements elements;
578 CHECK_NE(constraint_index, kUnspecifiedConstraint);
579 elements.linear_constraints.push_back(constraint_index);
580 if (slack_index != kUnspecifiedIndex) {
581 elements.variables.push_back(slack_index);
582 }
583 return elements;
584}
585
586GurobiSolver::GurobiModelElements
587GurobiSolver::SecondOrderConeConstraintData::DependentElements() const {
588 const auto index_is_valid = [](const auto index) { return index >= 0; };
589 CHECK(absl::c_all_of(slack_variables, index_is_valid));
590 CHECK(absl::c_all_of(slack_constraints, index_is_valid));
591 CHECK_NE(constraint_index, kUnspecifiedConstraint);
592 GurobiModelElements elements{.variables = slack_variables,
593 .linear_constraints = slack_constraints};
594 elements.quadratic_constraints.push_back(constraint_index);
595 return elements;
596}
597
598GurobiSolver::GurobiModelElements
599GurobiSolver::SosConstraintData::DependentElements() const {
600 const auto index_is_valid = [](const auto index) { return index >= 0; };
601 CHECK(absl::c_all_of(slack_variables, index_is_valid));
602 CHECK(absl::c_all_of(slack_constraints, index_is_valid));
603 CHECK_NE(constraint_index, kUnspecifiedConstraint);
604 GurobiModelElements elements{.variables = slack_variables,
605 .linear_constraints = slack_constraints};
606 elements.sos_constraints.push_back(constraint_index);
607 return elements;
608}
609
610absl::StatusOr<TerminationProto> GurobiSolver::ConvertTerminationReason(
611 const int gurobi_status, const SolutionClaims solution_claims,
612 const double best_primal_bound, const double best_dual_bound) {
613 ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize());
614 switch (gurobi_status) {
615 case GRB_OPTIMAL:
616 // TODO(b/290359402): it appears Gurobi could return an infinite
617 // best_dual_bound (e.g in Qp/Qc/Socp/multi-obj tests). If so, we could
618 // improve the bound by using the target absolute and relative GAPs (and
619 // best_primal_bound).
620 return OptimalTerminationProto(best_primal_bound, best_dual_bound);
621 case GRB_INFEASIBLE:
623 is_maximize, solution_claims.dual_feasible_solution_exists
624 ? FEASIBILITY_STATUS_FEASIBLE
625 : FEASIBILITY_STATUS_UNDETERMINED);
626 case GRB_UNBOUNDED:
627 // GRB_UNBOUNDED does necessarily imply the primal is feasible
628 // https://www.gurobi.com/documentation/9.1/refman/optimization_status_codes.html
629 if (solution_claims.primal_feasible_solution_exists) {
630 return UnboundedTerminationProto(is_maximize);
631 }
633 is_maximize,
634 /*dual_feasibility_status=*/FEASIBILITY_STATUS_INFEASIBLE,
635 "Gurobi status GRB_UNBOUNDED");
636 case GRB_INF_OR_UNBD:
638 is_maximize,
639 /*dual_feasibility_status=*/FEASIBILITY_STATUS_UNDETERMINED,
640 "Gurobi status GRB_UNBOUNDED");
641 case GRB_CUTOFF:
642 return CutoffTerminationProto(is_maximize, "Gurobi status GRB_CUTOFF");
645 LIMIT_ITERATION, best_primal_bound, best_dual_bound,
646 solution_claims.dual_feasible_solution_exists);
647 case GRB_NODE_LIMIT:
649 LIMIT_NODE, best_primal_bound, best_dual_bound,
650 solution_claims.dual_feasible_solution_exists);
651 case GRB_TIME_LIMIT:
653 LIMIT_TIME, best_primal_bound, best_dual_bound,
654 solution_claims.dual_feasible_solution_exists);
657 LIMIT_SOLUTION, best_primal_bound, best_dual_bound,
658 solution_claims.dual_feasible_solution_exists);
659 case GRB_INTERRUPTED:
661 LIMIT_INTERRUPTED, best_primal_bound, best_dual_bound,
662 solution_claims.dual_feasible_solution_exists);
663 case GRB_NUMERIC:
664 return TerminateForReason(is_maximize,
665 TERMINATION_REASON_NUMERICAL_ERROR);
666 case GRB_SUBOPTIMAL: {
667 if (is_multi_objective_mode()) {
668 // Note: We state authoritatively that suboptimal means an objective
669 // time out, but we don't really know for sure that there are no other
670 // situations in multi-objective mode where this can occur.
672 LIMIT_TIME, best_primal_bound, best_dual_bound,
673 solution_claims.dual_feasible_solution_exists,
674 "Gurobi returned GRB_SUBOPTIMAL for a multi-objective model, which "
675 "indicates that one or more objectives hit their per-objective "
676 "time limit");
677 }
678 return TerminateForReason(is_maximize, TERMINATION_REASON_IMPRECISE);
679 }
681 // Note: maybe we should override
682 // solution_claims.primal_feasible_solution_exists to true or false
683 // depending on whether objective_limit and best_bound_limit triggered
684 // this. Not sure if it's possible to detect this though.
686 LIMIT_OBJECTIVE, best_primal_bound, best_dual_bound,
687 solution_claims.dual_feasible_solution_exists);
688 case GRB_LOADED:
689 return absl::InternalError(
690 "Error creating termination reason, unexpected gurobi status code "
691 "GRB_LOADED.");
692 case GRB_INPROGRESS:
693 return absl::InternalError(
694 "Error creating termination reason, unexpected gurobi status code "
695 "GRB_INPROGRESS.");
696 default:
697 return absl::InternalError(absl::StrCat(
698 "Missing Gurobi optimization status code case: ", gurobi_status));
699 }
700}
701
702absl::StatusOr<bool> GurobiSolver::IsMaximize() const {
703 ASSIGN_OR_RETURN(const int obj_sense,
704 gurobi_->GetIntAttr(GRB_INT_ATTR_MODELSENSE));
705 return obj_sense == GRB_MAXIMIZE;
706}
707
708absl::StatusOr<bool> GurobiSolver::IsMIP() const {
709 ASSIGN_OR_RETURN(const int is_mip, gurobi_->GetIntAttr(GRB_INT_ATTR_IS_MIP));
710 return static_cast<bool>(is_mip);
711}
712
713// TODO(b/204595455): Revisit logic when nonconvex QP support is decided upon
714absl::StatusOr<bool> GurobiSolver::IsQP() const {
715 ASSIGN_OR_RETURN(const int is_qp, gurobi_->GetIntAttr(GRB_INT_ATTR_IS_QP));
716 return static_cast<bool>(is_qp);
717}
718
719absl::StatusOr<bool> GurobiSolver::IsQCP() const {
720 ASSIGN_OR_RETURN(const int is_qcp, gurobi_->GetIntAttr(GRB_INT_ATTR_IS_QCP));
721 return static_cast<bool>(is_qcp);
722}
723
724// TODO(user): switch the use of this function to something closer to
725// GetGurobiDualRay()
726template <typename T>
727void GurobiSolver::GurobiVectorToSparseDoubleVector(
728 const absl::Span<const double> gurobi_values, const T& map,
729 SparseDoubleVectorProto& result,
730 const SparseVectorFilterProto& filter) const {
731 SparseVectorFilterPredicate predicate(filter);
732 for (auto [id, gurobi_data] : map) {
733 const double value = gurobi_values[get_model_index(gurobi_data)];
734 if (predicate.AcceptsAndUpdate(id, value)) {
735 result.add_ids(id);
736 result.add_values(value);
737 }
738 }
739}
740
741absl::Status GurobiSolver::SetGurobiBasis(const BasisProto& basis) {
742 std::vector<int> gurobi_variable_basis_status(num_gurobi_variables_);
743 for (const auto [id, value] : MakeView(basis.variable_status())) {
744 gurobi_variable_basis_status[variables_map_.at(id)] =
745 GrbVariableStatus(static_cast<BasisStatusProto>(value));
746 }
747
748 std::vector<int> gurobi_constraint_basis_status;
749 gurobi_constraint_basis_status.reserve(num_gurobi_lin_cons_);
750 for (const auto [id, value] : MakeView(basis.constraint_status())) {
751 const LinearConstraintData& constraint_data =
752 linear_constraints_map_.at(id);
753 // Non-ranged constraints
754 if (constraint_data.slack_index == kUnspecifiedIndex) {
755 if (value == BASIS_STATUS_BASIC) {
756 gurobi_constraint_basis_status.push_back(kGrbBasicConstraint);
757 } else {
758 gurobi_constraint_basis_status.push_back(kGrbNonBasicConstraint);
759 }
760 // Ranged constraints
761 } else if (value == BASIS_STATUS_BASIC) {
762 // Either constraint or MathOpt slack is basic, but not both (because
763 // columns for MathOpt slack and internal Gurobi slack are linearly
764 // dependent). We choose the MathOpt slack to be basic.
765 gurobi_variable_basis_status[constraint_data.slack_index] = GRB_BASIC;
766 gurobi_constraint_basis_status.push_back(kGrbNonBasicConstraint);
767 } else {
768 gurobi_variable_basis_status[constraint_data.slack_index] =
769 GrbVariableStatus(static_cast<BasisStatusProto>(value));
770 gurobi_constraint_basis_status.push_back(kGrbNonBasicConstraint);
771 }
772 }
773 RETURN_IF_ERROR(gurobi_->SetIntAttrArray(GRB_INT_ATTR_VBASIS,
774 gurobi_variable_basis_status));
775 RETURN_IF_ERROR(gurobi_->SetIntAttrArray(GRB_INT_ATTR_CBASIS,
776 gurobi_constraint_basis_status));
777 return absl::OkStatus();
778}
779
780absl::StatusOr<BasisProto> GurobiSolver::GetGurobiBasis() {
781 BasisProto basis;
783 const std::vector<int> gurobi_variable_basis_status,
784 gurobi_->GetIntAttrArray(GRB_INT_ATTR_VBASIS, num_gurobi_variables_));
785
786 for (auto [variable_id, gurobi_variable_index] : variables_map_) {
787 basis.mutable_variable_status()->add_ids(variable_id);
788 const BasisStatusProto variable_status = ConvertVariableStatus(
789 gurobi_variable_basis_status[gurobi_variable_index]);
790 if (variable_status == BASIS_STATUS_UNSPECIFIED) {
791 return absl::InternalError(
792 absl::StrCat("Invalid Gurobi variable basis status: ",
793 gurobi_variable_basis_status[gurobi_variable_index]));
794 }
795 basis.mutable_variable_status()->add_values(variable_status);
796 }
797
799 const std::vector<int> gurobi_constraint_basis_status,
800 gurobi_->GetIntAttrArray(GRB_INT_ATTR_CBASIS, num_gurobi_lin_cons_));
801 for (auto [constraint_id, gurobi_data] : linear_constraints_map_) {
802 basis.mutable_constraint_status()->add_ids(constraint_id);
803 const int gurobi_constraint_status =
804 gurobi_constraint_basis_status[gurobi_data.constraint_index];
805 if ((gurobi_constraint_status != kGrbBasicConstraint) &&
806 (gurobi_constraint_status != kGrbNonBasicConstraint)) {
807 return absl::InternalError(
808 absl::StrCat("Invalid Gurobi constraint basis status: ",
809 gurobi_constraint_status));
810 }
811 // linear_terms <= upper_bound
812 if (gurobi_data.lower_bound <= -GRB_INFINITY &&
813 gurobi_data.upper_bound < GRB_INFINITY) {
814 if (gurobi_constraint_status == kGrbBasicConstraint) {
815 basis.mutable_constraint_status()->add_values(BASIS_STATUS_BASIC);
816 } else {
817 basis.mutable_constraint_status()->add_values(
818 BASIS_STATUS_AT_UPPER_BOUND);
819 }
820 // linear_terms >= lower_bound
821 } else if (gurobi_data.lower_bound > -GRB_INFINITY &&
822 gurobi_data.upper_bound >= GRB_INFINITY) {
823 if (gurobi_constraint_status == kGrbBasicConstraint) {
824 basis.mutable_constraint_status()->add_values(BASIS_STATUS_BASIC);
825 } else {
826 basis.mutable_constraint_status()->add_values(
827 BASIS_STATUS_AT_LOWER_BOUND);
828 }
829 // linear_terms == xxxxx_bound
830 } else if (gurobi_data.lower_bound == gurobi_data.upper_bound) {
831 if (gurobi_constraint_status == kGrbBasicConstraint) {
832 basis.mutable_constraint_status()->add_values(BASIS_STATUS_BASIC);
833 } else {
834 // TODO(user): consider refining this to
835 // AT_LOWER_BOUND/AT_UPPER_BOUND using the sign of the dual variable.
836 basis.mutable_constraint_status()->add_values(BASIS_STATUS_FIXED_VALUE);
837 }
838 // linear_term - slack == 0 (ranged constraint)
839 } else {
840 const BasisStatusProto slack_status = ConvertVariableStatus(
841 gurobi_variable_basis_status[gurobi_data.slack_index]);
842 if (slack_status == BASIS_STATUS_UNSPECIFIED) {
843 return absl::InternalError(absl::StrCat(
844 "Invalid Gurobi slack variable basis status: ", slack_status));
845 }
846 if ((gurobi_constraint_status == kGrbBasicConstraint) ||
847 (slack_status == BASIS_STATUS_BASIC)) {
848 basis.mutable_constraint_status()->add_values(BASIS_STATUS_BASIC);
849 } else {
850 basis.mutable_constraint_status()->add_values(slack_status);
851 }
852 }
853 }
854 return basis;
855}
856
857absl::StatusOr<DualRayProto> GurobiSolver::GetGurobiDualRay(
858 const SparseVectorFilterProto& linear_constraints_filter,
859 const SparseVectorFilterProto& variables_filter, const bool is_maximize) {
860 // farkas_dual = lambda
861 ASSIGN_OR_RETURN(const std::vector<double> farkas_dual,
862 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_FARKASDUAL,
863 num_gurobi_lin_cons_));
864
865 DualRayProto dual_ray;
866
867 // Compute y = -lambda
868 {
869 SparseVectorFilterPredicate predicate(linear_constraints_filter);
870 for (auto [constraint_id, gurobi_data] : linear_constraints_map_) {
871 // constraint_dual_value = y[constraint_id]
872 const double value = -farkas_dual[gurobi_data.constraint_index];
873 if (predicate.AcceptsAndUpdate(constraint_id, value)) {
874 dual_ray.mutable_dual_values()->add_ids(constraint_id);
875 if (is_maximize) {
876 dual_ray.mutable_dual_values()->add_values(-value);
877 } else {
878 dual_ray.mutable_dual_values()->add_values(value);
879 }
880 }
881 }
882 }
883
884 // Compute r = \bar{a} = A^T lambda
885 {
886 SparseVectorFilterPredicate predicate(variables_filter);
887 for (auto [var_id, gurobi_variable_index] : variables_map_) {
888 // reduced_cost_value = r[gurobi_variable_index]
889 // = \bar{a}[gurobi_variable_index]
890 double reduced_cost_value = 0.0;
891 ASSIGN_OR_RETURN(Gurobi::SparseMat column,
892 gurobi_->GetVars(gurobi_variable_index, 1));
893 for (int i = 0; i < column.inds.size(); ++i) {
894 reduced_cost_value += farkas_dual[column.inds[i]] * column.vals[i];
895 }
896 if (predicate.AcceptsAndUpdate(var_id, reduced_cost_value)) {
897 dual_ray.mutable_reduced_costs()->add_ids(var_id);
898 if (is_maximize) {
899 dual_ray.mutable_reduced_costs()->add_values(-reduced_cost_value);
900 } else {
901 dual_ray.mutable_reduced_costs()->add_values(reduced_cost_value);
902 }
903 }
904 }
905 }
906 return dual_ray;
907}
908
909absl::StatusOr<SolveResultProto> GurobiSolver::ExtractSolveResultProto(
910 const absl::Time start, const ModelSolveParametersProto& model_parameters) {
911 SolveResultProto result;
912 ASSIGN_OR_RETURN(const int grb_termination,
913 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
914 SolutionClaims solution_claims;
915 ASSIGN_OR_RETURN(SolutionsAndClaims solution_and_claims,
916 GetSolutions(model_parameters));
917 if (grb_termination == GRB_CUTOFF) {
918 // Cutoff will not be triggered by bounds e.g. for LP dual feasible
919 // solutions. In particular, if the problem is both primal and dual
920 // infeasible, we will not get a bound and should not be returning CUTOFF.
921 //
922 // TODO(b/272268188): test that this has no bad interactions with primal +
923 // dual infeasible problems.
924 solution_and_claims.solutions.clear();
925 solution_and_claims.solution_claims.primal_feasible_solution_exists = false;
926 solution_and_claims.solution_claims.dual_feasible_solution_exists = true;
927 }
928 ASSIGN_OR_RETURN(const double best_primal_bound,
929 GetBestPrimalBound(solution_and_claims.solutions));
930 ASSIGN_OR_RETURN(const double best_dual_bound,
931 GetBestDualBound(solution_and_claims.solutions));
932 solution_claims = solution_and_claims.solution_claims;
933
934 // TODO(b/195295177): Add tests for rays in unbounded MIPs
935 RETURN_IF_ERROR(FillRays(model_parameters, solution_claims, result));
936
937 for (SolutionProto& solution : solution_and_claims.solutions) {
938 *result.add_solutions() = std::move(solution);
939 }
940
942 *result.mutable_termination(),
943 ConvertTerminationReason(grb_termination, solution_claims,
944 best_primal_bound, best_dual_bound));
945
946 ASSIGN_OR_RETURN(*result.mutable_solve_stats(), GetSolveStats(start));
947 return result;
948}
949
950absl::StatusOr<bool> GurobiSolver::AnyElementInIIS(
951 const GurobiModelElements& grb_elements) {
952 OrAccumulator any_in_iis;
953 for (const GurobiVariableIndex grb_index : grb_elements.variables) {
954 RETURN_IF_ERROR(any_in_iis.Or(VariableInIIS(grb_index)));
955 }
956 for (const GurobiLinearConstraintIndex grb_index :
957 grb_elements.linear_constraints) {
958 RETURN_IF_ERROR(any_in_iis.Or(
959 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_CONSTR, grb_index)));
960 }
961 for (const GurobiQuadraticConstraintIndex grb_index :
962 grb_elements.quadratic_constraints) {
963 RETURN_IF_ERROR(any_in_iis.Or(GetIntAttrElementAsBool(
964 *gurobi_, GRB_INT_ATTR_IIS_QCONSTR, grb_index)));
965 }
966 for (const GurobiSosConstraintIndex grb_index :
967 grb_elements.sos_constraints) {
968 RETURN_IF_ERROR(any_in_iis.Or(
969 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_SOS, grb_index)));
970 }
971 for (const GurobiGeneralConstraintIndex grb_index :
972 grb_elements.general_constraints) {
973 RETURN_IF_ERROR(any_in_iis.Or(GetIntAttrElementAsBool(
974 *gurobi_, GRB_INT_ATTR_IIS_GENCONSTR, grb_index)));
975 }
976 return any_in_iis.value();
977}
978
979// If set, the returned value will have at least one of .lower and/or .upper set
980// to true.
981absl::StatusOr<std::optional<ModelSubsetProto::Bounds>>
982GurobiSolver::VariableBoundsInIIS(const GurobiVariableIndex grb_index) {
984 const bool var_lb_is_in_iis,
985 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_LB, grb_index));
987 const bool var_ub_is_in_iis,
988 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_UB, grb_index));
989 if (!var_lb_is_in_iis && !var_ub_is_in_iis) {
990 return std::nullopt;
991 }
992 ModelSubsetProto::Bounds bounds;
993 bounds.set_lower(var_lb_is_in_iis);
994 bounds.set_upper(var_ub_is_in_iis);
995 return bounds;
996}
997
998absl::StatusOr<bool> GurobiSolver::VariableInIIS(
999 const GurobiVariableIndex grb_index) {
1000 ASSIGN_OR_RETURN(const std::optional<ModelSubsetProto::Bounds> bounds,
1001 VariableBoundsInIIS(grb_index));
1002 return bounds.has_value();
1003}
1004
1005absl::StatusOr<std::optional<ModelSubsetProto::Bounds>>
1006GurobiSolver::LinearConstraintInIIS(const LinearConstraintData& grb_data) {
1007 // Attributing which part of an equality/ranged constraint is part of the
1008 // IIS can be tricky. So, we take a conservative approach: if anything
1009 // associated with this linear constraint is part of the IIS (including
1010 // slack variable/constraints), then we mark any finite bound as being part
1011 // of the IIS.
1012 ASSIGN_OR_RETURN(const bool constr_in_iis,
1013 AnyElementInIIS(grb_data.DependentElements()));
1014 if (!constr_in_iis) {
1015 return std::nullopt;
1016 }
1017 ModelSubsetProto::Bounds result;
1018 result.set_lower(grb_data.lower_bound != -kInf);
1019 result.set_upper(grb_data.upper_bound != kInf);
1020 return result;
1021}
1022
1023absl::StatusOr<std::optional<ModelSubsetProto::Bounds>>
1024GurobiSolver::QuadraticConstraintInIIS(
1025 const GurobiQuadraticConstraintIndex grb_index) {
1027 const bool constr_in_iis,
1028 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_QCONSTR, grb_index));
1029 if (!constr_in_iis) {
1030 return std::nullopt;
1031 }
1033 const char constr_sense,
1034 gurobi_->GetCharAttrElement(GRB_CHAR_ATTR_QCSENSE, grb_index));
1035 ModelSubsetProto::Bounds result;
1036 result.set_lower(constr_sense == GRB_EQUAL ||
1037 constr_sense == GRB_GREATER_EQUAL);
1038 result.set_upper(constr_sense == GRB_EQUAL || constr_sense == GRB_LESS_EQUAL);
1039 return result;
1040}
1041
1042absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
1043GurobiSolver::ExtractComputeInfeasibleSubsystemResultProto(
1044 const bool proven_infeasible) {
1045 ComputeInfeasibleSubsystemResultProto result;
1046 if (!proven_infeasible) {
1047 result.set_feasibility(FEASIBILITY_STATUS_UNDETERMINED);
1048 return result;
1049 }
1050 result.set_feasibility(FEASIBILITY_STATUS_INFEASIBLE);
1051 {
1052 ASSIGN_OR_RETURN(const bool is_minimal,
1053 gurobi_->GetIntAttr(GRB_INT_ATTR_IIS_MINIMAL));
1054 result.set_is_minimal(is_minimal);
1055 }
1056 for (const auto [id, grb_index] : variables_map_) {
1057 RETURN_IF_ERROR(AddMapEntryIfPresent(
1058 id, /*maybe_value=*/VariableBoundsInIIS(grb_index),
1059 *result.mutable_infeasible_subsystem()->mutable_variable_bounds()));
1060 {
1061 // Gurobi does not report integrality in its IIS, so we mark any integer
1062 // variable in the model as being involved.
1064 const char var_type,
1065 gurobi_->GetCharAttrElement(GRB_CHAR_ATTR_VTYPE, grb_index));
1066 if (var_type == GRB_BINARY || var_type == GRB_INTEGER) {
1067 result.mutable_infeasible_subsystem()->add_variable_integrality(
1068 grb_index);
1069 }
1070 }
1071 }
1072
1073 for (const auto [id, grb_data] : linear_constraints_map_) {
1074 RETURN_IF_ERROR(AddMapEntryIfPresent(
1075 id, /*maybe_value=*/LinearConstraintInIIS(grb_data),
1076 *result.mutable_infeasible_subsystem()->mutable_linear_constraints()));
1077 }
1078
1079 for (const auto [id, grb_index] : quadratic_constraints_map_) {
1080 RETURN_IF_ERROR(AddMapEntryIfPresent(
1081 id, /*maybe_value=*/QuadraticConstraintInIIS(grb_index),
1082 *result.mutable_infeasible_subsystem()
1083 ->mutable_quadratic_constraints()));
1084 }
1085
1086 for (const auto& [id, grb_data] : soc_constraints_map_) {
1087 RETURN_IF_ERROR(AppendEntryIfTrue(
1088 id, /*should_append=*/AnyElementInIIS(grb_data.DependentElements()),
1089 *result.mutable_infeasible_subsystem()
1090 ->mutable_second_order_cone_constraints()));
1091 }
1092
1093 for (const auto& [id, grb_data] : sos1_constraints_map_) {
1094 RETURN_IF_ERROR(AppendEntryIfTrue(
1095 id, /*should_append=*/AnyElementInIIS(grb_data.DependentElements()),
1096 *result.mutable_infeasible_subsystem()->mutable_sos1_constraints()));
1097 }
1098
1099 for (const auto& [id, grb_data] : sos2_constraints_map_) {
1100 RETURN_IF_ERROR(AppendEntryIfTrue(
1101 id, /*should_append=*/AnyElementInIIS(grb_data.DependentElements()),
1102 *result.mutable_infeasible_subsystem()->mutable_sos2_constraints()));
1103 }
1104
1105 for (const auto& [id, maybe_grb_data] : indicator_constraints_map_) {
1106 if (!maybe_grb_data.has_value()) {
1107 continue;
1108 }
1109 RETURN_IF_ERROR(AppendEntryIfTrue(
1110 id,
1111 /*should_append=*/
1112 GetIntAttrElementAsBool(*gurobi_, GRB_INT_ATTR_IIS_GENCONSTR,
1113 maybe_grb_data->constraint_index),
1114 *result.mutable_infeasible_subsystem()
1115 ->mutable_indicator_constraints()));
1116 }
1117
1118 // Since each xxx_map_ is in insertion order, we do not need to worry about
1119 // sorting the repeated fields in `result`.
1120 return result;
1121}
1122
1123absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetSolutions(
1124 const ModelSolveParametersProto& model_parameters) {
1125 // Note that all multi-objective models will have `IsMip()` return true.
1126 ASSIGN_OR_RETURN(const bool is_mip, IsMIP());
1127 ASSIGN_OR_RETURN(const bool is_qp, IsQP());
1128 ASSIGN_OR_RETURN(const bool is_qcp, IsQCP());
1129
1130 if (is_mip) {
1131 return GetMipSolutions(model_parameters);
1132 } else if (is_qcp) {
1133 return GetQcpSolution(model_parameters);
1134 } else if (is_qp) {
1135 return GetQpSolution(model_parameters);
1136 } else {
1137 return GetLpSolution(model_parameters);
1138 }
1139}
1140
1141// TODO: b/365762174 - Remove logging and clamping below when Gurobi fixes its
1142// behavior upstream (and we migrate onto that version).
1143absl::StatusOr<SolveStatsProto> GurobiSolver::GetSolveStats(
1144 const absl::Time start) const {
1145 SolveStatsProto solve_stats;
1146
1147 CHECK_OK(util_time::EncodeGoogleApiProto(absl::Now() - start,
1148 solve_stats.mutable_solve_time()));
1149
1150 if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_ITERCOUNT)) {
1151 ASSIGN_OR_RETURN(const double simplex_iters_double,
1152 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_ITERCOUNT));
1153 ASSIGN_OR_RETURN(const int64_t simplex_iters,
1154 SafeInt64FromDouble(simplex_iters_double));
1155 LOG_IF(ERROR, simplex_iters < 0)
1156 << "Expected GRB_DBL_ATTR_ITERCOUNT to be non-negative, got: "
1157 << simplex_iters << "; clamping to 0";
1158 solve_stats.set_simplex_iterations(std::max(int64_t{0}, simplex_iters));
1159 }
1160
1161 if (gurobi_->IsAttrAvailable(GRB_INT_ATTR_BARITERCOUNT)) {
1162 ASSIGN_OR_RETURN(const int barrier_iters,
1163 gurobi_->GetIntAttr(GRB_INT_ATTR_BARITERCOUNT));
1164 LOG_IF(ERROR, barrier_iters < 0)
1165 << "Expected GRB_INT_ATTR_BARITERCOUNT to be non-negative, got: "
1166 << barrier_iters << "; clamping to 0";
1167 solve_stats.set_barrier_iterations(std::max(0, barrier_iters));
1168 }
1169
1170 if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_NODECOUNT)) {
1171 ASSIGN_OR_RETURN(const double nodes_double,
1172 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_NODECOUNT));
1173 ASSIGN_OR_RETURN(const int64_t nodes, SafeInt64FromDouble(nodes_double));
1174 LOG_IF(ERROR, nodes < 0)
1175 << "Expected GRB_DBL_ATTR_NODECOUNT to be non-negative, got: " << nodes
1176 << "; clamping to 0";
1177 solve_stats.set_node_count(std::max(int64_t{0}, nodes));
1178 }
1179 return solve_stats;
1180}
1181
1182absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetMipSolutions(
1183 const ModelSolveParametersProto& model_parameters) {
1184 int num_solutions = 0;
1185 if (gurobi_->IsAttrAvailable(GRB_INT_ATTR_SOLCOUNT)) {
1186 ASSIGN_OR_RETURN(num_solutions, gurobi_->GetIntAttr(GRB_INT_ATTR_SOLCOUNT));
1187 }
1188 std::vector<SolutionProto> solutions;
1189 solutions.reserve(num_solutions);
1190
1191 for (int i = 0; i < num_solutions; ++i) {
1192 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_SOLUTIONNUMBER, i));
1193
1194 PrimalSolutionProto primal_solution;
1195 ASSIGN_OR_RETURN(const double sol_val,
1196 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_POOLOBJVAL));
1197 primal_solution.set_objective_value(sol_val);
1198 if (is_multi_objective_mode()) {
1199 for (const auto [id, grb_index] : multi_objectives_map_) {
1200 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_OBJNUMBER, grb_index));
1201 ASSIGN_OR_RETURN(const double obj_val,
1202 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_OBJNVAL));
1203 // If unset, this is the primary objective. We have already queried its
1204 // value via PoolObjVal above.
1205 if (id.has_value()) {
1206 (*primal_solution.mutable_auxiliary_objective_values())[*id] =
1207 obj_val;
1208 }
1209 }
1210 }
1211 // Gurobi v9 provides a feasibility status for the instance as a whole but
1212 // not for each solution, and pool entries may be infeasible. To be
1213 // conservative, we only label the first ("best") solution as primal
1214 // feasible.
1215 primal_solution.set_feasibility_status(
1216 i == 0 ? SOLUTION_STATUS_FEASIBLE : SOLUTION_STATUS_UNDETERMINED);
1218 const std::vector<double> grb_var_values,
1219 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_XN, num_gurobi_variables_));
1220 GurobiVectorToSparseDoubleVector(grb_var_values, variables_map_,
1221 *primal_solution.mutable_variable_values(),
1222 model_parameters.variable_values_filter());
1223 *solutions.emplace_back(SolutionProto()).mutable_primal_solution() =
1224 std::move(primal_solution);
1225 }
1226
1227 ASSIGN_OR_RETURN(const int grb_termination,
1228 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
1229 // Set solution claims
1230 ASSIGN_OR_RETURN(const double best_dual_bound, GetGurobiBestDualBound());
1231 // Note: here the existence of a dual solution refers to a dual solution to
1232 // some convex relaxation of the MIP. This convex relaxation can likely be
1233 // interpreted as an LP between the LP relaxation of the MIP and the convex
1234 // hull of feasible solutions of the MIP. However, here we only use the fact
1235 // that best_dual_bound being finite implies the existence of the trivial
1236 // convex relaxation given by (assuming a minimization problem with objective
1237 // function c^T x): min{c^T x : c^T x >= best_dual_bound}.
1238 //
1239 // If this is a multi-objective model, Gurobi v10 does not expose ObjBound.
1240 // Instead, we fake its existence for optimal solves only.
1241 // By convention infeasible MIPs are always dual feasible.
1242 const SolutionClaims solution_claims = {
1243 .primal_feasible_solution_exists = num_solutions > 0,
1244 .dual_feasible_solution_exists =
1245 std::isfinite(best_dual_bound) || grb_termination == GRB_INFEASIBLE ||
1246 (is_multi_objective_mode() && grb_termination == GRB_OPTIMAL)};
1247
1248 // Check consistency of solutions, bounds and statuses.
1249 if (grb_termination == GRB_OPTIMAL && num_solutions == 0) {
1250 return absl::InternalError(
1251 "GRB_INT_ATTR_STATUS == GRB_OPTIMAL, but solution pool is empty.");
1252 }
1253 // As set above, in multi-objective mode the dual bound is not informative and
1254 // it will not pass this validation.
1255 if (!is_multi_objective_mode() && grb_termination == GRB_OPTIMAL &&
1256 !std::isfinite(best_dual_bound)) {
1257 return absl::InternalError(
1258 "GRB_INT_ATTR_STATUS == GRB_OPTIMAL, but GRB_DBL_ATTR_OBJBOUND is "
1259 "unavailable or infinite.");
1260 }
1261
1262 return SolutionsAndClaims{.solutions = std::move(solutions),
1263 .solution_claims = solution_claims};
1264}
1265
1266absl::StatusOr<GurobiSolver::SolutionAndClaim<PrimalSolutionProto>>
1267GurobiSolver::GetConvexPrimalSolutionIfAvailable(
1268 const ModelSolveParametersProto& model_parameters) {
1269 if (!gurobi_->IsAttrAvailable(GRB_DBL_ATTR_X)) {
1270 return SolutionAndClaim<PrimalSolutionProto>{
1271 .solution = std::nullopt, .feasible_solution_exists = false};
1272 }
1273 ASSIGN_OR_RETURN(const int grb_termination,
1274 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
1275
1276 // Get primal solutions if available.
1278 const std::vector<double> grb_var_values,
1279 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_X, num_gurobi_variables_));
1280
1281 PrimalSolutionProto primal_solution;
1282 // The objective value may be missing for primal feasible solutions for
1283 // unbounded problems.
1284 // TODO(b/195295177): for GRB_ITERATION_LIMIT an objective value of 0.0 is
1285 // returned which breaks LpIncompleteSolveTest.PrimalSimplexAlgorithm. Explore
1286 // more and make simple example to file a bug.
1287 if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_OBJVAL) &&
1288 grb_termination != GRB_ITERATION_LIMIT) {
1289 ASSIGN_OR_RETURN(const double sol_val,
1290 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_OBJVAL));
1291 primal_solution.set_objective_value(sol_val);
1292 } else {
1293 double objective_value = 0.0;
1295 const std::vector<double> linear_obj_coefs,
1296 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_OBJ, num_gurobi_variables_));
1297 for (int i = 0; i < num_gurobi_variables_; ++i) {
1298 objective_value += linear_obj_coefs[i] * grb_var_values[i];
1299 }
1300 primal_solution.set_objective_value(objective_value);
1301 }
1302
1303 primal_solution.set_feasibility_status(SOLUTION_STATUS_UNDETERMINED);
1304 if (grb_termination == GRB_OPTIMAL) {
1305 primal_solution.set_feasibility_status(SOLUTION_STATUS_FEASIBLE);
1306 } else if (grb_termination == GRB_INFEASIBLE) {
1307 primal_solution.set_feasibility_status(SOLUTION_STATUS_INFEASIBLE);
1308 } else if (PrimalSolutionQualityAvailable()) {
1309 ASSIGN_OR_RETURN(const double solution_quality, GetPrimalSolutionQuality());
1310 ASSIGN_OR_RETURN(const double tolerance,
1311 gurobi_->GetDoubleParam(GRB_DBL_PAR_FEASIBILITYTOL));
1312 if (solution_quality <= tolerance) {
1313 primal_solution.set_feasibility_status(SOLUTION_STATUS_FEASIBLE);
1314 } else {
1315 primal_solution.set_feasibility_status(SOLUTION_STATUS_INFEASIBLE);
1316 }
1317 }
1318
1319 GurobiVectorToSparseDoubleVector(grb_var_values, variables_map_,
1320 *primal_solution.mutable_variable_values(),
1321 model_parameters.variable_values_filter());
1322 const bool primal_feasible_solution_exists =
1323 (primal_solution.feasibility_status() == SOLUTION_STATUS_FEASIBLE);
1324 return SolutionAndClaim<PrimalSolutionProto>{
1325 .solution = std::move(primal_solution),
1326 .feasible_solution_exists = primal_feasible_solution_exists};
1327}
1328
1329bool GurobiSolver::PrimalSolutionQualityAvailable() const {
1330 return gurobi_->IsAttrAvailable(GRB_DBL_ATTR_CONSTR_RESIDUAL) &&
1331 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_CONSTR_VIO) &&
1332 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_BOUND_VIO) &&
1333 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_CONSTR_SRESIDUAL) &&
1334 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_CONSTR_SVIO) &&
1335 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_BOUND_SVIO);
1336}
1337
1338absl::StatusOr<double> GurobiSolver::GetPrimalSolutionQuality() const {
1339 ASSIGN_OR_RETURN(const double constraint_residual,
1340 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_CONSTR_RESIDUAL));
1341 ASSIGN_OR_RETURN(const double constraint_violation,
1342 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_CONSTR_VIO));
1343 ASSIGN_OR_RETURN(const double bound_violation,
1344 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_BOUND_VIO));
1345 ASSIGN_OR_RETURN(const double constraint_scaled_residual,
1346 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_CONSTR_SRESIDUAL));
1347 ASSIGN_OR_RETURN(const double constraint_scaled_violation,
1348 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_CONSTR_SVIO));
1349 ASSIGN_OR_RETURN(const double bound_scaled_violation,
1350 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_BOUND_SVIO));
1351 return std::max({constraint_residual, constraint_violation, bound_violation,
1352 constraint_scaled_residual, constraint_scaled_violation,
1353 bound_scaled_violation});
1354}
1355
1356absl::StatusOr<double> GurobiSolver::GetBestPrimalBound(
1357 absl::Span<const SolutionProto> solutions) const {
1358 // We avoid using GRB_DBL_ATTR_OBJVAL because it may be incorrect on early
1359 // termination and for infeasible solutions (as of Gurobi 9.0.1).
1360 // Note that for (primal) unbounded problems the primal_bound is correctly
1361 // filled by ConvertTerminationReason() possibly overriding the value
1362 // returned by GetBestPrimalBound().
1363 ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize());
1364 double best_objective_value = is_maximize ? -kInf : kInf;
1365 for (const SolutionProto& solution : solutions) {
1366 if (solution.has_primal_solution() &&
1367 solution.primal_solution().feasibility_status() ==
1368 SOLUTION_STATUS_FEASIBLE) {
1369 const double sol_val = solution.primal_solution().objective_value();
1370 best_objective_value = is_maximize
1371 ? std::max(best_objective_value, sol_val)
1372 : std::min(best_objective_value, sol_val);
1373 }
1374 }
1375 return best_objective_value;
1376}
1377
1378absl::StatusOr<double> GurobiSolver::GetBestDualBound(
1379 absl::Span<const SolutionProto> solutions) const {
1380 ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize());
1381 // GetGurobiBestDualBound() returns the correct bound for problems without
1382 // dual solutions (e.g. MIP).
1383 ASSIGN_OR_RETURN(double best_dual_bound, GetGurobiBestDualBound());
1384 // However, GetGurobiBestDualBound() may be incorrect for QP problems.
1385 for (const SolutionProto& solution : solutions) {
1386 if (solution.has_dual_solution() &&
1387 solution.dual_solution().feasibility_status() ==
1388 SOLUTION_STATUS_FEASIBLE &&
1389 solution.dual_solution().has_objective_value()) {
1390 const double sol_val = solution.dual_solution().objective_value();
1391 best_dual_bound = is_maximize ? std::min(best_dual_bound, sol_val)
1392 : std::max(best_dual_bound, sol_val);
1393 }
1394 }
1395 return best_dual_bound;
1396}
1397
1398absl::StatusOr<double> GurobiSolver::GetGurobiBestDualBound() const {
1399 // As of v9.0.2, on multi objective models Gurobi incorrectly reports that
1400 // ObjBound is available. We work around this by adding a check if we are in
1401 // multi objective mode.
1402 if (gurobi_->IsAttrAvailable(GRB_DBL_ATTR_OBJBOUND) &&
1403 !is_multi_objective_mode()) {
1404 ASSIGN_OR_RETURN(const double obj_bound,
1405 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_OBJBOUND));
1406 // Note: Unbounded models return GRB_DBL_ATTR_OBJBOUND = GRB_INFINITY so
1407 // the conversion of +/-GRB_INFINITY to +/-kInf is needed and consistent.
1408 if (std::abs(obj_bound) < GRB_INFINITY) {
1409 return obj_bound;
1410 }
1411 }
1412 ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize());
1413 return is_maximize ? kInf : -kInf;
1414}
1415
1416absl::StatusOr<std::optional<BasisProto>> GurobiSolver::GetBasisIfAvailable() {
1417 if (gurobi_->IsAttrAvailable(GRB_INT_ATTR_VBASIS) &&
1418 gurobi_->IsAttrAvailable(GRB_INT_ATTR_CBASIS)) {
1419 ASSIGN_OR_RETURN(BasisProto basis, GetGurobiBasis());
1420 ASSIGN_OR_RETURN(const int grb_termination,
1421 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
1422 basis.set_basic_dual_feasibility(SOLUTION_STATUS_UNDETERMINED);
1423 if (grb_termination == GRB_OPTIMAL) {
1424 basis.set_basic_dual_feasibility(SOLUTION_STATUS_FEASIBLE);
1425 } else if (grb_termination == GRB_UNBOUNDED) {
1426 basis.set_basic_dual_feasibility(SOLUTION_STATUS_INFEASIBLE);
1427 }
1428 // TODO(b/195295177): double check if the move is needed
1429 return std::move(basis);
1430 }
1431 return std::nullopt;
1432}
1433
1434absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetLpSolution(
1435 const ModelSolveParametersProto& model_parameters) {
1436 ASSIGN_OR_RETURN(auto primal_solution_and_claim,
1437 GetConvexPrimalSolutionIfAvailable(model_parameters));
1438 ASSIGN_OR_RETURN(auto dual_solution_and_claim,
1439 GetConvexDualSolutionIfAvailable(model_parameters));
1440 ASSIGN_OR_RETURN(auto basis, GetBasisIfAvailable());
1441 const SolutionClaims solution_claims = {
1442 .primal_feasible_solution_exists =
1443 primal_solution_and_claim.feasible_solution_exists,
1444 .dual_feasible_solution_exists =
1445 dual_solution_and_claim.feasible_solution_exists};
1446
1447 if (!primal_solution_and_claim.solution.has_value() &&
1448 !dual_solution_and_claim.solution.has_value() && !basis.has_value()) {
1449 return SolutionsAndClaims{.solution_claims = solution_claims};
1450 }
1451 SolutionsAndClaims solution_and_claims{.solution_claims = solution_claims};
1452 SolutionProto& solution =
1453 solution_and_claims.solutions.emplace_back(SolutionProto());
1454 if (primal_solution_and_claim.solution.has_value()) {
1455 *solution.mutable_primal_solution() =
1456 std::move(*primal_solution_and_claim.solution);
1457 }
1458 if (dual_solution_and_claim.solution.has_value()) {
1459 *solution.mutable_dual_solution() =
1460 std::move(*dual_solution_and_claim.solution);
1461 }
1462 if (basis.has_value()) {
1463 *solution.mutable_basis() = std::move(*basis);
1464 }
1465 return solution_and_claims;
1466}
1467
1468absl::StatusOr<GurobiSolver::SolutionAndClaim<DualSolutionProto>>
1469GurobiSolver::GetConvexDualSolutionIfAvailable(
1470 const ModelSolveParametersProto& model_parameters) {
1471 if (!gurobi_->IsAttrAvailable(GRB_DBL_ATTR_PI) ||
1472 !gurobi_->IsAttrAvailable(GRB_DBL_ATTR_RC)) {
1473 return SolutionAndClaim<DualSolutionProto>{
1474 .solution = std::nullopt, .feasible_solution_exists = false};
1475 }
1476
1477 // Note that we can ignore the reduced costs of the slack variables for
1478 // ranged constraints.
1479 DualSolutionProto dual_solution;
1480 bool dual_feasible_solution_exists = false;
1482 const std::vector<double> grb_constraint_duals,
1483 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_PI, num_gurobi_lin_cons_));
1484 GurobiVectorToSparseDoubleVector(grb_constraint_duals,
1485 linear_constraints_map_,
1486 *dual_solution.mutable_dual_values(),
1487 model_parameters.dual_values_filter());
1488
1490 const std::vector<double> grb_reduced_cost_values,
1491 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_RC, num_gurobi_variables_));
1492 GurobiVectorToSparseDoubleVector(grb_reduced_cost_values, variables_map_,
1493 *dual_solution.mutable_reduced_costs(),
1494 model_parameters.reduced_costs_filter());
1495
1496 if (!quadratic_constraints_map_.empty() &&
1497 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_QCPI)) {
1499 const std::vector<double> grb_quad_constraint_duals,
1500 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_QCPI, num_gurobi_quad_cons_));
1501 GurobiVectorToSparseDoubleVector(
1502 grb_quad_constraint_duals, quadratic_constraints_map_,
1503 *dual_solution.mutable_quadratic_dual_values(),
1504 model_parameters.quadratic_dual_values_filter());
1505 }
1506
1507 ASSIGN_OR_RETURN(const int grb_termination,
1508 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
1509 if (grb_termination == GRB_OPTIMAL &&
1510 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_OBJVAL)) {
1511 ASSIGN_OR_RETURN(const double obj_val,
1512 gurobi_->GetDoubleAttr(GRB_DBL_ATTR_OBJVAL));
1513 dual_solution.set_objective_value(obj_val);
1514 }
1515
1516 dual_solution.set_feasibility_status(SOLUTION_STATUS_UNDETERMINED);
1517 if (grb_termination == GRB_OPTIMAL) {
1518 dual_solution.set_feasibility_status(SOLUTION_STATUS_FEASIBLE);
1519 dual_feasible_solution_exists = true;
1520 } else if (grb_termination == GRB_UNBOUNDED) {
1521 dual_solution.set_feasibility_status(SOLUTION_STATUS_INFEASIBLE);
1522 }
1523 // TODO(b/290359402): We could use gurobi's dual solution quality measures
1524 // for further upgrade the dual feasibility but it likely is only useful
1525 // for phase II of dual simplex because:
1526 // * the quality measures seem to evaluate if the basis is dual feasible
1527 // so for primal simplex we would not improve over checking
1528 // GRB_OPTIMAL.
1529 // * for phase I dual simplex we cannot rely on the quality measures.
1530 // We could also use finiteness of GRB_DBL_ATTR_OBJBOUND to deduce dual
1531 // feasibility.
1532
1533 // Note: GRB_DBL_ATTR_OBJBOUND can sometimes provide the objective value of a
1534 // sub-optimal dual feasible solution.
1535 // Here we only use it to possibly update dual_feasible_solution_exists.
1536 ASSIGN_OR_RETURN(const double best_dual_bound, GetGurobiBestDualBound());
1537 if (dual_feasible_solution_exists || std::isfinite(best_dual_bound)) {
1538 dual_feasible_solution_exists = true;
1539 } else if (grb_termination == GRB_OPTIMAL) {
1540 return absl::InternalError(
1541 "GRB_INT_ATTR_STATUS == GRB_OPTIMAL, but GRB_DBL_ATTR_OBJBOUND is "
1542 "unavailable or infinite, and no dual feasible solution is returned");
1543 }
1544 return SolutionAndClaim<DualSolutionProto>{
1545 .solution = std::move(dual_solution),
1546 .feasible_solution_exists = dual_feasible_solution_exists};
1547}
1548
1549absl::Status GurobiSolver::FillRays(
1550 const ModelSolveParametersProto& model_parameters,
1551 const SolutionClaims solution_claims, SolveResultProto& result) {
1552 ASSIGN_OR_RETURN(const bool is_maximize, IsMaximize());
1553 // GRB_DBL_ATTR_UNBDRAY is sometimes incorrectly available for problems
1554 // without variables. We also give priority to the conclusions obtained from
1555 // dual solutions or bounds.
1556 if (!solution_claims.dual_feasible_solution_exists &&
1557 num_gurobi_variables_ > 0 &&
1558 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_UNBDRAY)) {
1559 ASSIGN_OR_RETURN(const std::vector<double> grb_ray_var_values,
1560 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_UNBDRAY,
1561 num_gurobi_variables_));
1562 PrimalRayProto* const primal_ray = result.add_primal_rays();
1563 GurobiVectorToSparseDoubleVector(grb_ray_var_values, variables_map_,
1564 *primal_ray->mutable_variable_values(),
1565 model_parameters.variable_values_filter());
1566 }
1567 // GRB_DBL_ATTR_FARKASDUAL is sometimes incorrectly available for problems
1568 // without constraints. We also give priority to the conclusions obtained from
1569 // primal solutions.
1570 if (!solution_claims.primal_feasible_solution_exists &&
1571 num_gurobi_lin_cons_ > 0 &&
1572 gurobi_->IsAttrAvailable(GRB_DBL_ATTR_FARKASDUAL)) {
1574 DualRayProto dual_ray,
1575 GetGurobiDualRay(model_parameters.dual_values_filter(),
1576 model_parameters.reduced_costs_filter(), is_maximize));
1577 result.mutable_dual_rays()->Add(std::move(dual_ray));
1578 }
1579 return absl::OkStatus();
1580}
1581
1582absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetQpSolution(
1583 const ModelSolveParametersProto& model_parameters) {
1584 ASSIGN_OR_RETURN((auto [primal_solution, found_primal_feasible_solution]),
1585 GetConvexPrimalSolutionIfAvailable(model_parameters));
1586 // TODO(b/225189115): Expand QpDualsTest to check maximization problems and
1587 // other edge cases.
1588 ASSIGN_OR_RETURN((auto [dual_solution, found_dual_feasible_solution]),
1589 GetConvexDualSolutionIfAvailable(model_parameters));
1590 // TODO(b/280353996): we want to extract the basis here (when we solve via
1591 // simplex), but the existing code extracts a basis which fails our validator.
1592
1593 const SolutionClaims solution_claims = {
1594 .primal_feasible_solution_exists = found_primal_feasible_solution,
1595 .dual_feasible_solution_exists = found_dual_feasible_solution};
1596
1597 if (!primal_solution.has_value()) {
1598 return GurobiSolver::SolutionsAndClaims{.solution_claims = solution_claims};
1599 }
1600 SolutionsAndClaims solution_and_claims{.solution_claims = solution_claims};
1601 SolutionProto& solution =
1602 solution_and_claims.solutions.emplace_back(SolutionProto());
1603 if (primal_solution.has_value()) {
1604 *solution.mutable_primal_solution() = *std::move(primal_solution);
1605 }
1606 if (dual_solution.has_value()) {
1607 *solution.mutable_dual_solution() = *std::move(dual_solution);
1608 }
1609 return solution_and_claims;
1610}
1611
1612absl::StatusOr<GurobiSolver::SolutionsAndClaims> GurobiSolver::GetQcpSolution(
1613 const ModelSolveParametersProto& model_parameters) {
1614 ASSIGN_OR_RETURN((auto [primal_solution, found_primal_feasible_solution]),
1615 GetConvexPrimalSolutionIfAvailable(model_parameters));
1616 ASSIGN_OR_RETURN((auto [dual_solution, found_dual_feasible_solution]),
1617 GetConvexDualSolutionIfAvailable(model_parameters));
1618 ASSIGN_OR_RETURN(const int grb_termination,
1619 gurobi_->GetIntAttr(GRB_INT_ATTR_STATUS));
1620 // By default, Gurobi will not return duals for optimally solved QCPs.
1621 const bool proven_feasible = grb_termination == GRB_OPTIMAL;
1622 const SolutionClaims solution_claims = {
1623 .primal_feasible_solution_exists = found_primal_feasible_solution,
1624 .dual_feasible_solution_exists =
1625 found_dual_feasible_solution || proven_feasible};
1626
1627 SolutionsAndClaims solution_and_claims{.solution_claims = solution_claims};
1628 if (primal_solution.has_value() || dual_solution.has_value()) {
1629 SolutionProto& solution =
1630 solution_and_claims.solutions.emplace_back(SolutionProto());
1631 if (primal_solution.has_value()) {
1632 *solution.mutable_primal_solution() = *std::move(primal_solution);
1633 }
1634 if (dual_solution.has_value()) {
1635 *solution.mutable_dual_solution() = *std::move(dual_solution);
1636 }
1637 }
1638 return solution_and_claims;
1639}
1640
1641absl::Status GurobiSolver::SetParameters(
1642 const SolveParametersProto& parameters,
1643 const ModelSolveParametersProto& model_parameters) {
1644 ASSIGN_OR_RETURN(const bool is_mip, IsMIP());
1645 ASSIGN_OR_RETURN(const GurobiParametersProto gurobi_parameters,
1646 MergeParameters(parameters, model_parameters, is_mip,
1647 is_multi_objective_mode()));
1648 std::vector<std::string> parameter_errors;
1649 for (const GurobiParametersProto::Parameter& parameter :
1650 gurobi_parameters.parameters()) {
1651 absl::Status param_status =
1652 gurobi_->SetParam(parameter.name().c_str(), parameter.value());
1653 if (!param_status.ok()) {
1654 parameter_errors.emplace_back(std::move(param_status).message());
1655 }
1656 }
1657 if (!parameter_errors.empty()) {
1658 return absl::InvalidArgumentError(absl::StrJoin(parameter_errors, "; "));
1659 }
1660 return absl::OkStatus();
1661}
1662
1663absl::Status GurobiSolver::AddNewVariables(
1664 const VariablesProto& new_variables) {
1665 const int num_new_variables = new_variables.lower_bounds().size();
1666 std::vector<char> variable_type(num_new_variables);
1667 for (int j = 0; j < num_new_variables; ++j) {
1668 const VariableId id = new_variables.ids(j);
1669 gtl::InsertOrDie(&variables_map_, id, j + num_gurobi_variables_);
1670 variable_type[j] = new_variables.integers(j) ? GRB_INTEGER : GRB_CONTINUOUS;
1671 }
1672 // We need to copy the names, RepeatedPtrField cannot be converted to
1673 // absl::Span<std::string>.
1674 const std::vector<std::string> variable_names =
1675 TruncateNames(new_variables.names());
1676 RETURN_IF_ERROR(gurobi_->AddVars(
1677 /*obj=*/{},
1678 /*lb=*/new_variables.lower_bounds(),
1679 /*ub=*/new_variables.upper_bounds(),
1680 /*vtype=*/variable_type, variable_names));
1681 num_gurobi_variables_ += num_new_variables;
1682
1683 return absl::OkStatus();
1684}
1685
1686absl::Status GurobiSolver::AddSingleObjective(const ObjectiveProto& objective) {
1687 const int model_sense = objective.maximize() ? GRB_MAXIMIZE : GRB_MINIMIZE;
1688 RETURN_IF_ERROR(gurobi_->SetIntAttr(GRB_INT_ATTR_MODELSENSE, model_sense));
1690 gurobi_->SetDoubleAttr(GRB_DBL_ATTR_OBJCON, objective.offset()));
1691 RETURN_IF_ERROR(UpdateDoubleListAttribute(objective.linear_coefficients(),
1692 GRB_DBL_ATTR_OBJ, variables_map_));
1694 ResetQuadraticObjectiveTerms(objective.quadratic_coefficients()));
1695 return absl::OkStatus();
1696}
1697
1698absl::Status GurobiSolver::AddMultiObjectives(
1699 const ObjectiveProto& primary_objective,
1700 const google::protobuf::Map<int64_t, ObjectiveProto>&
1701 auxiliary_objectives) {
1702 absl::flat_hash_set<int64_t> priorities = {primary_objective.priority()};
1703 for (const auto& [id, objective] : auxiliary_objectives) {
1704 const int64_t priority = objective.priority();
1705 if (!priorities.insert(priority).second) {
1707 << "repeated objective priority: " << priority;
1708 }
1709 }
1710 const bool is_maximize = primary_objective.maximize();
1711 RETURN_IF_ERROR(gurobi_->SetIntAttr(
1713 RETURN_IF_ERROR(AddNewMultiObjective(
1714 primary_objective, /*objective_id=*/std::nullopt, is_maximize));
1715 for (const int64_t id : SortedMapKeys(auxiliary_objectives)) {
1717 AddNewMultiObjective(auxiliary_objectives.at(id), id, is_maximize));
1718 }
1719 return absl::OkStatus();
1720}
1721
1722absl::Status GurobiSolver::AddNewMultiObjective(
1723 const ObjectiveProto& objective,
1724 const std::optional<AuxiliaryObjectiveId> objective_id,
1725 const bool is_maximize) {
1726 std::vector<GurobiVariableIndex> var_indices;
1727 var_indices.reserve(objective.linear_coefficients().ids_size());
1728 for (const int64_t var_id : objective.linear_coefficients().ids()) {
1729 var_indices.push_back(variables_map_.at(var_id));
1730 }
1731 const GurobiMultiObjectiveIndex grb_index =
1732 static_cast<int>(multi_objectives_map_.size());
1733 // * MathOpt and Gurobi have different priority orderings (lower and higher
1734 // are more important, respectively). Therefore, we negate priorities from
1735 // MathOpt (which is OK as they are restricted to be nonnegative in MathOpt,
1736 // but not in Gurobi).
1737 // * Tolerances are set to default values, as of Gurobi v9.5.
1738 // * Gurobi exposes only a single objective sense for the entire model. We use
1739 // the objective weight to handle mixing senses across objectives (weight of
1740 // 1 if objective sense agrees with model sense, -1 otherwise).
1741 RETURN_IF_ERROR(gurobi_->SetNthObjective(
1742 /*index=*/grb_index, /*priority=*/static_cast<int>(-objective.priority()),
1743 /*weight=*/objective.maximize() == is_maximize ? +1.0 : -1.0,
1744 /*abs_tol=*/1.0e-6,
1745 /*rel_tol=*/0.0, /*name=*/objective.name(),
1746 /*constant=*/objective.offset(), /*lind=*/var_indices,
1747 /*lval=*/objective.linear_coefficients().values()));
1748 multi_objectives_map_.insert({objective_id, grb_index});
1749 return absl::OkStatus();
1750}
1751
1752// Given a vector of pairs<LinearConstraintId, LinearConstraintData&> add a
1753// slack variable for each of the constraints in the underlying `gurobi_` using
1754// the referenced bounds.
1755absl::Status GurobiSolver::AddNewSlacks(
1756 const std::vector<LinearConstraintData*>& new_slacks) {
1757 // Note that we are really adding the sub-matrix
1758 // D * slack
1759 // to the set of linear constraints, and the D matrix is stored in compressed
1760 // sparse column (CSC) format. In our particular case, D is a diagonal matrix
1761 // with -1.0 coefficients for each new slack in the row indicated in the
1762 // row_indices vector.
1763 const int num_slacks = new_slacks.size();
1764 if (num_slacks == 0) {
1765 return absl::OkStatus();
1766 }
1767 // Build the D matrix in CSC format.
1768 const std::vector<double> column_non_zeros(num_slacks, -1.0);
1769 std::vector<double> lower_bounds;
1770 std::vector<double> upper_bounds;
1771 const std::vector<char> vtypes(num_slacks, GRB_CONTINUOUS);
1772 std::vector<GurobiLinearConstraintIndex> row_indices;
1773 std::vector<int> column_non_zero_begin;
1774 column_non_zero_begin.reserve(num_slacks);
1775 row_indices.reserve(num_slacks);
1776 lower_bounds.reserve(num_slacks);
1777 upper_bounds.reserve(num_slacks);
1778 for (int k = 0; k < num_slacks; ++k) {
1779 CHECK_NE(new_slacks[k], nullptr);
1780 const LinearConstraintData& constraint_data = *new_slacks[k];
1781 row_indices.push_back(constraint_data.constraint_index);
1782 lower_bounds.push_back(constraint_data.lower_bound);
1783 upper_bounds.push_back(constraint_data.upper_bound);
1784 column_non_zero_begin.push_back(k);
1785 }
1786 // Add variables to the underlying model.
1787 RETURN_IF_ERROR(gurobi_->AddVars(/*vbegin=*/column_non_zero_begin,
1788 /*vind=*/row_indices,
1789 /*vval=*/column_non_zeros, /*obj=*/{},
1790 /*lb=*/lower_bounds, /*ub=*/upper_bounds,
1791 /*vtype=*/vtypes, /*names=*/{}));
1792 num_gurobi_variables_ += num_slacks;
1793 return absl::OkStatus();
1794}
1795
1796absl::Status GurobiSolver::AddNewLinearConstraints(
1797 const LinearConstraintsProto& constraints) {
1798 const int num_new_constraints = constraints.lower_bounds().size();
1799
1800 // We need to copy the names, RepeatedPtrField cannot be converted to
1801 // absl::Span<std::string>.
1802 const std::vector<std::string> constraint_names =
1803 TruncateNames(constraints.names());
1804 // Constraints are translated into:
1805 // 1. ax <= upper_bound (if lower bound <= -GRB_INFINITY, and upper_bound
1806 // is finite and less than GRB_INFINITY)
1807 // 2. ax >= lower_bound (if upper bound >= GRB_INFINITY, and lower_bound is
1808 // finite and greater than -GRB_INFINITY)
1809 // 3. ax == xxxxx_bound (if both bounds are finite, equal, and their
1810 // absolute values less than GRB_INFINITY)
1811 // 4. ax - slack = 0.0 (otherwise,
1812 // slack bounds == [lower_bound, upper_bound])
1813 std::vector<double> constraint_rhs;
1814 std::vector<char> constraint_sense;
1815 std::vector<LinearConstraintData*> new_slacks;
1816 constraint_rhs.reserve(num_new_constraints);
1817 constraint_sense.reserve(num_new_constraints);
1818 new_slacks.reserve(num_new_constraints);
1819 for (int i = 0; i < num_new_constraints; ++i) {
1820 const int64_t id = constraints.ids(i);
1821 LinearConstraintData& constraint_data =
1822 gtl::InsertKeyOrDie(&linear_constraints_map_, id);
1823 const double lb = constraints.lower_bounds(i);
1824 const double ub = constraints.upper_bounds(i);
1825 RETURN_IF_ERROR(SafeGurobiDouble(lb))
1826 << "lower bound for linear constraint " << id << ": "
1827 << EscapedNameForLogging(
1828 constraints.names().empty() ? "" : constraints.names(i));
1829 RETURN_IF_ERROR(SafeGurobiDouble(ub))
1830 << "upper bound for linear constraint " << id << ": "
1831 << EscapedNameForLogging(
1832 constraints.names().empty() ? "" : constraints.names(i));
1833 constraint_data.lower_bound = lb;
1834 constraint_data.upper_bound = ub;
1835 constraint_data.constraint_index = i + num_gurobi_lin_cons_;
1836 char sense = GRB_EQUAL;
1837 double rhs = 0.0;
1838 const bool lb_is_grb_neg_inf = lb <= -GRB_INFINITY;
1839 const bool ub_is_grb_pos_inf = ub >= GRB_INFINITY;
1840 if (lb_is_grb_neg_inf && !ub_is_grb_pos_inf) {
1841 sense = GRB_LESS_EQUAL;
1842 rhs = ub;
1843 } else if (!lb_is_grb_neg_inf && ub_is_grb_pos_inf) {
1844 sense = GRB_GREATER_EQUAL;
1845 rhs = lb;
1846 } else if (lb == ub) {
1847 sense = GRB_EQUAL;
1848 rhs = lb;
1849 } else {
1850 // Note that constraints where the lower bound and the upper bound are
1851 // -+infinity translate into a range constraint with an unbounded slack.
1852 constraint_data.slack_index = new_slacks.size() + num_gurobi_variables_;
1853 new_slacks.push_back(&constraint_data);
1854 }
1855 constraint_rhs.emplace_back(rhs);
1856 constraint_sense.emplace_back(sense);
1857 }
1858 // Add all constraints in one call.
1860 gurobi_->AddConstrs(constraint_sense, constraint_rhs, constraint_names));
1861 num_gurobi_lin_cons_ += num_new_constraints;
1862 // Add slacks for true ranged constraints (if needed)
1863 if (!new_slacks.empty()) {
1864 RETURN_IF_ERROR(AddNewSlacks(new_slacks));
1865 }
1866 return absl::OkStatus();
1867}
1868
1869absl::Status GurobiSolver::AddNewQuadraticConstraints(
1870 const google::protobuf::Map<QuadraticConstraintId,
1871 QuadraticConstraintProto>& constraints) {
1872 // Constraints are translated into:
1873 // 1. ax <= upper_bound (if lower bound <= -GRB_INFINITY, and upper_bound
1874 // is finite and less than GRB_INFINITY)
1875 // 2. ax >= lower_bound (if upper bound >= GRB_INFINITY, and lower_bound is
1876 // finite and greater than -GRB_INFINITY)
1877 // 3. ax == xxxxx_bound (if both bounds are finite, equal, and their
1878 // absolute values less than GRB_INFINITY)
1879 // 4. Return an error otherwise, we do not currently support ranged quadratic
1880 // constraints.
1881 for (const auto& [id, constraint] : constraints) {
1882 char sense = GRB_EQUAL;
1883 double rhs = 0.0;
1884 const double lb = constraint.lower_bound();
1885 const double ub = constraint.upper_bound();
1886 RETURN_IF_ERROR(SafeGurobiDouble(lb))
1887 << "lower bound for quadratic constraint " << id << ": "
1888 << EscapedNameForLogging(constraint.name());
1889 RETURN_IF_ERROR(SafeGurobiDouble(ub))
1890 << "upper bound for quadratic constraint " << id << ": "
1891 << EscapedNameForLogging(constraint.name());
1892 const bool lb_is_grb_neg_inf = lb <= -GRB_INFINITY;
1893 const bool ub_is_grb_pos_inf = ub >= GRB_INFINITY;
1894 if (lb_is_grb_neg_inf && ub_is_grb_pos_inf) {
1895 // The constraint is vacuous, so we just skip it.
1896 // TODO(b/227217735): Ensure duals properly account for this constraint.
1897 continue;
1898 } else if (lb_is_grb_neg_inf && !ub_is_grb_pos_inf) {
1899 sense = GRB_LESS_EQUAL;
1900 rhs = ub;
1901 } else if (!lb_is_grb_neg_inf && ub_is_grb_pos_inf) {
1902 sense = GRB_GREATER_EQUAL;
1903 rhs = lb;
1904 } else if (lb == ub) {
1905 sense = GRB_EQUAL;
1906 rhs = lb;
1907 } else {
1908 // We do not currently support ranged quadratic constraints, though it is
1909 // possible to support this if there is a need.
1910 return absl::UnimplementedError(
1911 "ranged quadratic constraints are not currently supported in Gurobi "
1912 "interface");
1913 }
1914 const SparseDoubleVectorProto& linear_coeffs = constraint.linear_terms();
1915 const int num_linear_coeffs = linear_coeffs.ids_size();
1916 std::vector<GurobiVariableIndex> linear_col_index(num_linear_coeffs);
1917 for (int k = 0; k < num_linear_coeffs; ++k) {
1918 linear_col_index[k] = variables_map_.at(linear_coeffs.ids(k));
1919 }
1920 const SparseDoubleMatrixProto& quad_coeffs = constraint.quadratic_terms();
1921 const int num_quad_coeffs = quad_coeffs.row_ids_size();
1922 std::vector<GurobiVariableIndex> quad_row_index(num_quad_coeffs);
1923 std::vector<GurobiVariableIndex> quad_col_index(num_quad_coeffs);
1924 for (int k = 0; k < num_quad_coeffs; ++k) {
1925 quad_row_index[k] = variables_map_.at(quad_coeffs.row_ids(k));
1926 quad_col_index[k] = variables_map_.at(quad_coeffs.column_ids(k));
1927 }
1928 RETURN_IF_ERROR(gurobi_->AddQConstr(
1929 linear_col_index, linear_coeffs.values(), quad_row_index,
1930 quad_col_index, quad_coeffs.coefficients(), sense, rhs,
1931 TruncateName(constraint.name())));
1932 gtl::InsertOrDie(&quadratic_constraints_map_, id, num_gurobi_quad_cons_);
1933 ++num_gurobi_quad_cons_;
1934 }
1935 return absl::OkStatus();
1936}
1937
1938std::optional<GurobiSolver::VariableId> GurobiSolver::TryExtractVariable(
1939 const LinearExpressionProto& expression) {
1940 if (expression.offset() == 0 && expression.ids_size() == 1 &&
1941 expression.coefficients(0) == 1) {
1942 return expression.ids(0);
1943 }
1944 return std::nullopt;
1945}
1946
1947absl::StatusOr<GurobiSolver::VariableEqualToExpression>
1948GurobiSolver::CreateSlackVariableEqualToExpression(
1949 const LinearExpressionProto& expression) {
1950 const GurobiVariableIndex slack_variable = num_gurobi_variables_;
1951 std::vector<GurobiVariableIndex> slack_col_indices = {slack_variable};
1952 std::vector<double> slack_coeffs = {-1.0};
1953 for (int j = 0; j < expression.ids_size(); ++j) {
1954 slack_col_indices.push_back(variables_map_.at(expression.ids(j)));
1955 slack_coeffs.push_back(expression.coefficients(j));
1956 }
1957 RETURN_IF_ERROR(gurobi_->AddVar(0, -kInf, kInf, GRB_CONTINUOUS, ""));
1958 RETURN_IF_ERROR(gurobi_->AddConstr(slack_col_indices, slack_coeffs, GRB_EQUAL,
1959 -expression.offset(), ""));
1960 return VariableEqualToExpression{.variable_index = num_gurobi_variables_++,
1961 .constraint_index = num_gurobi_lin_cons_++};
1962}
1963
1964absl::Status GurobiSolver::AddNewSecondOrderConeConstraints(
1965 const google::protobuf::Map<SecondOrderConeConstraintId,
1966 SecondOrderConeConstraintProto>& constraints) {
1967 for (const auto& [id, constraint] : constraints) {
1968 // The MathOpt proto representation for a second-order cone constraint is:
1969 // ||`constraint`.arguments_to_norm||_2 <= `constraint`.upper_bound.
1970 // Gurobi requires second-order cone constraints to be passed via the
1971 // quadratic constraint API as:
1972 // arg_var[0]^2 + ... + arg_var[d]^2 <= ub_var^2
1973 // ub_var >= 0,
1974 // for variables arg_var[0], ..., arg_var[d], ub_var. To get to this form,
1975 // we add slack variables:
1976 // ub_var = `constraint`.upper_bound
1977 // arg_var[i] = `constraint`.arguments_to_norm[i] for each i
1978 // Note that we elide adding a slack variable/constraint if the expression
1979 // we are setting it equal to is just a variable already in the model.
1980 SecondOrderConeConstraintData& constraint_data =
1981 gtl::InsertKeyOrDie(&soc_constraints_map_, id);
1982 constraint_data.constraint_index = num_gurobi_quad_cons_;
1983 // We force a new variable to be added so that we can add a lower bound on
1984 // it. Otherwise, we must update the model to flush bounds, or risk either
1985 // a Gurobi error, or stomping on a potentially stronger bound.
1987 (const auto [ub_var, ub_cons]),
1988 CreateSlackVariableEqualToExpression(constraint.upper_bound()));
1990 gurobi_->SetDoubleAttrElement(GRB_DBL_ATTR_LB, ub_var, 0.0));
1991 constraint_data.slack_variables.push_back(ub_var);
1992 constraint_data.slack_constraints.push_back(ub_cons);
1993 std::vector<GurobiVariableIndex> quad_var_indices = {ub_var};
1994 std::vector<double> quad_coeffs = {-1.0};
1995 for (const LinearExpressionProto& expression :
1996 constraint.arguments_to_norm()) {
1997 quad_coeffs.push_back(1.0);
1998 if (const std::optional<VariableId> maybe_variable =
1999 TryExtractVariable(expression);
2000 maybe_variable.has_value()) {
2001 quad_var_indices.push_back(variables_map_.at(*maybe_variable));
2002 continue;
2003 }
2004 ASSIGN_OR_RETURN((const auto [arg_var, arg_cons]),
2005 CreateSlackVariableEqualToExpression(expression));
2006 quad_var_indices.push_back(arg_var);
2007 constraint_data.slack_variables.push_back(arg_var);
2008 constraint_data.slack_constraints.push_back(arg_cons);
2009 }
2010 RETURN_IF_ERROR(gurobi_->AddQConstr({}, {}, quad_var_indices,
2011 quad_var_indices, quad_coeffs,
2012 GRB_LESS_EQUAL, 0.0, ""));
2013 ++num_gurobi_quad_cons_;
2014 }
2015 return absl::OkStatus();
2016}
2017
2018absl::Status GurobiSolver::AddNewSosConstraints(
2019 const google::protobuf::Map<AnyConstraintId, SosConstraintProto>&
2020 constraints,
2021 const int sos_type,
2022 absl::flat_hash_map<int64_t, SosConstraintData>& constraints_map) {
2023 for (const auto& [id, constraint] : constraints) {
2024 SosConstraintData& constraint_data =
2025 gtl::InsertKeyOrDie(&constraints_map, id);
2026 constraint_data.constraint_index = num_gurobi_sos_cons_;
2027 std::vector<GurobiVariableIndex> sos_var_indices;
2028 std::vector<double> weights;
2029 // As of v9.0.2, Gurobi does not allow SOS constraints to contain repeated
2030 // variables. So, we track the variables we "reuse" (i.e., were already
2031 // present in the model). Slack variables introduced in
2032 // `ExtractVariableEqualToExpression()` will not be present in the proto
2033 // inputs, so we can safely ignore tracking them.
2034 absl::flat_hash_set<VariableId> reused_variables;
2035 for (int i = 0; i < constraint.expressions_size(); ++i) {
2036 weights.push_back(constraint.weights().empty() ? i + 1
2037 : constraint.weights(i));
2038 if (const std::optional<VariableId> maybe_variable =
2039 TryExtractVariable(constraint.expressions(i));
2040 maybe_variable.has_value() &&
2041 !reused_variables.contains(*maybe_variable)) {
2042 sos_var_indices.push_back(variables_map_.at(*maybe_variable));
2043 reused_variables.insert(*maybe_variable);
2044 if (sos_type == 2) {
2045 // If this variable is deleted, Gurobi will drop the corresponding
2046 // term from the SOS constraint, potentially changing the meaning of
2047 // the constraint.
2048 undeletable_variables_.insert(*maybe_variable);
2049 }
2050 continue;
2051 }
2053 (const auto [var_index, cons_index]),
2054 CreateSlackVariableEqualToExpression(constraint.expressions(i)));
2055 sos_var_indices.push_back(var_index);
2056 constraint_data.slack_variables.push_back(var_index);
2057 constraint_data.slack_constraints.push_back(cons_index);
2058 }
2059 RETURN_IF_ERROR(gurobi_->AddSos({sos_type}, {0}, sos_var_indices, weights));
2060 ++num_gurobi_sos_cons_;
2061 }
2062 return absl::OkStatus();
2063}
2064
2065absl::Status GurobiSolver::AddNewIndicatorConstraints(
2066 const google::protobuf::Map<IndicatorConstraintId,
2067 IndicatorConstraintProto>& constraints) {
2068 for (const auto& [id, constraint] : constraints) {
2069 if (!constraint.has_indicator_id()) {
2070 if (constraint.activate_on_zero()) {
2072 << "MathOpt does not currently support Gurobi models with "
2073 "indicator constraints that activate on zero with unset "
2074 "indicator variables; encountered one at id: "
2075 << id;
2076 }
2077 gtl::InsertOrDie(&indicator_constraints_map_, id, std::nullopt);
2078 continue;
2079 }
2080 const int num_terms = constraint.expression().ids_size();
2081 std::vector<GurobiVariableIndex> grb_ids(num_terms);
2082 for (int k = 0; k < num_terms; ++k) {
2083 grb_ids[k] = variables_map_.at(constraint.expression().ids(k));
2084 }
2085 char sense = GRB_EQUAL;
2086 double rhs = 0.0;
2087 const double lb = constraint.lower_bound();
2088 const double ub = constraint.upper_bound();
2089 RETURN_IF_ERROR(SafeGurobiDouble(lb))
2090 << "lower bound for indicator constraint " << id << ": "
2091 << EscapedNameForLogging(constraint.name());
2092 RETURN_IF_ERROR(SafeGurobiDouble(ub))
2093 << "upper bound for indicator constraint " << id << ": "
2094 << EscapedNameForLogging(constraint.name());
2095 const bool lb_is_grb_neg_inf = lb <= -GRB_INFINITY;
2096 const bool ub_is_grb_pos_inf = ub >= GRB_INFINITY;
2097 if (lb_is_grb_neg_inf && ub_is_grb_pos_inf) {
2098 // The constraint is vacuous, so we just skip it.
2099 continue;
2100 } else if (lb_is_grb_neg_inf && !ub_is_grb_pos_inf) {
2101 sense = GRB_LESS_EQUAL;
2102 rhs = ub;
2103 } else if (!lb_is_grb_neg_inf && ub_is_grb_pos_inf) {
2104 sense = GRB_GREATER_EQUAL;
2105 rhs = lb;
2106 } else if (lb == ub) {
2107 sense = GRB_EQUAL;
2108 rhs = lb;
2109 } else {
2110 // We do not currently support ranged indicator constraints, though it
2111 // is possible to support this if there is a need.
2112 return absl::UnimplementedError(
2113 "ranged indicator constraints are not currently supported in "
2114 "Gurobi "
2115 "interface");
2116 }
2117 RETURN_IF_ERROR(gurobi_->AddIndicator(
2118 /*name=*/constraint.name(),
2119 /*binvar=*/variables_map_.at(constraint.indicator_id()),
2120 /*binval=*/constraint.activate_on_zero() ? 0 : 1,
2121 /*ind=*/grb_ids, /*val=*/constraint.expression().values(),
2122 /*sense=*/sense, /*rhs=*/rhs));
2123 gtl::InsertOrDie(&indicator_constraints_map_, id,
2124 IndicatorConstraintData{
2125 .constraint_index = num_gurobi_gen_cons_,
2126 .indicator_variable_id = constraint.indicator_id()});
2127 ++num_gurobi_gen_cons_;
2128 // Deleting the indicator variable, but not the associated indicator
2129 // constraint, will lead to a Gurobi error.
2130 undeletable_variables_.insert(constraint.indicator_id());
2131 }
2132 return absl::OkStatus();
2133}
2134
2135absl::Status GurobiSolver::ChangeCoefficients(
2136 const SparseDoubleMatrixProto& matrix) {
2137 const int num_coefficients = matrix.row_ids().size();
2138 std::vector<GurobiLinearConstraintIndex> row_index(num_coefficients);
2139 std::vector<GurobiVariableIndex> col_index(num_coefficients);
2140 for (int k = 0; k < num_coefficients; ++k) {
2141 row_index[k] =
2142 linear_constraints_map_.at(matrix.row_ids(k)).constraint_index;
2143 col_index[k] = variables_map_.at(matrix.column_ids(k));
2144 }
2145 return gurobi_->ChgCoeffs(row_index, col_index, matrix.coefficients());
2146}
2147
2148absl::Status GurobiSolver::UpdateDoubleListAttribute(
2149 const SparseDoubleVectorProto& update,
2150 const char* absl_nonnull attribute_name, const IdHashMap& id_hash_map) {
2151 if (update.ids_size() == 0) {
2152 return absl::OkStatus();
2153 }
2154 std::vector<int> index;
2155 index.reserve(update.ids_size());
2156 for (const int64_t id : update.ids()) {
2157 index.push_back(id_hash_map.at(id));
2158 }
2159 return gurobi_->SetDoubleAttrList(attribute_name, index, update.values());
2160}
2161
2162absl::Status GurobiSolver::UpdateInt32ListAttribute(
2163 const SparseInt32VectorProto& update,
2164 const char* absl_nonnull attribute_name, const IdHashMap& id_hash_map) {
2165 if (update.ids_size() == 0) {
2166 return absl::OkStatus();
2167 }
2168 std::vector<int> index;
2169 index.reserve(update.ids_size());
2170 for (const int64_t id : update.ids()) {
2171 index.push_back(id_hash_map.at(id));
2172 }
2173 return gurobi_->SetIntAttrList(attribute_name, index, update.values());
2174}
2175
2176absl::Status GurobiSolver::LoadModel(const ModelProto& input_model) {
2177 CHECK(gurobi_ != nullptr);
2178 RETURN_IF_ERROR(gurobi_->SetStringAttr(GRB_STR_ATTR_MODELNAME,
2179 TruncateName(input_model.name())));
2180 RETURN_IF_ERROR(AddNewVariables(input_model.variables()));
2181
2182 RETURN_IF_ERROR(AddNewLinearConstraints(input_model.linear_constraints()));
2184 AddNewQuadraticConstraints(input_model.quadratic_constraints()));
2185 RETURN_IF_ERROR(AddNewSecondOrderConeConstraints(
2186 input_model.second_order_cone_constraints()));
2187 RETURN_IF_ERROR(AddNewSosConstraints(input_model.sos1_constraints(),
2188 GRB_SOS_TYPE1, sos1_constraints_map_));
2189 RETURN_IF_ERROR(AddNewSosConstraints(input_model.sos2_constraints(),
2190 GRB_SOS_TYPE2, sos2_constraints_map_));
2192 AddNewIndicatorConstraints(input_model.indicator_constraints()));
2193
2194 RETURN_IF_ERROR(ChangeCoefficients(input_model.linear_constraint_matrix()));
2195
2196 if (input_model.auxiliary_objectives().empty()) {
2197 RETURN_IF_ERROR(AddSingleObjective(input_model.objective()));
2198 } else {
2199 RETURN_IF_ERROR(AddMultiObjectives(input_model.objective(),
2200 input_model.auxiliary_objectives()));
2201 }
2202 return absl::OkStatus();
2203}
2204
2205absl::Status GurobiSolver::ResetQuadraticObjectiveTerms(
2206 const SparseDoubleMatrixProto& terms) {
2207 quadratic_objective_coefficients_.clear();
2208 RETURN_IF_ERROR(gurobi_->DelQ());
2209 const int num_terms = terms.row_ids().size();
2210 if (num_terms > 0) {
2211 std::vector<GurobiVariableIndex> first_var_index(num_terms);
2212 std::vector<GurobiVariableIndex> second_var_index(num_terms);
2213 for (int k = 0; k < num_terms; ++k) {
2214 const VariableId row_id = terms.row_ids(k);
2215 const VariableId column_id = terms.column_ids(k);
2216 first_var_index[k] = variables_map_.at(row_id);
2217 second_var_index[k] = variables_map_.at(column_id);
2218 quadratic_objective_coefficients_[{row_id, column_id}] =
2219 terms.coefficients(k);
2220 }
2221 RETURN_IF_ERROR(gurobi_->AddQpTerms(first_var_index, second_var_index,
2222 terms.coefficients()));
2223 }
2224 return absl::OkStatus();
2225}
2226
2227absl::Status GurobiSolver::UpdateQuadraticObjectiveTerms(
2228 const SparseDoubleMatrixProto& terms) {
2229 CHECK(gurobi_ != nullptr);
2230 const int num_terms = terms.row_ids().size();
2231 if (num_terms > 0) {
2232 std::vector<GurobiVariableIndex> first_var_index(num_terms);
2233 std::vector<GurobiVariableIndex> second_var_index(num_terms);
2234 std::vector<double> coefficient_updates(num_terms);
2235 for (int k = 0; k < num_terms; ++k) {
2236 const VariableId row_id = terms.row_ids(k);
2237 const VariableId column_id = terms.column_ids(k);
2238 first_var_index[k] = variables_map_.at(row_id);
2239 second_var_index[k] = variables_map_.at(column_id);
2240 const std::pair<VariableId, VariableId> qp_term_key(row_id, column_id);
2241 const double new_coefficient = terms.coefficients(k);
2242 // Gurobi will maintain any existing quadratic coefficients unless we
2243 // call GRBdelq (which we don't). So, since stored entries in terms
2244 // specify the target coefficients, we need to compute the difference
2245 // from the existing coefficient with Gurobi, if any.
2246 coefficient_updates[k] =
2247 new_coefficient - quadratic_objective_coefficients_[qp_term_key];
2248 quadratic_objective_coefficients_[qp_term_key] = new_coefficient;
2249 }
2250 RETURN_IF_ERROR(gurobi_->AddQpTerms(first_var_index, second_var_index,
2251 coefficient_updates));
2252 }
2253 return absl::OkStatus();
2254}
2255
2256// Bound changes in constraints can induce new variables, and also remove
2257// some slacks. We first add all new variables, and queue all deletions to be
2258// dealt with later on.
2259absl::Status GurobiSolver::UpdateLinearConstraints(
2260 const LinearConstraintUpdatesProto& constraints_update,
2261 std::vector<GurobiVariableIndex>& deleted_variables_index) {
2262 const SparseDoubleVectorProto& constraint_lower_bounds =
2263 constraints_update.lower_bounds();
2264 const SparseDoubleVectorProto& constraint_upper_bounds =
2265 constraints_update.upper_bounds();
2266
2267 // If no update, just return.
2268 if (constraint_lower_bounds.ids().empty() &&
2269 constraint_upper_bounds.ids().empty()) {
2270 return absl::OkStatus();
2271 }
2272
2273 // We want to avoid changing the right-hand-side, sense, or slacks of each
2274 // constraint more than once. Since we can refer to the same constraint ID
2275 // both in the `constraint_upper_bounds` and `constraint_lower_bounds`
2276 // sparse vectors, we collect all changes into a single structure:
2277 struct UpdateConstraintData {
2278 LinearConstraintId constraint_id;
2279 LinearConstraintData& source;
2280 double new_lower_bound;
2281 double new_upper_bound;
2282 UpdateConstraintData(const LinearConstraintId id,
2283 LinearConstraintData& reference)
2284 : constraint_id(id),
2285 source(reference),
2286 new_lower_bound(reference.lower_bound),
2287 new_upper_bound(reference.upper_bound) {}
2288 };
2289 const int upper_bounds_size = constraint_upper_bounds.ids().size();
2290 const int lower_bounds_size = constraint_lower_bounds.ids().size();
2291 std::vector<UpdateConstraintData> update_vector;
2292 update_vector.reserve(upper_bounds_size + lower_bounds_size);
2293 // We exploit the fact that IDs are sorted in increasing order to merge
2294 // changes into a vector of aggregated changes.
2295 for (int lower_index = 0, upper_index = 0;
2296 lower_index < lower_bounds_size || upper_index < upper_bounds_size;) {
2297 VariableId lower_id = std::numeric_limits<int64_t>::max();
2298 if (lower_index < lower_bounds_size) {
2299 lower_id = constraint_lower_bounds.ids(lower_index);
2300 }
2301 VariableId upper_id = std::numeric_limits<int64_t>::max();
2302 if (upper_index < upper_bounds_size) {
2303 upper_id = constraint_upper_bounds.ids(upper_index);
2304 }
2305 const VariableId id = std::min(lower_id, upper_id);
2306 DCHECK(id < std::numeric_limits<int64_t>::max());
2307 UpdateConstraintData update(id, linear_constraints_map_.at(id));
2308 if (lower_id == upper_id) {
2309 update.new_lower_bound = constraint_lower_bounds.values(lower_index++);
2310 update.new_upper_bound = constraint_upper_bounds.values(upper_index++);
2311 } else if (lower_id < upper_id) {
2312 update.new_lower_bound = constraint_lower_bounds.values(lower_index++);
2313 } else { /* upper_id < lower_id */
2314 update.new_upper_bound = constraint_upper_bounds.values(upper_index++);
2315 }
2316 update_vector.emplace_back(update);
2317 }
2318
2319 // We have grouped all changes in update_vector, now generate changes in
2320 // slack bounds, rhs, senses, new slacks, and deleted_slacks (to be dealt
2321 // with later, outside this function).
2322 // These three vectors keep changes to right-hand-side and senses.
2323 std::vector<char> sense_data;
2324 std::vector<double> rhs_data;
2325 std::vector<GurobiLinearConstraintIndex> rhs_index;
2326 // These three vectors keep changes to bounds on existing slack.
2327 std::vector<double> lower_bound_data;
2328 std::vector<double> upper_bound_data;
2329 std::vector<GurobiVariableIndex> bound_index;
2330 // This vector keep newly introduced slacks.
2331 std::vector<LinearConstraintData*> new_slacks;
2332 // Iterate on the changes, and populate the three possible changes.
2333 for (UpdateConstraintData& update_data : update_vector) {
2334 const bool same_lower_bound =
2335 (update_data.source.lower_bound == update_data.new_lower_bound) ||
2336 ((update_data.source.lower_bound <= -GRB_INFINITY) &&
2337 (update_data.new_lower_bound <= -GRB_INFINITY));
2338 const bool same_upper_bound =
2339 (update_data.source.upper_bound == update_data.new_upper_bound) ||
2340 ((update_data.source.upper_bound >= GRB_INFINITY) &&
2341 (update_data.new_upper_bound >= GRB_INFINITY));
2342 if (same_upper_bound && same_lower_bound) continue;
2343 // Save into linear_constraints_map_[id] the new bounds for the linear
2344 // constraint.
2345 update_data.source.lower_bound = update_data.new_lower_bound;
2346 update_data.source.upper_bound = update_data.new_upper_bound;
2347 bool delete_slack = false;
2348 // Detect the type of constraint to add and store RHS and bounds.
2349 if (update_data.new_lower_bound <= -GRB_INFINITY &&
2350 update_data.new_upper_bound < GRB_INFINITY) {
2351 delete_slack = true;
2352 rhs_index.emplace_back(update_data.source.constraint_index);
2353 rhs_data.emplace_back(update_data.new_upper_bound);
2354 sense_data.emplace_back(GRB_LESS_EQUAL);
2355 } else if (update_data.new_lower_bound > -GRB_INFINITY &&
2356 update_data.new_upper_bound >= GRB_INFINITY) {
2357 delete_slack = true;
2358 rhs_index.emplace_back(update_data.source.constraint_index);
2359 rhs_data.emplace_back(update_data.new_lower_bound);
2360 sense_data.emplace_back(GRB_GREATER_EQUAL);
2361 } else if (update_data.new_lower_bound == update_data.new_upper_bound) {
2362 delete_slack = true;
2363 rhs_index.emplace_back(update_data.source.constraint_index);
2364 rhs_data.emplace_back(update_data.new_lower_bound);
2365 sense_data.emplace_back(GRB_EQUAL);
2366 } else {
2367 // Note that constraints where the lower bound and the upper bound are
2368 // -+infinity translated into a range constraint with an unbounded
2369 // slack.
2370 if (update_data.source.slack_index != kUnspecifiedIndex) {
2371 bound_index.emplace_back(update_data.source.slack_index);
2372 lower_bound_data.emplace_back(update_data.new_lower_bound);
2373 upper_bound_data.emplace_back(update_data.new_upper_bound);
2374 } else {
2375 // Note that if we add a new slack, we must both reset the sense and
2376 // right hand side for the inequality.
2377 rhs_index.emplace_back(update_data.source.constraint_index);
2378 rhs_data.emplace_back(0.0);
2379 sense_data.emplace_back(GRB_EQUAL);
2380 // Update the slack_index in the linear_constraints_map_[id]
2381 update_data.source.slack_index =
2382 new_slacks.size() + num_gurobi_variables_;
2383 // Save the data needed to add the new slack.
2384 new_slacks.push_back(&update_data.source);
2385 }
2386 }
2387 // If the constraint had a slack, and now is marked for deletion, we reset
2388 // the stored slack_index in linear_constraints_map_[id], save the index
2389 // in the list of variables to be deleted later on and remove the
2390 // constraint from slack_map_.
2391 if (delete_slack && update_data.source.slack_index != kUnspecifiedIndex) {
2392 deleted_variables_index.emplace_back(update_data.source.slack_index);
2393 update_data.source.slack_index = kUnspecifiedIndex;
2394 }
2395 }
2396
2397 // Pass down changes to Gurobi.
2398 if (!rhs_index.empty()) {
2399 RETURN_IF_ERROR(
2400 gurobi_->SetDoubleAttrList(GRB_DBL_ATTR_RHS, rhs_index, rhs_data));
2401 RETURN_IF_ERROR(
2402 gurobi_->SetCharAttrList(GRB_CHAR_ATTR_SENSE, rhs_index, sense_data));
2403 } // rhs changes
2404 if (!bound_index.empty()) {
2405 RETURN_IF_ERROR(gurobi_->SetDoubleAttrList(GRB_DBL_ATTR_LB, bound_index,
2406 lower_bound_data));
2407 RETURN_IF_ERROR(gurobi_->SetDoubleAttrList(GRB_DBL_ATTR_UB, bound_index,
2408 upper_bound_data));
2409 } // Slack bound changes.
2410
2411 if (!new_slacks.empty()) {
2412 RETURN_IF_ERROR(AddNewSlacks(new_slacks));
2413 }
2414 return absl::OkStatus();
2415}
2416
2417// This function re-assign indices for variables and constraints after
2418// deletion. The updated indices are computed from the previous indices,
2419// sorted in incremental form, but re-assigned so that all indices are
2420// contiguous between [0, num_variables-1], [0, num_linear_constraints-1], and
2421// [0, num_quad_constraints-1], etc.
2422void GurobiSolver::UpdateGurobiIndices(const DeletedIndices& deleted_indices) {
2423 // Recover the updated indices of variables.
2424 if (!deleted_indices.variables.empty()) {
2425 const std::vector<GurobiVariableIndex> old_to_new =
2426 IndexUpdateMap(num_gurobi_variables_, deleted_indices.variables);
2427 for (auto& [_, grb_index] : variables_map_) {
2428 grb_index = old_to_new[grb_index];
2429 CHECK_NE(grb_index, kDeletedIndex);
2430 }
2431 for (auto& [_, lin_con_data] : linear_constraints_map_) {
2432 if (lin_con_data.slack_index != kUnspecifiedIndex) {
2433 lin_con_data.slack_index = old_to_new[lin_con_data.slack_index];
2434 CHECK_NE(lin_con_data.slack_index, kDeletedIndex);
2435 }
2436 }
2437 for (auto& [_, soc_con_data] : soc_constraints_map_) {
2438 for (GurobiVariableIndex& index : soc_con_data.slack_variables) {
2439 index = old_to_new[index];
2440 CHECK_NE(index, kDeletedIndex);
2441 }
2442 }
2443 for (auto& [_, sos1_con_data] : sos1_constraints_map_) {
2444 for (GurobiVariableIndex& index : sos1_con_data.slack_variables) {
2445 index = old_to_new[index];
2446 CHECK_NE(index, kDeletedIndex);
2447 }
2448 }
2449 for (auto& [_, sos2_con_data] : sos2_constraints_map_) {
2450 for (GurobiVariableIndex& index : sos2_con_data.slack_variables) {
2451 index = old_to_new[index];
2452 CHECK_NE(index, kDeletedIndex);
2453 }
2454 }
2455 }
2456 // Recover the updated indices of linear constraints.
2457 if (!deleted_indices.linear_constraints.empty()) {
2458 const std::vector<GurobiLinearConstraintIndex> old_to_new = IndexUpdateMap(
2459 num_gurobi_lin_cons_, deleted_indices.linear_constraints);
2460 for (auto& [_, lin_con_data] : linear_constraints_map_) {
2461 lin_con_data.constraint_index = old_to_new[lin_con_data.constraint_index];
2462 CHECK_NE(lin_con_data.constraint_index, kDeletedIndex);
2463 }
2464 for (auto& [_, soc_con_data] : soc_constraints_map_) {
2465 for (GurobiLinearConstraintIndex& index :
2466 soc_con_data.slack_constraints) {
2467 index = old_to_new[index];
2468 CHECK_NE(index, kDeletedIndex);
2469 }
2470 }
2471 for (auto& [_, sos1_con_data] : sos1_constraints_map_) {
2472 for (GurobiLinearConstraintIndex& index :
2473 sos1_con_data.slack_constraints) {
2474 index = old_to_new[index];
2475 CHECK_NE(index, kDeletedIndex);
2476 }
2477 }
2478 for (auto& [_, sos2_con_data] : sos2_constraints_map_) {
2479 for (GurobiLinearConstraintIndex& index :
2480 sos2_con_data.slack_constraints) {
2481 index = old_to_new[index];
2482 CHECK_NE(index, kDeletedIndex);
2483 }
2484 }
2485 }
2486 // Recover the updated indices of quadratic constraints.
2487 if (!deleted_indices.quadratic_constraints.empty()) {
2488 const std::vector<GurobiQuadraticConstraintIndex> old_to_new =
2489 IndexUpdateMap(num_gurobi_quad_cons_,
2490 deleted_indices.quadratic_constraints);
2491 for (auto& [_, grb_index] : quadratic_constraints_map_) {
2492 grb_index = old_to_new[grb_index];
2493 CHECK_NE(grb_index, kDeletedIndex);
2494 }
2495 for (auto& [_, soc_con_data] : soc_constraints_map_) {
2496 GurobiQuadraticConstraintIndex& grb_index = soc_con_data.constraint_index;
2497 grb_index = old_to_new[soc_con_data.constraint_index];
2498 CHECK_NE(grb_index, kDeletedIndex);
2499 }
2500 }
2501 // Recover the updated indices of SOS constraints.
2502 if (!deleted_indices.sos_constraints.empty()) {
2503 const std::vector<GurobiSosConstraintIndex> old_to_new =
2504 IndexUpdateMap(num_gurobi_sos_cons_, deleted_indices.sos_constraints);
2505 for (auto& [_, sos1_data] : sos1_constraints_map_) {
2506 GurobiSosConstraintIndex& grb_index = sos1_data.constraint_index;
2507 grb_index = old_to_new[grb_index];
2508 CHECK_NE(grb_index, kDeletedIndex);
2509 }
2510 for (auto& [_, sos2_data] : sos2_constraints_map_) {
2511 GurobiSosConstraintIndex& grb_index = sos2_data.constraint_index;
2512 grb_index = old_to_new[grb_index];
2513 CHECK_NE(grb_index, kDeletedIndex);
2514 }
2515 }
2516 // Recover the updated indices of general constraints.
2517 if (!deleted_indices.general_constraints.empty()) {
2518 const std::vector<GurobiGeneralConstraintIndex> old_to_new = IndexUpdateMap(
2519 num_gurobi_gen_cons_, deleted_indices.general_constraints);
2520 for (auto& [_, indicator_data] : indicator_constraints_map_) {
2521 if (!indicator_data.has_value()) {
2522 continue;
2523 }
2524 GurobiGeneralConstraintIndex& grb_index =
2525 indicator_data->constraint_index;
2526 grb_index = old_to_new[grb_index];
2527 CHECK_NE(grb_index, kDeletedIndex);
2528 }
2529 }
2530}
2531
2532absl::StatusOr<bool> GurobiSolver::Update(
2533 const ModelUpdateProto& model_update) {
2534 if (!undeletable_variables_.empty()) {
2535 for (const VariableId id : model_update.deleted_variable_ids()) {
2536 if (undeletable_variables_.contains(id)) {
2537 return false;
2538 }
2539 }
2540 }
2541 if (!UpdateIsSupported(model_update, kGurobiSupportedStructures)) {
2542 return false;
2543 }
2544 // As of 2022-12-06 we do not support incrementalism for multi-objective
2545 // models: not adding/deleting/modifying the auxiliary objectives...
2546 if (const AuxiliaryObjectivesUpdatesProto& objs_update =
2547 model_update.auxiliary_objectives_updates();
2548 !objs_update.deleted_objective_ids().empty() ||
2549 !objs_update.new_objectives().empty() ||
2550 !objs_update.objective_updates().empty()) {
2551 return false;
2552 }
2553 // ...or modifying the primary objective of a multi-objective model.
2554 if (is_multi_objective_mode() && model_update.has_objective_updates()) {
2555 return false;
2556 }
2557
2558 RETURN_IF_ERROR(AddNewVariables(model_update.new_variables()));
2559
2561 AddNewLinearConstraints(model_update.new_linear_constraints()));
2562
2563 RETURN_IF_ERROR(AddNewQuadraticConstraints(
2565 RETURN_IF_ERROR(AddNewSecondOrderConeConstraints(
2567 RETURN_IF_ERROR(AddNewSosConstraints(
2569 sos1_constraints_map_));
2570 RETURN_IF_ERROR(AddNewSosConstraints(
2572 sos2_constraints_map_));
2573 RETURN_IF_ERROR(AddNewIndicatorConstraints(
2575
2577 ChangeCoefficients(model_update.linear_constraint_matrix_updates()));
2578
2579 if (model_update.objective_updates().has_direction_update()) {
2580 const int model_sense = model_update.objective_updates().direction_update()
2581 ? GRB_MAXIMIZE
2582 : GRB_MINIMIZE;
2583 RETURN_IF_ERROR(gurobi_->SetIntAttr(GRB_INT_ATTR_MODELSENSE, model_sense));
2584 }
2585
2586 if (model_update.objective_updates().has_offset_update()) {
2587 RETURN_IF_ERROR(gurobi_->SetDoubleAttr(
2589 }
2590
2591 RETURN_IF_ERROR(UpdateDoubleListAttribute(
2593 variables_map_));
2594
2595 RETURN_IF_ERROR(UpdateQuadraticObjectiveTerms(
2596 model_update.objective_updates().quadratic_coefficients()));
2597
2599 UpdateDoubleListAttribute(model_update.variable_updates().lower_bounds(),
2600 GRB_DBL_ATTR_LB, variables_map_));
2601
2603 UpdateDoubleListAttribute(model_update.variable_updates().upper_bounds(),
2604 GRB_DBL_ATTR_UB, variables_map_));
2605
2606 if (model_update.variable_updates().has_integers()) {
2607 const SparseBoolVectorProto& update =
2608 model_update.variable_updates().integers();
2609 std::vector<GurobiVariableIndex> index;
2610 index.reserve(update.ids_size());
2611 for (const int64_t id : update.ids()) {
2612 index.push_back(variables_map_.at(id));
2613 }
2614 std::vector<char> value;
2615 value.reserve(update.values_size());
2616 for (const bool val : update.values()) {
2617 value.push_back(val ? GRB_INTEGER : GRB_CONTINUOUS);
2618 }
2620 gurobi_->SetCharAttrList(GRB_CHAR_ATTR_VTYPE, index, value));
2621 }
2622
2623 // Now we update quadratic_objective_coefficients_, removing any terms where
2624 // either one or both of the involved variables are about to be deleted.
2625 const absl::flat_hash_set<VariableId> variable_ids_to_be_deleted(
2626 model_update.deleted_variable_ids().begin(),
2627 model_update.deleted_variable_ids().end());
2628 // NOTE: Introducing more state and complexity should speed this up, but we
2629 // opt for the simpler approach for now.
2630 for (auto it = quadratic_objective_coefficients_.cbegin();
2631 it != quadratic_objective_coefficients_.cend();
2632 /*incremented in loop*/) {
2633 if (variable_ids_to_be_deleted.contains(it->first.first) ||
2634 variable_ids_to_be_deleted.contains(it->first.second)) {
2635 quadratic_objective_coefficients_.erase(it++);
2636 } else {
2637 ++it;
2638 }
2639 }
2640 // We cache all Gurobi variables and constraint indices that must be
2641 // deleted, and perform deletions at the end of the update call.
2642 DeletedIndices deleted_indices;
2643
2644 RETURN_IF_ERROR(UpdateLinearConstraints(
2645 model_update.linear_constraint_updates(), deleted_indices.variables));
2646
2647 for (const VariableId id : model_update.deleted_variable_ids()) {
2648 deleted_indices.variables.emplace_back(variables_map_.at(id));
2649 variables_map_.erase(id);
2650 }
2651
2652 for (const LinearConstraintId id :
2653 model_update.deleted_linear_constraint_ids()) {
2654 LinearConstraintData& constraint_data = linear_constraints_map_.at(id);
2655 deleted_indices.linear_constraints.push_back(
2656 constraint_data.constraint_index);
2657 if (constraint_data.slack_index != kUnspecifiedIndex) {
2658 deleted_indices.variables.push_back(constraint_data.slack_index);
2659 constraint_data.slack_index = kUnspecifiedIndex;
2660 }
2661 linear_constraints_map_.erase(id);
2662 }
2663
2664 for (const QuadraticConstraintId id :
2666 const GurobiQuadraticConstraintIndex grb_index =
2667 quadratic_constraints_map_.at(id);
2668 deleted_indices.quadratic_constraints.push_back(grb_index);
2669 quadratic_constraints_map_.erase(id);
2670 }
2671
2672 for (const SecondOrderConeConstraintId id :
2675 deleted_indices.quadratic_constraints.push_back(
2676 soc_constraints_map_.at(id).constraint_index);
2677 for (const GurobiVariableIndex index :
2678 soc_constraints_map_.at(id).slack_variables) {
2679 deleted_indices.variables.push_back(index);
2680 }
2681 for (const GurobiLinearConstraintIndex index :
2682 soc_constraints_map_.at(id).slack_constraints) {
2683 deleted_indices.linear_constraints.push_back(index);
2684 }
2685 soc_constraints_map_.erase(id);
2686 }
2687
2688 const auto sos_updater = [&](const SosConstraintData& sos_constraint) {
2689 deleted_indices.sos_constraints.push_back(sos_constraint.constraint_index);
2690 for (const GurobiVariableIndex index : sos_constraint.slack_variables) {
2691 deleted_indices.variables.push_back(index);
2692 }
2693 for (const GurobiLinearConstraintIndex index :
2694 sos_constraint.slack_constraints) {
2695 deleted_indices.linear_constraints.push_back(index);
2696 }
2697 };
2698 for (const Sos1ConstraintId id :
2700 sos_updater(sos1_constraints_map_.at(id));
2701 sos1_constraints_map_.erase(id);
2702 }
2703
2704 for (const Sos2ConstraintId id :
2706 sos_updater(sos2_constraints_map_.at(id));
2707 sos2_constraints_map_.erase(id);
2708 }
2709
2710 for (const IndicatorConstraintId id :
2712 // Otherwise the constraint is not actually registered with Gurobi.
2713 const auto it = indicator_constraints_map_.find(id);
2714 CHECK(it != indicator_constraints_map_.end()) << "id: " << id;
2715 if (it->second.has_value()) {
2716 deleted_indices.general_constraints.push_back(
2717 it->second->constraint_index);
2718 }
2719 indicator_constraints_map_.erase(it);
2720 }
2721
2722 UpdateGurobiIndices(deleted_indices);
2723
2724 // If we are removing variables or constraints we remove them after adding
2725 // any variable or constraint. This is to avoid problems with
2726 // the numbering of possibly new variables and constraints.
2727 // After that we must update the model so that sequence of updates don't
2728 // interfere with one-another.
2729 if (!deleted_indices.linear_constraints.empty()) {
2730 RETURN_IF_ERROR(gurobi_->DelConstrs(deleted_indices.linear_constraints));
2731 num_gurobi_lin_cons_ -= deleted_indices.linear_constraints.size();
2732 }
2733
2734 if (!deleted_indices.quadratic_constraints.empty()) {
2736 gurobi_->DelQConstrs(deleted_indices.quadratic_constraints));
2737 num_gurobi_quad_cons_ -= deleted_indices.quadratic_constraints.size();
2738 }
2739
2740 if (!deleted_indices.sos_constraints.empty()) {
2741 RETURN_IF_ERROR(gurobi_->DelSos(deleted_indices.sos_constraints));
2742 }
2743
2744 if (!deleted_indices.general_constraints.empty()) {
2746 gurobi_->DelGenConstrs(deleted_indices.general_constraints));
2747 }
2748
2749 if (!deleted_indices.variables.empty()) {
2750 RETURN_IF_ERROR(gurobi_->DelVars(deleted_indices.variables));
2751 num_gurobi_variables_ -= deleted_indices.variables.size();
2752 }
2753
2754 // Synchronize all pending changes.
2755 RETURN_IF_ERROR(gurobi_->UpdateModel());
2756
2757 return true;
2758}
2759
2760absl::StatusOr<std::unique_ptr<GurobiSolver>> GurobiSolver::New(
2761 const ModelProto& input_model, const SolverInterface::InitArgs& init_args) {
2762 // TODO(user): Correctly load the gurobi library in open source.
2764 ModelIsSupported(input_model, kGurobiSupportedStructures, "Gurobi"));
2765 if (!input_model.auxiliary_objectives().empty() &&
2766 !input_model.objective().quadratic_coefficients().row_ids().empty()) {
2768 << "Gurobi does not support multiple objective models with "
2769 "quadratic objectives";
2770 }
2771 for (const auto& [id, obj] : input_model.auxiliary_objectives()) {
2772 if (!obj.quadratic_coefficients().row_ids().empty()) {
2774 << "Gurobi does not support multiple objective models with "
2775 "quadratic objectives";
2776 }
2777 }
2778 ASSIGN_OR_RETURN(std::unique_ptr<Gurobi> gurobi,
2779 GurobiFromInitArgs(init_args));
2780 auto gurobi_solver = absl::WrapUnique(new GurobiSolver(std::move(gurobi)));
2781 RETURN_IF_ERROR(gurobi_solver->LoadModel(input_model));
2782 return gurobi_solver;
2783}
2784
2785absl::StatusOr<std::unique_ptr<GurobiSolver::GurobiCallbackData>>
2786GurobiSolver::RegisterCallback(
2787 const CallbackRegistrationProto& registration, const Callback cb,
2788 const MessageCallback message_cb, const absl::Time start,
2789 SolveInterrupter* absl_nullable const local_interrupter) {
2790 const absl::flat_hash_set<CallbackEventProto> events = EventSet(registration);
2791
2792 // Note that IS_MIP does not necessarily mean the problem has integer
2793 // variables. Please refer to Gurobi's doc for details:
2794 // https://www.gurobi.com/documentation/9.1/refman/ismip.html.
2795 //
2796 // Here we assume that we get MIP related events and use a MIP solving
2797 // strategy when IS_MIP is true.
2798 ASSIGN_OR_RETURN(const int is_mip, gurobi_->GetIntAttr(GRB_INT_ATTR_IS_MIP));
2799
2801 registration, is_mip ? SupportedMIPEvents() : SupportedLPEvents()))
2802 << "for a " << (is_mip ? "MIP" : "LP") << " model";
2803
2804 // Set Gurobi parameters.
2805 if (message_cb != nullptr) {
2806 // Disable logging messages to the console the user wants to handle
2807 // messages.
2808 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_LOGTOCONSOLE, 0));
2809 }
2810 if (registration.add_cuts() || registration.add_lazy_constraints()) {
2811 // This is to signal the solver presolve to limit primal transformations
2812 // that precludes crushing cuts to the presolved model.
2813 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_PRECRUSH, 1));
2814 }
2815 if (registration.add_lazy_constraints()) {
2816 // This is needed so that the solver knows that some presolve reductions
2817 // can not be performed safely.
2818 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_LAZYCONSTRAINTS, 1));
2819 }
2820 return std::make_unique<GurobiCallbackData>(
2821 GurobiCallbackInput{
2822 .user_cb = cb,
2823 .message_cb = message_cb,
2824 .variable_ids = variables_map_,
2825 .num_gurobi_vars = num_gurobi_variables_,
2826 .events = EventToGurobiWhere(events),
2827 .mip_solution_filter = registration.mip_solution_filter(),
2828 .mip_node_filter = registration.mip_node_filter(),
2829 .start = start},
2830 local_interrupter);
2831}
2832
2833absl::StatusOr<InvertedBounds> GurobiSolver::ListInvertedBounds() const {
2834 InvertedBounds inverted_bounds;
2835 {
2837 const std::vector<double> var_lbs,
2838 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_LB, num_gurobi_variables_));
2840 const std::vector<double> var_ubs,
2841 gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_UB, num_gurobi_variables_));
2842 for (const auto& [id, index] : variables_map_) {
2843 if (var_lbs[index] > var_ubs[index]) {
2844 inverted_bounds.variables.push_back(id);
2845 }
2846 }
2847 }
2848 for (const auto& [id, cstr_data] : linear_constraints_map_) {
2849 if (cstr_data.lower_bound > cstr_data.upper_bound) {
2850 inverted_bounds.linear_constraints.push_back(id);
2851 }
2852 }
2853
2854 // Above code have inserted ids in non-stable order.
2855 std::sort(inverted_bounds.variables.begin(), inverted_bounds.variables.end());
2856 std::sort(inverted_bounds.linear_constraints.begin(),
2857 inverted_bounds.linear_constraints.end());
2858 return inverted_bounds;
2859}
2860
2861absl::StatusOr<InvalidIndicators> GurobiSolver::ListInvalidIndicators() const {
2862 InvalidIndicators invalid_indicators;
2863 for (const auto& [constraint_id, indicator_data] :
2864 indicator_constraints_map_) {
2865 if (!indicator_data.has_value()) {
2866 continue;
2867 }
2868 const int64_t indicator_id = indicator_data->indicator_variable_id;
2869 const GurobiVariableIndex variable_index = variables_map_.at(indicator_id);
2870 ASSIGN_OR_RETURN(const double var_lb, gurobi_->GetDoubleAttrElement(
2871 GRB_DBL_ATTR_LB, variable_index));
2872 ASSIGN_OR_RETURN(const double var_ub, gurobi_->GetDoubleAttrElement(
2873 GRB_DBL_ATTR_UB, variable_index));
2875 const char var_type,
2876 gurobi_->GetCharAttrElement(GRB_CHAR_ATTR_VTYPE, variable_index));
2877 if (!(var_type == GRB_BINARY ||
2878 (var_type == GRB_INTEGER && var_lb >= 0.0 && var_ub <= 1.0))) {
2879 invalid_indicators.invalid_indicators.push_back(
2880 {.variable = indicator_id, .constraint = constraint_id});
2881 }
2882 }
2883 // Above code may have inserted ids in non-stable order.
2884 invalid_indicators.Sort();
2885 return invalid_indicators;
2886}
2887
2888bool GurobiSolver::is_multi_objective_mode() const {
2889 return !multi_objectives_map_.empty();
2890}
2891
2892absl::Status GurobiSolver::SetMultiObjectiveParameters(
2893 const ModelSolveParametersProto& model_parameters) {
2894 const auto set_tolerances =
2895 [&](const GurobiMultiObjectiveIndex index,
2896 const ObjectiveParametersProto& objective_parameters)
2897 -> absl::Status {
2898 RETURN_IF_ERROR(gurobi_->SetIntParam("ObjNumber", index));
2899 if (objective_parameters.has_objective_degradation_absolute_tolerance()) {
2900 RETURN_IF_ERROR(gurobi_->SetDoubleAttr(
2901 "ObjNAbsTol",
2902 objective_parameters.objective_degradation_absolute_tolerance()));
2903 }
2904 if (objective_parameters.has_objective_degradation_relative_tolerance()) {
2905 RETURN_IF_ERROR(gurobi_->SetDoubleAttr(
2906 "ObjNRelTol",
2907 objective_parameters.objective_degradation_relative_tolerance()));
2908 }
2909 return absl::OkStatus();
2910 };
2911 const auto set_time_limit =
2912 [&](const GurobiMultiObjectiveIndex index,
2913 const ObjectiveParametersProto& objective_parameters)
2914 -> absl::Status {
2915 if (!objective_parameters.has_time_limit()) {
2916 // Unset time_limit defaults to infinite, so we don't need to do anything.
2917 return absl::OkStatus();
2918 }
2920 const absl::Duration time_limit,
2921 util_time::DecodeGoogleApiProto(objective_parameters.time_limit()));
2922 return gurobi_->SetMultiObjectiveDoubleParam(
2923 GRB_DBL_PAR_TIMELIMIT, index, absl::ToDoubleSeconds(time_limit));
2924 };
2925 if (model_parameters.has_primary_objective_parameters()) {
2926 const GurobiMultiObjectiveIndex obj_index =
2927 multi_objectives_map_.at(std::nullopt);
2928 RETURN_IF_ERROR(set_tolerances(
2929 obj_index, model_parameters.primary_objective_parameters()))
2930 << " for primary objective";
2931 RETURN_IF_ERROR(set_time_limit(
2932 obj_index, model_parameters.primary_objective_parameters()))
2933 << " for primary objective";
2934 }
2935 for (const auto& [id, objective_parameters] :
2936 model_parameters.auxiliary_objective_parameters()) {
2937 const GurobiMultiObjectiveIndex obj_index = multi_objectives_map_.at(id);
2938 RETURN_IF_ERROR(set_tolerances(obj_index, objective_parameters))
2939 << " for auxiliary objective " << id;
2940 RETURN_IF_ERROR(set_time_limit(obj_index, objective_parameters))
2941 << " for auxiliary objective " << id;
2942 }
2943 return absl::OkStatus();
2944}
2945
2946absl::Status GurobiSolver::ResetModelParameters(
2947 const ModelSolveParametersProto& model_parameters) {
2948 for (int i = 0; i < model_parameters.branching_priorities().ids_size(); ++i) {
2949 const int64_t var_id = model_parameters.branching_priorities().ids(i);
2950 const GurobiVariableIndex grb_index = variables_map_.at(var_id);
2952 gurobi_->SetIntAttrElement(GRB_INT_ATTR_BRANCHPRIORITY, grb_index, 0))
2953 << "failed to reset branching priority for variable ID " << var_id
2954 << " (Gurobi index = " << grb_index << ")";
2955 }
2956 for (const int64_t lazy_constraint_id :
2957 model_parameters.lazy_linear_constraint_ids()) {
2958 const GurobiLinearConstraintIndex lazy_constraint_index =
2959 linear_constraints_map_.at(lazy_constraint_id).constraint_index;
2961 gurobi_->SetIntAttrElement(GRB_INT_ATTR_LAZY, lazy_constraint_index, 0))
2962 << "failed to reset lazy constraint for lazy constraint ID "
2963 << lazy_constraint_id << " (Gurobi index = " << lazy_constraint_index
2964 << ")";
2965 }
2966 return absl::OkStatus();
2967}
2968
2969absl::StatusOr<SolveResultProto> GurobiSolver::Solve(
2970 const SolveParametersProto& parameters,
2971 const ModelSolveParametersProto& model_parameters,
2972 const MessageCallback message_cb,
2973 const CallbackRegistrationProto& callback_registration, const Callback cb,
2974 const SolveInterrupter* absl_nullable const interrupter) {
2976 model_parameters, kGurobiSupportedStructures, "Gurobi"));
2977 const absl::Time start = absl::Now();
2978
2979 // Need to run GRBupdatemodel before:
2980 // 1. setting parameters (to test if the problem is a MIP)
2981 // 2. registering callbacks (to test if the problem is a MIP),
2982 // 3. setting basis and getting the obj sense.
2983 // We just run it first.
2984 RETURN_IF_ERROR(gurobi_->UpdateModel());
2985
2986 // Gurobi returns "infeasible" when bounds are inverted.
2987 {
2988 ASSIGN_OR_RETURN(const InvertedBounds inverted_bounds,
2989 ListInvertedBounds());
2990 RETURN_IF_ERROR(inverted_bounds.ToStatus());
2991 }
2992
2993 // Gurobi will silently impose that indicator variables are binary even if
2994 // not so specified by the user in the model. We return an error here if
2995 // this is the case to be consistent across solvers.
2996 {
2997 ASSIGN_OR_RETURN(const InvalidIndicators invalid_indicators,
2998 ListInvalidIndicators());
2999 RETURN_IF_ERROR(invalid_indicators.ToStatus());
3000 }
3001
3002 // We must set the parameters before calling RegisterCallback since it
3003 // changes some parameters depending on the callback registration.
3004 RETURN_IF_ERROR(SetParameters(parameters, model_parameters));
3005
3006 // We use a local interrupter that will triggers the calls to GRBterminate()
3007 // when either the user interrupter is triggered or when a callback returns
3008 // a true `terminate`.
3009 std::unique_ptr<SolveInterrupter> local_interrupter;
3010 if (cb != nullptr || interrupter != nullptr) {
3011 local_interrupter = std::make_unique<SolveInterrupter>();
3012 }
3013 const ScopedSolveInterrupterCallback scoped_terminate_callback(
3014 local_interrupter.get(), [&]() {
3015 // Make an immediate call to GRBterminate() as soon as this
3016 // interrupter is triggered (which may immediately happen in the code
3017 // below when it is chained with the optional user interrupter).
3018 //
3019 // This call may happen too early. This is not an issue since we will
3020 // repeat this call at each call of the Gurobi callback. See the
3021 // comment in GurobiCallbackImpl() for details.
3022 gurobi_->Terminate();
3023 });
3024
3025 // Chain the user interrupter to the local interrupter. If/when the user
3026 // interrupter is triggered, this triggers the local interrupter. This may
3027 // happen immediately if the user interrupter is already triggered.
3028 //
3029 // The local interrupter can also be triggered by a callback returning a
3030 // true `terminate`.
3031 const ScopedSolveInterrupterCallback scoped_chaining_callback(
3032 interrupter, [&]() { local_interrupter->Interrupt(); });
3033
3034 if (model_parameters.has_initial_basis()) {
3035 RETURN_IF_ERROR(SetGurobiBasis(model_parameters.initial_basis()));
3036 }
3037 RETURN_IF_ERROR(gurobi_->SetIntAttr(GRB_INT_ATTR_NUMSTART,
3038 model_parameters.solution_hints_size()));
3039 for (int i = 0; i < model_parameters.solution_hints_size(); ++i) {
3040 RETURN_IF_ERROR(gurobi_->SetIntParam(GRB_INT_PAR_STARTNUMBER, i));
3041 RETURN_IF_ERROR(UpdateDoubleListAttribute(
3042 model_parameters.solution_hints(i).variable_values(),
3043 GRB_DBL_ATTR_START, variables_map_));
3044 }
3046 UpdateInt32ListAttribute(model_parameters.branching_priorities(),
3047 GRB_INT_ATTR_BRANCHPRIORITY, variables_map_));
3048 if (is_multi_objective_mode()) {
3049 RETURN_IF_ERROR(SetMultiObjectiveParameters(model_parameters));
3050 }
3051 for (const int64_t lazy_constraint_id :
3052 model_parameters.lazy_linear_constraint_ids()) {
3053 const GurobiLinearConstraintIndex lazy_constraint_index =
3054 linear_constraints_map_.at(lazy_constraint_id).constraint_index;
3055 // We select a value of "1" here, which means that the lazy constraints will
3056 // be separated at feasible solutions, and that Gurobi has latitude to
3057 // select which violated constraints to add to the model if multiple are
3058 // violated. This seems like a reasonable default choice for us, but we may
3059 // want to revisit later (or expose this choice to the user).
3060 RETURN_IF_ERROR(gurobi_->SetIntAttrElement(GRB_INT_ATTR_LAZY,
3061 lazy_constraint_index, 1));
3062 }
3063
3064 // Here we register the callback when we either have a user callback or a
3065 // local interrupter. The rationale for doing so when we have only an
3066 // interrupter is explained in GurobiCallbackImpl().
3067 Gurobi::Callback grb_cb = nullptr;
3068 std::unique_ptr<GurobiCallbackData> gurobi_cb_data;
3069 if (cb != nullptr || local_interrupter != nullptr || message_cb != nullptr) {
3070 ASSIGN_OR_RETURN(gurobi_cb_data,
3071 RegisterCallback(callback_registration, cb, message_cb,
3072 start, local_interrupter.get()));
3073 grb_cb = [&gurobi_cb_data](
3074 const Gurobi::CallbackContext& cb_context) -> absl::Status {
3075 return GurobiCallbackImpl(cb_context, gurobi_cb_data->callback_input,
3076 gurobi_cb_data->message_callback_data,
3077 gurobi_cb_data->local_interrupter);
3078 };
3079 }
3080
3081 RETURN_IF_ERROR(gurobi_->Optimize(grb_cb));
3082
3083 // We flush message callbacks before testing for Gurobi error in case where
3084 // the unfinished line of message would help with the error.
3085 if (gurobi_cb_data != nullptr) {
3086 GurobiCallbackImplFlush(gurobi_cb_data->callback_input,
3087 gurobi_cb_data->message_callback_data);
3088 }
3089
3091 ExtractSolveResultProto(start, model_parameters));
3092 // Reset Gurobi parameters.
3093 // TODO(b/277246682): ensure that resetting parameters does not degrade
3094 // incrementalism performance.
3095 RETURN_IF_ERROR(gurobi_->ResetParameters());
3096 RETURN_IF_ERROR(ResetModelParameters(model_parameters));
3097
3098 return solve_result;
3099}
3100
3101// TODO(b/277339044): Remove code duplication with GurobiSolver::Solve().
3102absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
3104 const SolveParametersProto& parameters, MessageCallback message_cb,
3105 const SolveInterrupter* absl_nullable const interrupter) {
3106 const absl::Time start = absl::Now();
3107
3108 // Need to run GRBupdatemodel before:
3109 // 1. setting parameters (to test if the problem is a MIP)
3110 // 2. registering callbacks (to test if the problem is a MIP),
3111 RETURN_IF_ERROR(gurobi_->UpdateModel());
3112
3113 // Gurobi will silently impose that indicator variables are binary even if
3114 // not so specified by the user in the model. We return an error here if
3115 // this is the case to be consistent across solvers.
3116 {
3117 ASSIGN_OR_RETURN(const InvalidIndicators invalid_indicators,
3118 ListInvalidIndicators());
3119 RETURN_IF_ERROR(invalid_indicators.ToStatus());
3120 }
3121
3122 // We must set the parameters before calling RegisterCallback since it
3123 // changes some parameters depending on the callback registration.
3124 RETURN_IF_ERROR(SetParameters(parameters));
3125
3126 // We use a local interrupter that will triggers the calls to
3127 // GRBterminate() when either the user interrupter is triggered or when a
3128 // callback returns a true `terminate`.
3129 std::unique_ptr<SolveInterrupter> local_interrupter;
3130 if (interrupter != nullptr) {
3131 local_interrupter = std::make_unique<SolveInterrupter>();
3132 }
3133 const ScopedSolveInterrupterCallback scoped_terminate_callback(
3134 local_interrupter.get(), [&]() {
3135 // Make an immediate call to GRBterminate() as soon as this
3136 // interrupter is triggered (which may immediately happen in the
3137 // code below when it is chained with the optional user
3138 // interrupter).
3139 //
3140 // This call may happen too early. This is not an issue since we
3141 // will repeat this call at each call of the Gurobi callback. See
3142 // the comment in GurobiCallbackImpl() for details.
3143 gurobi_->Terminate();
3144 });
3145
3146 // Chain the user interrupter to the local interrupter. If/when the user
3147 // interrupter is triggered, this triggers the local interrupter. This may
3148 // happen immediately if the user interrupter is already triggered.
3149 //
3150 // The local interrupter can also be triggered by a callback returning a
3151 // true `terminate`.
3152 const ScopedSolveInterrupterCallback scoped_chaining_callback(
3153 interrupter, [&]() { local_interrupter->Interrupt(); });
3154
3155 // Here we register the callback when we either have a message callback or a
3156 // local interrupter. The rationale for doing so when we have only an
3157 // interrupter is explained in GurobiCallbackImpl().
3158 Gurobi::Callback grb_cb = nullptr;
3159 std::unique_ptr<GurobiCallbackData> gurobi_cb_data;
3160 if (local_interrupter != nullptr || message_cb != nullptr) {
3161 ASSIGN_OR_RETURN(gurobi_cb_data,
3162 RegisterCallback({}, nullptr, message_cb, start,
3163 local_interrupter.get()));
3164 grb_cb = [&gurobi_cb_data](
3165 const Gurobi::CallbackContext& cb_context) -> absl::Status {
3166 return GurobiCallbackImpl(cb_context, gurobi_cb_data->callback_input,
3167 gurobi_cb_data->message_callback_data,
3168 gurobi_cb_data->local_interrupter);
3169 };
3170 }
3171
3172 ASSIGN_OR_RETURN(const bool proven_infeasible, gurobi_->ComputeIIS(grb_cb));
3173
3174 // We flush message callbacks before testing for Gurobi error in case
3175 // where the unfinished line of message would help with the error.
3176 if (gurobi_cb_data != nullptr) {
3177 GurobiCallbackImplFlush(gurobi_cb_data->callback_input,
3178 gurobi_cb_data->message_callback_data);
3179 }
3180
3183 ExtractComputeInfeasibleSubsystemResultProto(proven_infeasible));
3184 // Reset Gurobi parameters.
3185 // TODO(b/277246682): ensure that resetting parameters does not degrade
3186 // incrementalism performance.
3187 RETURN_IF_ERROR(gurobi_->ResetParameters());
3188
3189 return iis_result;
3190}
3191
3193
3194} // namespace math_opt
3195} // namespace operations_research
#define ASSIGN_OR_RETURN(lhs, rexpr)
#define RETURN_IF_ERROR(expr)
const ::operations_research::math_opt::SparseVectorFilterProto & mip_solution_filter() const
const ::operations_research::math_opt::SparseVectorFilterProto & mip_node_filter() const
GurobiParametersProto_Parameter Parameter
Definition gurobi.pb.h:684
absl::StatusOr< SolveResultProto > Solve(const SolveParametersProto &parameters, const ModelSolveParametersProto &model_parameters, MessageCallback message_cb, const CallbackRegistrationProto &callback_registration, Callback cb, const SolveInterrupter *absl_nullable interrupter) override
absl::StatusOr< ComputeInfeasibleSubsystemResultProto > ComputeInfeasibleSubsystem(const SolveParametersProto &parameters, MessageCallback message_cb, const SolveInterrupter *absl_nullable interrupter) override
absl::StatusOr< bool > Update(const ModelUpdateProto &model_update) override
static absl::StatusOr< std::unique_ptr< GurobiSolver > > New(const ModelProto &input_model, const SolverInterface::InitArgs &init_args)
static absl::StatusOr< std::unique_ptr< Gurobi > > NewWithSharedPrimaryEnv(GRBenv *primary_env)
Definition g_gurobi.cc:183
std::function< absl::Status(const CallbackContext &)> Callback
Definition g_gurobi.h:222
static absl::StatusOr< std::unique_ptr< Gurobi > > New(GRBenvUniquePtr primary_env=nullptr)
Definition g_gurobi.cc:189
const ::google::protobuf::Map<::int64_t, ::operations_research::math_opt::IndicatorConstraintProto > & new_constraints() const
const ::operations_research::math_opt::ObjectiveProto & objective() const
Definition model.pb.h:4563
const ::google::protobuf::Map<::int64_t, ::operations_research::math_opt::ObjectiveProto > & auxiliary_objectives() const
Definition model.pb.h:4662
const ::operations_research::math_opt::SparseInt32VectorProto & branching_priorities() const
const ::operations_research::math_opt::SolutionHintProto & solution_hints(int index) const
const ::operations_research::math_opt::BasisProto & initial_basis() const
const ::operations_research::math_opt::SosConstraintUpdatesProto & sos1_constraint_updates() const
::int64_t deleted_linear_constraint_ids(int index) const
const ::operations_research::math_opt::SparseDoubleMatrixProto & linear_constraint_matrix_updates() const
const ::operations_research::math_opt::LinearConstraintsProto & new_linear_constraints() const
const ::operations_research::math_opt::AuxiliaryObjectivesUpdatesProto & auxiliary_objectives_updates() const
const ::operations_research::math_opt::VariableUpdatesProto & variable_updates() const
const ::operations_research::math_opt::QuadraticConstraintUpdatesProto & quadratic_constraint_updates() const
const ::operations_research::math_opt::SosConstraintUpdatesProto & sos2_constraint_updates() const
const ::operations_research::math_opt::IndicatorConstraintUpdatesProto & indicator_constraint_updates() const
const ::operations_research::math_opt::ObjectiveUpdatesProto & objective_updates() const
const ::operations_research::math_opt::SecondOrderConeConstraintUpdatesProto & second_order_cone_constraint_updates() const
const ::operations_research::math_opt::LinearConstraintUpdatesProto & linear_constraint_updates() const
const ::operations_research::math_opt::VariablesProto & new_variables() const
const ::operations_research::math_opt::SparseDoubleMatrixProto & quadratic_coefficients() const
Definition model.pb.h:3004
const ::operations_research::math_opt::SparseDoubleMatrixProto & quadratic_coefficients() const
const ::operations_research::math_opt::SparseDoubleVectorProto & linear_coefficients() const
const ::google::protobuf::Map<::int64_t, ::operations_research::math_opt::QuadraticConstraintProto > & new_constraints() const
const ::google::protobuf::Map<::int64_t, ::operations_research::math_opt::SecondOrderConeConstraintProto > & new_constraints() const
const ::operations_research::math_opt::SparseDoubleVectorProto & variable_values() const
std::function< void(const std::vector< std::string > &)> MessageCallback
std::function< absl::StatusOr< CallbackResultProto >( const CallbackDataProto &)> Callback
const ::google::protobuf::Map<::int64_t, ::operations_research::math_opt::SosConstraintProto > & new_constraints() const
const ::operations_research::math_opt::SparseDoubleVectorProto & upper_bounds() const
const ::operations_research::math_opt::SparseDoubleVectorProto & lower_bounds() const
const ::operations_research::math_opt::SparseBoolVectorProto & integers() const
#define GRB_INT_PAR_BARITERLIMIT
#define GRB_SUPERBASIC
#define GRB_INT_PAR_LOGTOCONSOLE
#define GRB_DBL_ATTR_UB
#define GRB_INT_ATTR_BRANCHPRIORITY
#define GRB_DBL_ATTR_START
#define GRB_DBL_PAR_MIPGAP
#define GRB_DBL_PAR_FEASIBILITYTOL
#define GRB_SOLUTION_LIMIT
#define GRB_MAXIMIZE
#define GRB_INT_PAR_SOLUTIONLIMIT
#define GRB_NONBASIC_LOWER
#define GRB_INT_ATTR_MODELSENSE
#define GRB_INT_PAR_CUTS
#define GRB_INFINITY
#define GRB_INT_ATTR_VBASIS
#define GRB_GREATER_EQUAL
#define GRB_INT_PAR_SEED
#define GRB_DBL_ATTR_NODECOUNT
#define GRB_INT_PAR_PRESOLVE
#define GRB_DBL_PAR_MIPGAPABS
#define GRB_DBL_ATTR_ITERCOUNT
#define GRB_INT_PAR_THREADS
#define GRB_DBL_ATTR_BOUND_SVIO
#define GRB_OPTIMAL
#define GRB_DBL_PAR_CUTOFF
#define GRB_INT_ATTR_IIS_QCONSTR
#define GRB_INTEGER
#define GRB_DBL_PAR_ITERATIONLIMIT
#define GRB_INT_PAR_METHOD
#define GRB_INT_ATTR_IS_QP
#define GRB_DBL_ATTR_PI
#define GRB_DBL_ATTR_OBJVAL
#define GRB_NODE_LIMIT
#define GRB_INT_PAR_LAZYCONSTRAINTS
#define GRB_DBL_ATTR_XN
#define GRB_INT_PAR_OBJNUMBER
#define GRB_CONTINUOUS
#define GRB_INT_PAR_SCALEFLAG
#define GRB_DBL_ATTR_OBJ
#define GRB_DBL_PAR_HEURISTICS
#define GRB_METHOD_BARRIER
#define GRB_INT_ATTR_IS_MIP
#define GRB_DBL_ATTR_CONSTR_RESIDUAL
#define GRB_DBL_ATTR_OBJNVAL
#define GRB_SOS_TYPE1
#define GRB_INT_PAR_POOLSOLUTIONS
#define GRB_CHAR_ATTR_VTYPE
#define GRB_INT_ATTR_NUMSTART
#define GRB_DBL_ATTR_BOUND_VIO
#define GRB_TIME_LIMIT
#define GRB_NONBASIC_UPPER
#define GRB_INT_ATTR_IIS_UB
#define GRB_DBL_ATTR_OBJCON
#define GRB_DBL_PAR_BESTBDSTOP
#define GRB_STR_ATTR_MODELNAME
#define GRB_DBL_ATTR_CONSTR_VIO
#define GRB_DBL_ATTR_RC
#define GRB_DBL_ATTR_CONSTR_SRESIDUAL
#define GRB_LOADED
#define GRB_DBL_ATTR_CONSTR_SVIO
#define GRB_INPROGRESS
#define GRB_INT_ATTR_IS_QCP
#define GRB_DBL_PAR_BESTOBJSTOP
#define GRB_INT_ATTR_IIS_GENCONSTR
#define GRB_INF_OR_UNBD
#define GRB_DBL_ATTR_X
#define GRB_SUBOPTIMAL
#define GRB_INFEASIBLE
#define GRB_INT_ATTR_IIS_MINIMAL
#define GRB_MAXINT
#define GRB_INT_PAR_STARTNUMBER
#define GRB_EQUAL
#define GRB_SOS_TYPE2
#define GRB_CHAR_ATTR_QCSENSE
#define GRB_CUTOFF
#define GRB_UNBOUNDED
#define GRB_INT_ATTR_IIS_SOS
#define GRB_DBL_ATTR_QCPI
#define GRB_INT_ATTR_CBASIS
#define GRB_METHOD_DUAL
#define GRB_INT_ATTR_IIS_CONSTR
#define GRB_INT_ATTR_BARITERCOUNT
#define GRB_BASIC
#define GRB_MINIMIZE
#define GRB_INT_ATTR_STATUS
#define GRB_DBL_ATTR_FARKASDUAL
#define GRB_LESS_EQUAL
#define GRB_INTERRUPTED
#define GRB_DBL_ATTR_POOLOBJVAL
#define GRB_DBL_ATTR_LB
#define GRB_INT_PAR_SOLUTIONNUMBER
#define GRB_ITERATION_LIMIT
#define GRB_INT_ATTR_SOLCOUNT
#define GRB_NUMERIC
#define GRB_METHOD_PRIMAL
#define GRB_INT_ATTR_IIS_LB
#define GRB_DBL_PAR_TIMELIMIT
#define GRB_BINARY
#define GRB_DBL_PAR_NODELIMIT
#define GRB_USER_OBJ_LIMIT
#define GRB_DBL_ATTR_UNBDRAY
#define GRB_DBL_ATTR_OBJBOUND
#define GRB_INT_PAR_PRECRUSH
#define GRB_INT_ATTR_LAZY
void InsertOrDie(Collection *const collection, const typename Collection::value_type &value)
Definition map_util.h:159
auto & InsertKeyOrDie(Collection *const collection, const typename Collection::value_type::first_type &key)
Definition map_util.h:178
TerminationProto TerminateForReason(const TerminationReasonProto reason, const absl::string_view detail)
absl::Status ModelIsSupported(const ModelProto &model, const SupportedProblemStructures &support_menu, const absl::string_view solver_name)
std::vector< bool > EventToGurobiWhere(const absl::flat_hash_set< CallbackEventProto > &events)
absl::Status ModelSolveParametersAreSupported(const ModelSolveParametersProto &model_parameters, const SupportedProblemStructures &support_menu, const absl::string_view solver_name)
absl::StatusOr< GRBenvUniquePtr > NewPrimaryEnvironment(std::optional< GurobiInitializerProto::ISVKey > proto_isv_key)
std::function< CallbackResult(const CallbackData &)> Callback
Definition callback.h:94
absl::flat_hash_set< CallbackEventProto > EventSet(const CallbackRegistrationProto &callback_registration)
TerminationProto OptimalTerminationProto(const double finite_primal_objective, const double dual_objective, const absl::string_view detail)
TerminationProto LimitTerminationProto(const bool is_maximize, const LimitProto limit, const std::optional< double > optional_finite_primal_objective, const std::optional< double > optional_dual_objective, const absl::string_view detail)
TerminationProto InfeasibleOrUnboundedTerminationProto(bool is_maximize, const FeasibilityStatusProto dual_feasibility_status, const absl::string_view detail)
bool UpdateIsSupported(const ModelUpdateProto &update, const SupportedProblemStructures &support_menu)
std::function< void(const std::vector< std::string > &)> MessageCallback
ElementId< ElementType::kVariable > VariableId
Definition elements.h:80
void GurobiCallbackImplFlush(const GurobiCallbackInput &callback_input, MessageCallbackData &message_callback_data)
absl::Status CheckRegisteredCallbackEvents(const CallbackRegistrationProto &registration, const absl::flat_hash_set< CallbackEventProto > &supported_events)
SparseVectorView< T > MakeView(absl::Span< const int64_t > ids, const Collection &values)
ElementId< ElementType::kLinearConstraint > LinearConstraintId
Definition elements.h:81
absl::Status GurobiCallbackImpl(const Gurobi::CallbackContext &context, const GurobiCallbackInput &callback_input, MessageCallbackData &message_callback_data, SolveInterrupter *const local_interrupter)
std::unique_ptr< GRBenv, GurobiFreeEnv > GRBenvUniquePtr
Definition g_gurobi.h:59
TerminationProto InfeasibleTerminationProto(bool is_maximize, const FeasibilityStatusProto dual_feasibility_status, const absl::string_view detail)
TerminationProto CutoffTerminationProto(bool is_maximize, const absl::string_view detail)
TerminationProto UnboundedTerminationProto(const bool is_maximize, const absl::string_view detail)
std::vector< K > SortedMapKeys(const absl::flat_hash_map< K, V > &in_map)
Definition sorted.h:55
OR-Tools root namespace.
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
std::string ProtoEnumToString(ProtoEnumType enum_value)
Definition proto_utils.h:63
constexpr bool kAnyXsanEnabled
STL namespace.
inline ::absl::StatusOr< absl::Duration > DecodeGoogleApiProto(const google::protobuf::Duration &proto)
Definition protoutil.h:42
inline ::absl::StatusOr< google::protobuf::Duration > EncodeGoogleApiProto(absl::Duration d)
Definition protoutil.h:27
StatusBuilder InternalErrorBuilder()
StatusBuilder UnimplementedErrorBuilder()
StatusBuilder InvalidArgumentErrorBuilder()
if(!yyg->yy_init)
Definition parser.yy.cc:966
#define MATH_OPT_REGISTER_SOLVER(solver_type, solver_factory)
const NonStreamableGurobiInitArguments * ToNonStreamableGurobiInitArguments() const override