1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 //
15
16 #include <cstdint>
17 #include <limits>
18 #include <memory>
19 #include <string>
20 #include <utility>
21 #include <vector>
22
23 #include "absl/base/attributes.h"
24 #include "absl/status/status.h"
25 #include "absl/strings/str_format.h"
26 #include "ortools/base/commandlineflags.h"
27 #include "ortools/base/hash.h"
28 #include "ortools/base/integral_types.h"
29 #include "ortools/base/logging.h"
30 #include "ortools/base/timer.h"
31 #include "ortools/linear_solver/linear_solver.h"
32
33 #if defined(USE_CBC)
34
35 #undef PACKAGE
36 #undef VERSION
37 #include "CbcConfig.h"
38 #include "CbcMessage.hpp"
39 #include "CbcModel.hpp"
40 #include "CoinModel.hpp"
41 #include "OsiClpSolverInterface.hpp"
42
43 // Heuristics
44
45 namespace operations_research {
46
47 class CBCInterface : public MPSolverInterface {
48 public:
49 // Constructor that takes a name for the underlying glpk solver.
50 explicit CBCInterface(MPSolver* const solver);
51 ~CBCInterface() override;
52
53 // ----- Reset -----
54 void Reset() override;
55
56 // Sets the optimization direction (min/max).
57 void SetOptimizationDirection(bool maximize) override;
58
59 // ----- Parameters -----
60
SetNumThreads(int num_threads)61 absl::Status SetNumThreads(int num_threads) override {
62 CHECK_GE(num_threads, 1);
63 num_threads_ = num_threads;
64 return absl::OkStatus();
65 }
66
67 // ----- Solve -----
68 // Solve the problem using the parameter values specified.
69 MPSolver::ResultStatus Solve(const MPSolverParameters& param) override;
70
71 // TODO(user): separate the solve from the model extraction.
ExtractModel()72 virtual void ExtractModel() {}
73
74 // Query problem type.
IsContinuous() const75 bool IsContinuous() const override { return false; }
IsLP() const76 bool IsLP() const override { return false; }
IsMIP() const77 bool IsMIP() const override { return true; }
78
79 // Modify bounds.
80 void SetVariableBounds(int var_index, double lb, double ub) override;
81 void SetVariableInteger(int var_index, bool integer) override;
82 void SetConstraintBounds(int row_index, double lb, double ub) override;
83
84 // Add constraint incrementally.
85 void AddRowConstraint(MPConstraint* const ct) override;
86 // Add variable incrementally.
87 void AddVariable(MPVariable* const var) override;
88 // Change a coefficient in a constraint.
SetCoefficient(MPConstraint * const constraint,const MPVariable * const variable,double new_value,double old_value)89 void SetCoefficient(MPConstraint* const constraint,
90 const MPVariable* const variable, double new_value,
91 double old_value) override {
92 sync_status_ = MUST_RELOAD;
93 }
94 // Clear a constraint from all its terms.
ClearConstraint(MPConstraint * const constraint)95 void ClearConstraint(MPConstraint* const constraint) override {
96 sync_status_ = MUST_RELOAD;
97 }
98
99 // Change a coefficient in the linear objective.
SetObjectiveCoefficient(const MPVariable * const variable,double coefficient)100 void SetObjectiveCoefficient(const MPVariable* const variable,
101 double coefficient) override {
102 sync_status_ = MUST_RELOAD;
103 }
104 // Change the constant term in the linear objective.
SetObjectiveOffset(double value)105 void SetObjectiveOffset(double value) override { sync_status_ = MUST_RELOAD; }
106 // Clear the objective from all its terms.
ClearObjective()107 void ClearObjective() override { sync_status_ = MUST_RELOAD; }
108
109 // Number of simplex iterations
110 int64_t iterations() const override;
111 // Number of branch-and-bound nodes. Only available for discrete problems.
112 int64_t nodes() const override;
113
114 // Returns the basis status of a row.
row_status(int constraint_index) const115 MPSolver::BasisStatus row_status(int constraint_index) const override {
116 LOG(FATAL) << "Basis status only available for continuous problems";
117 return MPSolver::FREE;
118 }
119 // Returns the basis status of a column.
column_status(int variable_index) const120 MPSolver::BasisStatus column_status(int variable_index) const override {
121 LOG(FATAL) << "Basis status only available for continuous problems";
122 return MPSolver::FREE;
123 }
124
ExtractNewVariables()125 void ExtractNewVariables() override {}
ExtractNewConstraints()126 void ExtractNewConstraints() override {}
ExtractObjective()127 void ExtractObjective() override {}
128
SolverVersion() const129 std::string SolverVersion() const override { return "Cbc " CBC_VERSION; }
130
131 // TODO(user): Maybe we should expose the CbcModel build from osi_
132 // instead, but a new CbcModel is built every time Solve is called,
133 // so it is not possible right now.
underlying_solver()134 void* underlying_solver() override { return reinterpret_cast<void*>(&osi_); }
135
136 private:
137 // Reset best objective bound to +/- infinity depending on the
138 // optimization direction.
139 void ResetBestObjectiveBound();
140
141 // Set all parameters in the underlying solver.
142 void SetParameters(const MPSolverParameters& param) override;
143 // Set each parameter in the underlying solver.
144 void SetRelativeMipGap(double value) override;
145 void SetPrimalTolerance(double value) override;
146 void SetDualTolerance(double value) override;
147 void SetPresolveMode(int value) override;
148 void SetScalingMode(int value) override;
149 void SetLpAlgorithm(int value) override;
150
151 OsiClpSolverInterface osi_;
152 // TODO(user): remove and query number of iterations directly from CbcModel
153 int64_t iterations_;
154 int64_t nodes_;
155 // Special way to handle the relative MIP gap parameter.
156 double relative_mip_gap_;
157 int num_threads_ = 1;
158 };
159
160 // ----- Solver -----
161
162 // Creates a LP/MIP instance with the specified name and minimization objective.
CBCInterface(MPSolver * const solver)163 CBCInterface::CBCInterface(MPSolver* const solver)
164 : MPSolverInterface(solver),
165 iterations_(0),
166 nodes_(0),
167 relative_mip_gap_(MPSolverParameters::kDefaultRelativeMipGap) {
168 osi_.setStrParam(OsiProbName, solver_->name_);
169 osi_.setObjSense(1);
170 }
171
~CBCInterface()172 CBCInterface::~CBCInterface() {}
173
174 // Reset the solver.
Reset()175 void CBCInterface::Reset() {
176 osi_.reset();
177 osi_.setObjSense(maximize_ ? -1 : 1);
178 osi_.setStrParam(OsiProbName, solver_->name_);
179 ResetExtractionInformation();
180 }
181
ResetBestObjectiveBound()182 void CBCInterface::ResetBestObjectiveBound() {
183 if (maximize_) {
184 best_objective_bound_ = std::numeric_limits<double>::infinity();
185 } else {
186 best_objective_bound_ = -std::numeric_limits<double>::infinity();
187 }
188 }
189
SetOptimizationDirection(bool maximize)190 void CBCInterface::SetOptimizationDirection(bool maximize) {
191 InvalidateSolutionSynchronization();
192 if (sync_status_ == MODEL_SYNCHRONIZED) {
193 osi_.setObjSense(maximize ? -1 : 1);
194 } else {
195 sync_status_ = MUST_RELOAD;
196 }
197 }
198
199 namespace {
200 // CBC adds a "dummy" variable with index 0 to represent the objective offset.
MPSolverVarIndexToCbcVarIndex(int var_index)201 int MPSolverVarIndexToCbcVarIndex(int var_index) { return var_index + 1; }
202 } // namespace
203
SetVariableBounds(int var_index,double lb,double ub)204 void CBCInterface::SetVariableBounds(int var_index, double lb, double ub) {
205 InvalidateSolutionSynchronization();
206 if (sync_status_ == MODEL_SYNCHRONIZED) {
207 osi_.setColBounds(MPSolverVarIndexToCbcVarIndex(var_index), lb, ub);
208 } else {
209 sync_status_ = MUST_RELOAD;
210 }
211 }
212
SetVariableInteger(int var_index,bool integer)213 void CBCInterface::SetVariableInteger(int var_index, bool integer) {
214 InvalidateSolutionSynchronization();
215 // TODO(user) : Check if this is actually a change.
216 if (sync_status_ == MODEL_SYNCHRONIZED) {
217 if (integer) {
218 osi_.setInteger(MPSolverVarIndexToCbcVarIndex(var_index));
219 } else {
220 osi_.setContinuous(MPSolverVarIndexToCbcVarIndex(var_index));
221 }
222 } else {
223 sync_status_ = MUST_RELOAD;
224 }
225 }
226
SetConstraintBounds(int index,double lb,double ub)227 void CBCInterface::SetConstraintBounds(int index, double lb, double ub) {
228 InvalidateSolutionSynchronization();
229 if (sync_status_ == MODEL_SYNCHRONIZED) {
230 osi_.setRowBounds(index, lb, ub);
231 } else {
232 sync_status_ = MUST_RELOAD;
233 }
234 }
235
AddRowConstraint(MPConstraint * const ct)236 void CBCInterface::AddRowConstraint(MPConstraint* const ct) {
237 sync_status_ = MUST_RELOAD;
238 }
239
AddVariable(MPVariable * const var)240 void CBCInterface::AddVariable(MPVariable* const var) {
241 sync_status_ = MUST_RELOAD;
242 }
243
244 // Solve the LP/MIP. Returns true only if the optimal solution was revealed.
245 // Returns the status of the search.
Solve(const MPSolverParameters & param)246 MPSolver::ResultStatus CBCInterface::Solve(const MPSolverParameters& param) {
247 // CBC requires unique variable and constraint names. By using Lookup*, we
248 // generate variable and constraint indices and ensure the duplicate name
249 // crash will happen here with a readable error message.
250 if (!solver_->variables_.empty()) {
251 solver_->LookupVariableOrNull(solver_->variables_[0]->name());
252 }
253 if (!solver_->constraints_.empty()) {
254 solver_->LookupConstraintOrNull(solver_->constraints_[0]->name());
255 }
256
257 WallTimer timer;
258 timer.Start();
259
260 // Note that CBC does not provide any incrementality.
261 if (param.GetIntegerParam(MPSolverParameters::INCREMENTALITY) ==
262 MPSolverParameters::INCREMENTALITY_OFF) {
263 Reset();
264 }
265
266 // Special case if the model is empty since CBC is not able to
267 // handle this special case by itself.
268 if (solver_->variables_.empty() && solver_->constraints_.empty()) {
269 sync_status_ = SOLUTION_SYNCHRONIZED;
270 result_status_ = MPSolver::OPTIMAL;
271 objective_value_ = solver_->Objective().offset();
272 best_objective_bound_ = solver_->Objective().offset();
273 return result_status_;
274 }
275
276 // Finish preparing the problem.
277 // Define variables.
278 switch (sync_status_) {
279 case MUST_RELOAD: {
280 Reset();
281 CoinModel build;
282 // Create dummy variable for objective offset.
283 build.addColumn(0, nullptr, nullptr, 1.0, 1.0,
284 solver_->Objective().offset(), "dummy", false);
285 const int nb_vars = solver_->variables_.size();
286 for (int i = 0; i < nb_vars; ++i) {
287 MPVariable* const var = solver_->variables_[i];
288 set_variable_as_extracted(i, true);
289 const double obj_coeff = solver_->Objective().GetCoefficient(var);
290 if (var->name().empty()) {
291 build.addColumn(0, nullptr, nullptr, var->lb(), var->ub(), obj_coeff,
292 nullptr, var->integer());
293 } else {
294 build.addColumn(0, nullptr, nullptr, var->lb(), var->ub(), obj_coeff,
295 var->name().c_str(), var->integer());
296 }
297 }
298
299 // Define constraints.
300 int max_row_length = 0;
301 for (int i = 0; i < solver_->constraints_.size(); ++i) {
302 MPConstraint* const ct = solver_->constraints_[i];
303 set_constraint_as_extracted(i, true);
304 if (ct->coefficients_.size() > max_row_length) {
305 max_row_length = ct->coefficients_.size();
306 }
307 }
308 std::unique_ptr<int[]> indices(new int[max_row_length]);
309 std::unique_ptr<double[]> coefs(new double[max_row_length]);
310
311 for (int i = 0; i < solver_->constraints_.size(); ++i) {
312 MPConstraint* const ct = solver_->constraints_[i];
313 const int size = ct->coefficients_.size();
314 int j = 0;
315 for (const auto& entry : ct->coefficients_) {
316 const int index = MPSolverVarIndexToCbcVarIndex(entry.first->index());
317 indices[j] = index;
318 coefs[j] = entry.second;
319 j++;
320 }
321 if (ct->name().empty()) {
322 build.addRow(size, indices.get(), coefs.get(), ct->lb(), ct->ub());
323 } else {
324 build.addRow(size, indices.get(), coefs.get(), ct->lb(), ct->ub(),
325 ct->name().c_str());
326 }
327 }
328 osi_.loadFromCoinModel(build);
329 break;
330 }
331 case MODEL_SYNCHRONIZED: {
332 break;
333 }
334 case SOLUTION_SYNCHRONIZED: {
335 break;
336 }
337 }
338
339 // Changing optimization direction through OSI so that the model file
340 // (written through OSI) has the correct optimization duration.
341 osi_.setObjSense(maximize_ ? -1 : 1);
342
343 sync_status_ = MODEL_SYNCHRONIZED;
344 VLOG(1) << absl::StrFormat("Model built in %.3f seconds.", timer.Get());
345
346 ResetBestObjectiveBound();
347
348 // Solve
349 CbcModel model(osi_);
350
351 // Set log level.
352 CoinMessageHandler message_handler;
353 model.passInMessageHandler(&message_handler);
354 if (quiet_) {
355 message_handler.setLogLevel(0, 0); // Coin messages
356 message_handler.setLogLevel(1, 0); // Clp messages
357 message_handler.setLogLevel(2, 0); // Presolve messages
358 message_handler.setLogLevel(3, 0); // Cgl messages
359 } else {
360 message_handler.setLogLevel(0, 1); // Coin messages
361 message_handler.setLogLevel(1, 1); // Clp messages
362 message_handler.setLogLevel(2, 1); // Presolve messages
363 message_handler.setLogLevel(3, 1); // Cgl messages
364 }
365
366 // Time limit.
367 if (solver_->time_limit() != 0) {
368 VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
369 model.setMaximumSeconds(solver_->time_limit_in_secs());
370 }
371
372 // And solve.
373 timer.Restart();
374
375 // Here we use the default function from the command-line CBC solver.
376 // This enables to activate all the features and get the same performance
377 // as the CBC stand-alone executable. The syntax is ugly, however.
378 SetParameters(param);
379 // Always turn presolve on (it's the CBC default and it consistently
380 // improves performance).
381 model.setTypePresolve(0);
382 // Special way to set the relative MIP gap parameter as it cannot be set
383 // through callCbc.
384 model.setAllowableFractionGap(relative_mip_gap_);
385 // NOTE: Trailing space is required to avoid buffer overflow in cbc.
386 int return_status =
387 num_threads_ == 1
388 ? callCbc("-solve ", model)
389 : callCbc(absl::StrCat("-threads ", num_threads_, " -solve "), model);
390 const int kBadReturnStatus = 777;
391 CHECK_NE(kBadReturnStatus, return_status); // Should never happen according
392 // to the CBC source
393
394 VLOG(1) << absl::StrFormat("Solved in %.3f seconds.", timer.Get());
395
396 // Check the status: optimal, infeasible, etc.
397 int tmp_status = model.status();
398
399 VLOG(1) << "cbc result status: " << tmp_status;
400 /* Final status of problem
401 (info from cbc/.../CbcSolver.cpp,
402 See http://cs?q="cbc+status"+file:CbcSolver.cpp)
403 Some of these can be found out by is...... functions
404 -1 before branchAndBound
405 0 finished - check isProvenOptimal or isProvenInfeasible to see
406 if solution found
407 (or check value of best solution)
408 1 stopped - on maxnodes, maxsols, maxtime
409 2 difficulties so run was abandoned
410 (5 event user programmed event occurred)
411 */
412 switch (tmp_status) {
413 case 0:
414 // Order of tests counts; if model.isContinuousUnbounded() returns true,
415 // then so does model.isProvenInfeasible()!
416 if (model.isProvenOptimal()) {
417 result_status_ = MPSolver::OPTIMAL;
418 } else if (model.isContinuousUnbounded()) {
419 result_status_ = MPSolver::UNBOUNDED;
420 } else if (model.isProvenInfeasible()) {
421 result_status_ = MPSolver::INFEASIBLE;
422 } else if (model.isAbandoned()) {
423 result_status_ = MPSolver::ABNORMAL;
424 } else {
425 result_status_ = MPSolver::ABNORMAL;
426 }
427 break;
428 case 1:
429 if (model.bestSolution() != nullptr) {
430 result_status_ = MPSolver::FEASIBLE;
431 } else {
432 result_status_ = MPSolver::NOT_SOLVED;
433 }
434 break;
435 default:
436 result_status_ = MPSolver::ABNORMAL;
437 break;
438 }
439
440 if (result_status_ == MPSolver::OPTIMAL ||
441 result_status_ == MPSolver::FEASIBLE) {
442 // Get the results
443 objective_value_ = model.getObjValue();
444 VLOG(1) << "objective=" << objective_value_;
445 const double* const values = model.bestSolution();
446 if (values != nullptr) {
447 // if optimal or feasible solution is found.
448 for (int i = 0; i < solver_->variables_.size(); ++i) {
449 MPVariable* const var = solver_->variables_[i];
450 const int var_index = MPSolverVarIndexToCbcVarIndex(var->index());
451 const double val = values[var_index];
452 var->set_solution_value(val);
453 VLOG(3) << var->name() << "=" << val;
454 }
455 } else {
456 VLOG(1) << "No feasible solution found.";
457 }
458 }
459
460 iterations_ = model.getIterationCount();
461 nodes_ = model.getNodeCount();
462 best_objective_bound_ = model.getBestPossibleObjValue();
463 VLOG(1) << "best objective bound=" << best_objective_bound_;
464
465 sync_status_ = SOLUTION_SYNCHRONIZED;
466
467 return result_status_;
468 }
469
470 // ------ Query statistics on the solution and the solve ------
471
iterations() const472 int64_t CBCInterface::iterations() const {
473 if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes;
474 return iterations_;
475 }
476
nodes() const477 int64_t CBCInterface::nodes() const {
478 if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations;
479 return nodes_;
480 }
481
482 // ----- Parameters -----
483
484 // The support for parameters in CBC is intentionally sparse. There is
485 // a memory leak in callCbc that prevents to pass parameters through
486 // it, so handling parameters would require an comprehensive rewrite
487 // of the code. I will improve the parameter support only if there is
488 // a relevant use case.
489
SetParameters(const MPSolverParameters & param)490 void CBCInterface::SetParameters(const MPSolverParameters& param) {
491 SetCommonParameters(param);
492 SetMIPParameters(param);
493 }
494
SetRelativeMipGap(double value)495 void CBCInterface::SetRelativeMipGap(double value) {
496 relative_mip_gap_ = value;
497 }
498
SetPrimalTolerance(double value)499 void CBCInterface::SetPrimalTolerance(double value) {
500 // Skip the warning for the default value as it coincides with
501 // the default value in CBC.
502 if (value != MPSolverParameters::kDefaultPrimalTolerance) {
503 SetUnsupportedDoubleParam(MPSolverParameters::PRIMAL_TOLERANCE);
504 }
505 }
506
SetDualTolerance(double value)507 void CBCInterface::SetDualTolerance(double value) {
508 // Skip the warning for the default value as it coincides with
509 // the default value in CBC.
510 if (value != MPSolverParameters::kDefaultDualTolerance) {
511 SetUnsupportedDoubleParam(MPSolverParameters::DUAL_TOLERANCE);
512 }
513 }
514
SetPresolveMode(int value)515 void CBCInterface::SetPresolveMode(int value) {
516 switch (value) {
517 case MPSolverParameters::PRESOLVE_ON: {
518 // CBC presolve is always on.
519 break;
520 }
521 default: {
522 SetUnsupportedIntegerParam(MPSolverParameters::PRESOLVE);
523 }
524 }
525 }
526
SetScalingMode(int value)527 void CBCInterface::SetScalingMode(int value) {
528 SetUnsupportedIntegerParam(MPSolverParameters::SCALING);
529 }
530
SetLpAlgorithm(int value)531 void CBCInterface::SetLpAlgorithm(int value) {
532 SetUnsupportedIntegerParam(MPSolverParameters::LP_ALGORITHM);
533 }
534
BuildCBCInterface(MPSolver * const solver)535 MPSolverInterface* BuildCBCInterface(MPSolver* const solver) {
536 return new CBCInterface(solver);
537 }
538
539 } // namespace operations_research
540 #endif // #if defined(USE_CBC)
541