1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the HiGHS linear optimization suite */
4 /* */
5 /* Written and engineered 2008-2021 at the University of Edinburgh */
6 /* */
7 /* Available as open-source under the MIT License */
8 /* */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 #include "mip/HighsMipSolver.h"
11
12 #include "lp_data/HighsLpUtils.h"
13 #include "lp_data/HighsModelUtils.h"
14 #include "mip/HighsCliqueTable.h"
15 #include "mip/HighsCutPool.h"
16 #include "mip/HighsDomain.h"
17 #include "mip/HighsImplications.h"
18 #include "mip/HighsLpRelaxation.h"
19 #include "mip/HighsMipSolverData.h"
20 #include "mip/HighsPseudocost.h"
21 #include "mip/HighsSearch.h"
22 #include "mip/HighsSeparation.h"
23 #include "presolve/PresolveComponent.h"
24 #include "util/HighsCDouble.h"
25
26 #ifdef HIGHS_DEBUGSOL
27 std::vector<double> highsDebugSolution;
28 #endif
29
HighsMipSolver(const HighsOptions & options,const HighsLp & lp,bool submip)30 HighsMipSolver::HighsMipSolver(const HighsOptions& options, const HighsLp& lp,
31 bool submip)
32 : options_mip_(&options), model_(&lp), submip(submip), rootbasis(nullptr) {}
33
34 HighsMipSolver::~HighsMipSolver() = default;
35
runPresolve()36 HighsPresolveStatus HighsMipSolver::runPresolve() {
37 // todo: commented out parts or change Highs::runPresolve to operate on a
38 // parameter LP rather than Highs::lp_. Not sure which approach is preferable.
39 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
40 ML_MINIMAL, "\nrunning MIP presolve\n");
41 const HighsLp& lp_ = *(model_);
42
43 // Exit if the problem is empty or if presolve is set to off.
44 if (options_mip_->presolve == off_string)
45 return HighsPresolveStatus::NotPresolved;
46 if (lp_.numCol_ == 0 && lp_.numRow_ == 0)
47 return HighsPresolveStatus::NullError;
48
49 // Clear info from previous runs if lp_ has been modified.
50 if (presolve_.has_run_) presolve_.clear();
51 double start_presolve = timer_.readRunHighsClock();
52
53 // Set time limit.
54 if (options_mip_->time_limit > 0 &&
55 options_mip_->time_limit < HIGHS_CONST_INF) {
56 double left = options_mip_->time_limit - start_presolve;
57 if (left <= 0) {
58 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
59 ML_VERBOSE,
60 "Time limit reached while reading in matrix\n");
61 return HighsPresolveStatus::Timeout;
62 }
63
64 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
65 ML_VERBOSE,
66 "Time limit set: reading matrix took %.2g, presolve "
67 "time left: %.2g\n",
68 start_presolve, left);
69 presolve_.options_.time_limit = left;
70 }
71
72 // Presolve.
73 presolve_.init(lp_, timer_, true);
74
75 if (options_mip_->time_limit > 0 &&
76 options_mip_->time_limit < HIGHS_CONST_INF) {
77 double current = timer_.readRunHighsClock();
78 double time_init = current - start_presolve;
79 double left = presolve_.options_.time_limit - time_init;
80 if (left <= 0) {
81 HighsPrintMessage(
82 options_mip_->output, options_mip_->message_level, ML_VERBOSE,
83 "Time limit reached while copying matrix into presolve.\n");
84 return HighsPresolveStatus::Timeout;
85 }
86
87 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
88 ML_VERBOSE,
89 "Time limit set: copying matrix took %.2g, presolve "
90 "time left: %.2g\n",
91 time_init, left);
92 presolve_.options_.time_limit = options_mip_->time_limit;
93 }
94
95 presolve_.data_.presolve_[0].message_level = options_mip_->message_level;
96 presolve_.data_.presolve_[0].output = options_mip_->output;
97
98 HighsPresolveStatus presolve_return_status = presolve_.run();
99
100 // Handle max case.
101 if (presolve_return_status == HighsPresolveStatus::Reduced &&
102 lp_.sense_ == ObjSense::MAXIMIZE) {
103 presolve_.negateReducedLpCost();
104 presolve_.data_.reduced_lp_.sense_ = ObjSense::MAXIMIZE;
105 }
106
107 // Update reduction counts.
108 switch (presolve_.presolve_status_) {
109 case HighsPresolveStatus::Reduced: {
110 HighsLp& reduced_lp = presolve_.getReducedProblem();
111 presolve_.info_.n_cols_removed = lp_.numCol_ - reduced_lp.numCol_;
112 presolve_.info_.n_rows_removed = lp_.numRow_ - reduced_lp.numRow_;
113 presolve_.info_.n_nnz_removed =
114 (int)lp_.Avalue_.size() - (int)reduced_lp.Avalue_.size();
115 break;
116 }
117 case HighsPresolveStatus::ReducedToEmpty: {
118 presolve_.info_.n_cols_removed = lp_.numCol_;
119 presolve_.info_.n_rows_removed = lp_.numRow_;
120 presolve_.info_.n_nnz_removed = (int)lp_.Avalue_.size();
121 break;
122 }
123 default:
124 break;
125 }
126 return presolve_return_status;
127 }
128
runPostsolve()129 HighsPostsolveStatus HighsMipSolver::runPostsolve() {
130 assert(presolve_.has_run_);
131 if (!mipdata_) return HighsPostsolveStatus::ReducedSolutionEmpty;
132
133 const std::vector<double>& incumbent = mipdata_->getSolution();
134 bool solution_ok =
135 presolve_.getReducedProblem().numCol_ == (int)incumbent.size();
136 if (!solution_ok) return HighsPostsolveStatus::ReducedSolutionDimenionsError;
137
138 if (presolve_.presolve_status_ != HighsPresolveStatus::Reduced &&
139 presolve_.presolve_status_ != HighsPresolveStatus::ReducedToEmpty)
140 return HighsPostsolveStatus::NoPostsolve;
141
142 // Handle max case.
143 // if (lp_.sense_ == ObjSense::MAXIMIZE)
144 // presolve_.negateReducedLpColDuals(true);
145
146 // Run postsolve
147 HighsPostsolveStatus postsolve_status =
148 presolve_.data_.presolve_[0].primalPostsolve(
149 incumbent, presolve_.data_.recovered_solution_);
150
151 return postsolve_status;
152
153 // if (lp_.sense_ == ObjSense::MAXIMIZE)
154 // presolve_.negateReducedLpColDuals(false);
155
156 // return HighsPostsolveStatus::SolutionRecovered;
157 }
158
run()159 void HighsMipSolver::run() {
160 modelstatus_ = HighsModelStatus::NOTSET;
161 // std::cout << options_mip_->presolve << std::endl;
162 timer_.start(timer_.solve_clock);
163 if (options_mip_->presolve != "off") {
164 HighsPresolveStatus presolve_status = runPresolve();
165 switch (presolve_status) {
166 case HighsPresolveStatus::Reduced:
167 reportPresolveReductions(*options_mip_, *model_,
168 presolve_.getReducedProblem());
169 model_ = &presolve_.getReducedProblem();
170 break;
171 case HighsPresolveStatus::Unbounded:
172 modelstatus_ = HighsModelStatus::PRIMAL_UNBOUNDED;
173 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
174 ML_MINIMAL,
175 "Presolve: Model detected to be unbounded\n");
176 timer_.stop(timer_.solve_clock);
177 return;
178 case HighsPresolveStatus::Infeasible:
179 modelstatus_ = HighsModelStatus::PRIMAL_INFEASIBLE;
180 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
181 ML_MINIMAL,
182 "Presolve: Model detected to be infeasible\n");
183 timer_.stop(timer_.solve_clock);
184 return;
185 case HighsPresolveStatus::Timeout:
186 modelstatus_ = HighsModelStatus::REACHED_TIME_LIMIT;
187 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
188 ML_MINIMAL, "Time limit reached during presolve\n");
189 timer_.stop(timer_.solve_clock);
190 return;
191 case HighsPresolveStatus::ReducedToEmpty:
192 modelstatus_ = HighsModelStatus::OPTIMAL;
193 reportPresolveReductions(*options_mip_, *model_, true);
194 mipdata_ = decltype(mipdata_)(new HighsMipSolverData(*this));
195 mipdata_->init();
196 runPostsolve();
197 timer_.stop(timer_.solve_clock);
198 return;
199 case HighsPresolveStatus::NotReduced:
200 reportPresolveReductions(*options_mip_, *model_, false);
201 break;
202 default:
203 assert(false);
204 }
205 }
206
207 mipdata_ = decltype(mipdata_)(new HighsMipSolverData(*this));
208 mipdata_->init();
209 mipdata_->runSetup();
210 if (modelstatus_ == HighsModelStatus::NOTSET) {
211 mipdata_->evaluateRootNode();
212 }
213 if (mipdata_->nodequeue.empty()) {
214 mipdata_->printDisplayLine();
215 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
216 ML_MINIMAL, "\nSolving stopped with status: %s\n",
217 utilHighsModelStatusToString(modelstatus_).c_str());
218 bool haveSolution = mipdata_->upper_bound != HIGHS_CONST_INF;
219 if (mipdata_->modelcleanup) {
220 model_ = mipdata_->modelcleanup->origmodel;
221 if (haveSolution)
222 mipdata_->modelcleanup->recoverSolution(mipdata_->incumbent);
223 }
224 if (haveSolution) {
225 if (options_mip_->presolve != "off")
226 runPostsolve();
227 else if (!mipdata_->getSolution().empty()) {
228 presolve_.data_.recovered_solution_.col_value = mipdata_->getSolution();
229 calculateRowValues(*model_, presolve_.data_.recovered_solution_);
230 }
231 }
232 timer_.stop(timer_.solve_clock);
233 return;
234 }
235
236 std::shared_ptr<const HighsBasis> basis;
237 HighsSearch search{*this, mipdata_->pseudocost};
238 HighsSeparation sepa;
239
240 search.setLpRelaxation(&mipdata_->lp);
241 sepa.setLpRelaxation(&mipdata_->lp);
242
243 mipdata_->lower_bound = mipdata_->nodequeue.getBestLowerBound();
244 search.installNode(mipdata_->nodequeue.popBestBoundNode());
245
246 mipdata_->printDisplayLine();
247
248 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
249 ML_MINIMAL, "\nstarting tree search\n");
250
251 while (search.hasNode()) {
252 // set iteration limit for each lp solve during the dive to 10 times the
253 // average nodes
254
255 int iterlimit = 100 * int(mipdata_->lp.getNumLpIterations() /
256 (double)std::max(size_t{1}, mipdata_->num_nodes));
257 iterlimit = std::max(1000, iterlimit);
258
259 mipdata_->lp.setIterationLimit(iterlimit);
260
261 // perform the dive and put the open nodes to the queue
262 size_t plungestart = mipdata_->num_nodes;
263 bool limit_reached = false;
264 while (true) {
265 if (mipdata_->heuristic_lp_iterations <
266 mipdata_->total_lp_iterations * mipdata_->heuristic_effort) {
267 search.evaluateNode();
268 if (!search.currentNodePruned()) search.heuristicSearchNew();
269 }
270
271 search.dive();
272 ++mipdata_->num_leaves;
273
274 search.flushStatistics();
275 if (mipdata_->checkLimits()) {
276 limit_reached = true;
277 break;
278 }
279
280 if (!search.backtrack()) break;
281
282 if (search.getCurrentEstimate() >= mipdata_->upper_limit) break;
283
284 if (mipdata_->num_nodes - plungestart >=
285 std::min(size_t{1000}, mipdata_->num_nodes / 10))
286 break;
287
288 if (mipdata_->dispfreq != 0) {
289 if (mipdata_->num_leaves - mipdata_->last_displeave >=
290 mipdata_->dispfreq)
291 mipdata_->printDisplayLine();
292 }
293
294 // printf("continue plunging due to good esitmate\n");
295 }
296 search.openNodesToQueue(mipdata_->nodequeue);
297 mipdata_->lower_bound = std::min(mipdata_->upper_bound,
298 mipdata_->nodequeue.getBestLowerBound());
299
300 if (limit_reached) break;
301
302 if (mipdata_->dispfreq != 0) {
303 if (mipdata_->num_leaves - mipdata_->last_displeave >= mipdata_->dispfreq)
304 mipdata_->printDisplayLine();
305 }
306
307 // the search datastructure should have no installed node now
308 assert(!search.hasNode());
309
310 // propagate the global domain
311 mipdata_->domain.propagate();
312
313 #ifdef HIGHS_DEBUGSOL
314 assert(!mipdata_->domain.infeasible());
315 for (int i = 0; i != numCol(); ++i) {
316 assert(highsDebugSolution[i] + mipdata_->epsilon >=
317 mipdata_->domain.colLower_[i]);
318 assert(highsDebugSolution[i] - mipdata_->epsilon <=
319 mipdata_->domain.colUpper_[i]);
320 }
321 #endif
322
323 // if global propagation detected infeasibility, stop here
324 if (mipdata_->domain.infeasible()) {
325 mipdata_->nodequeue.clear();
326 mipdata_->pruned_treeweight = 1.0;
327 break;
328 }
329
330 // if global propagation found bound changes, we update the local domain
331 if (!mipdata_->domain.getChangedCols().empty()) {
332 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
333 ML_MINIMAL, "added %d global bound changes\n",
334 (int)mipdata_->domain.getChangedCols().size());
335 mipdata_->cliquetable.cleanupFixed(mipdata_->domain);
336 mipdata_->domain.setDomainChangeStack(std::vector<HighsDomainChange>());
337 search.resetLocalDomain();
338
339 mipdata_->domain.clearChangedCols();
340 }
341
342 // remove the iteration limit when installing a new node
343 mipdata_->lp.setIterationLimit();
344
345 // loop to install the next node for the search
346 while (!mipdata_->nodequeue.empty()) {
347 // printf("popping node from nodequeue (length = %d)\n",
348 // (int)nodequeue.size());
349 assert(!search.hasNode());
350 search.installNode(mipdata_->nodequeue.popBestNode());
351
352 assert(search.hasNode());
353
354 // set the current basis if available
355 if (basis) {
356 mipdata_->lp.setStoredBasis(basis);
357 mipdata_->lp.recoverBasis();
358 }
359
360 // we evaluate the node directly here instead of performing a dive
361 // because we first want to check if the node is not fathomed due to
362 // new global information before we perform separation rounds for the node
363 search.evaluateNode();
364
365 // if the node was pruned we remove it from the search and install the
366 // next node from the queue
367 if (search.currentNodePruned()) {
368 search.backtrack();
369 ++mipdata_->num_leaves;
370 ++mipdata_->num_nodes;
371 search.flushStatistics();
372
373 if (mipdata_->checkLimits()) {
374 limit_reached = true;
375 break;
376 }
377
378 mipdata_->lower_bound = std::min(
379 mipdata_->upper_bound, mipdata_->nodequeue.getBestLowerBound());
380
381 if (mipdata_->dispfreq != 0) {
382 if (mipdata_->num_leaves - mipdata_->last_displeave >=
383 mipdata_->dispfreq)
384 mipdata_->printDisplayLine();
385 }
386 continue;
387 }
388
389 // the node is still not fathomed, so perform separation
390 sepa.separate(search.getLocalDomain());
391
392 // after separation we store the new basis and proceed with the outer loop
393 // to perform a dive from this node
394 if (mipdata_->lp.getStatus() != HighsLpRelaxation::Status::Error &&
395 mipdata_->lp.getStatus() != HighsLpRelaxation::Status::NotSet)
396 mipdata_->lp.storeBasis();
397
398 basis = mipdata_->lp.getStoredBasis();
399
400 break;
401 }
402
403 if (limit_reached) break;
404 }
405
406 mipdata_->printDisplayLine();
407 bool havesolution = mipdata_->upper_bound != HIGHS_CONST_INF;
408
409 if (modelstatus_ == HighsModelStatus::NOTSET) {
410 if (havesolution)
411 modelstatus_ = HighsModelStatus::OPTIMAL;
412 else
413 modelstatus_ = HighsModelStatus::PRIMAL_INFEASIBLE;
414 }
415
416 HighsPrintMessage(options_mip_->output, options_mip_->message_level,
417 ML_MINIMAL, "\nSolving stopped with status: %s\n",
418 utilHighsModelStatusToString(modelstatus_).c_str());
419
420 if (mipdata_->modelcleanup) {
421 model_ = mipdata_->modelcleanup->origmodel;
422 if (havesolution)
423 mipdata_->modelcleanup->recoverSolution(mipdata_->incumbent);
424 }
425
426 if (havesolution) {
427 if (options_mip_->presolve != "off")
428 runPostsolve();
429 else if (!mipdata_->getSolution().empty()) {
430 presolve_.data_.recovered_solution_.col_value = mipdata_->getSolution();
431 calculateRowValues(*model_, presolve_.data_.recovered_solution_);
432 }
433 }
434
435 timer_.stop(timer_.solve_clock);
436
437 assert(modelstatus_ != HighsModelStatus::NOTSET);
438 }