Merge pull request #1992 from ERGO-Code/fix-1991

Improve performance and numerics when repairing infeasibilities in transformed integer solution

src/mip/HighsMipSolver.cpp

@@ -643,6 +643,7 @@ void HighsMipSolver::cleanupSolve() {
                "                    %.2f (presolve)\n"
                "                    %.2f (postsolve)\n"
                "  Nodes             %llu\n"
+               "  Repair LPs        %llu (%llu feasible; %llu iterations)\n"
                "  LP iterations     %llu (total)\n"
                "                    %llu (strong br.)\n"
                "                    %llu (separation)\n"
@@ -651,6 +652,9 @@ void HighsMipSolver::cleanupSolve() {
                timer_.read(timer_.presolve_clock),
                timer_.read(timer_.postsolve_clock),
                (long long unsigned)mipdata_->num_nodes,
+               (long long unsigned)mipdata_->total_repair_lp,
+               (long long unsigned)mipdata_->total_repair_lp_feasible,
+               (long long unsigned)mipdata_->total_repair_lp_iterations,
                (long long unsigned)mipdata_->total_lp_iterations,
                (long long unsigned)mipdata_->sb_lp_iterations,
                (long long unsigned)mipdata_->sepa_lp_iterations,

src/mip/HighsMipSolverData.cpp

@@ -366,6 +366,9 @@ void HighsMipSolverData::init() {
   num_nodes_before_run = 0;
   num_leaves = 0;
   num_leaves_before_run = 0;
+  total_repair_lp = 0;
+  total_repair_lp_feasible = 0;
+  total_repair_lp_iterations = 0;
   total_lp_iterations = 0;
   heuristic_lp_iterations = 0;
   sepa_lp_iterations = 0;
@@ -801,16 +804,31 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
         fixedModel.col_upper_[i] = std::min(fixedModel.col_upper_[i], solval);
       }
     }
+    this->total_repair_lp++;
+    double time_available =
+        std::max(mipsolver.options_mip_->time_limit -
+                     mipsolver.timer_.read(mipsolver.timer_.solve_clock),
+                 0.1);
     Highs tmpSolver;
-    tmpSolver.setOptionValue("output_flag", false);
-    tmpSolver.setOptionValue("simplex_scale_strategy", 0);
-    tmpSolver.setOptionValue("presolve", "off");
+    const bool debug_report = false;
+    if (debug_report) {
+      tmpSolver.setOptionValue("log_dev_level", 2);
+      tmpSolver.setOptionValue("highs_analysis_level", 4);
+    } else {
+      tmpSolver.setOptionValue("output_flag", false);
+    }
+    // tmpSolver.setOptionValue("simplex_scale_strategy", 0);
+    // tmpSolver.setOptionValue("presolve", "off");
+    tmpSolver.setOptionValue("time_limit", time_available);
     tmpSolver.setOptionValue("primal_feasibility_tolerance",
                              mipsolver.options_mip_->mip_feasibility_tolerance);
     tmpSolver.passModel(std::move(fixedModel));
     tmpSolver.run();
 
+    this->total_repair_lp_iterations =
+        tmpSolver.getInfo().simplex_iteration_count;
     if (tmpSolver.getInfo().primal_solution_status == kSolutionStatusFeasible) {
+      this->total_repair_lp_feasible++;
       solution = tmpSolver.getSolution();
       allow_try_again = false;
       goto try_again;

src/mip/HighsMipSolverData.h

@@ -94,6 +94,9 @@ struct HighsMipSolverData {
   int64_t num_leaves;
   int64_t num_leaves_before_run;
   int64_t num_nodes_before_run;
+  int64_t total_repair_lp;
+  int64_t total_repair_lp_feasible;
+  int64_t total_repair_lp_iterations;
   int64_t total_lp_iterations;
   int64_t heuristic_lp_iterations;
   int64_t sepa_lp_iterations;
@@ -154,6 +157,9 @@ struct HighsMipSolverData {
         num_leaves(0),
         num_leaves_before_run(0),
         num_nodes_before_run(0),
+        total_repair_lp(0),
+        total_repair_lp_feasible(0),
+        total_repair_lp_iterations(0),
         total_lp_iterations(0),
         heuristic_lp_iterations(0),
         sepa_lp_iterations(0),

src/mip/HighsPrimalHeuristics.cpp

@@ -37,6 +37,9 @@
 HighsPrimalHeuristics::HighsPrimalHeuristics(HighsMipSolver& mipsolver)
     : mipsolver(mipsolver),
       lp_iterations(0),
+      total_repair_lp(0),
+      total_repair_lp_feasible(0),
+      total_repair_lp_iterations(0),
       randgen(mipsolver.options_mip_->random_seed) {
   successObservations = 0;
   numSuccessObservations = 0;
@@ -150,7 +153,10 @@ bool HighsPrimalHeuristics::solveSubMip(
     int64_t adjusted_lp_iterations =
         (size_t)(adjustmentfactor * submipsolver.mipdata_->total_lp_iterations);
     lp_iterations += adjusted_lp_iterations;
-
+    total_repair_lp += submipsolver.mipdata_->total_repair_lp;
+    total_repair_lp_feasible += submipsolver.mipdata_->total_repair_lp_feasible;
+    total_repair_lp_iterations +=
+        submipsolver.mipdata_->total_repair_lp_iterations;
     if (mipsolver.submip)
       mipsolver.mipdata_->num_nodes += std::max(
           int64_t{1}, int64_t(adjustmentfactor * submipsolver.node_count_));
@@ -1199,6 +1205,12 @@ void HighsPrimalHeuristics::clique() {
 #endif
 
 void HighsPrimalHeuristics::flushStatistics() {
+  mipsolver.mipdata_->total_repair_lp += total_repair_lp;
+  mipsolver.mipdata_->total_repair_lp_feasible += total_repair_lp_feasible;
+  mipsolver.mipdata_->total_repair_lp_iterations += total_repair_lp_iterations;
+  total_repair_lp = 0;
+  total_repair_lp_feasible = 0;
+  total_repair_lp_iterations = 0;
   mipsolver.mipdata_->heuristic_lp_iterations += lp_iterations;
   mipsolver.mipdata_->total_lp_iterations += lp_iterations;
   lp_iterations = 0;

src/mip/HighsPrimalHeuristics.h

@@ -22,6 +22,9 @@ class HighsMipSolver;
 class HighsPrimalHeuristics {
  private:
   HighsMipSolver& mipsolver;
+  size_t total_repair_lp;
+  size_t total_repair_lp_feasible;
+  size_t total_repair_lp_iterations;
   size_t lp_iterations;
 
   double successObservations;