Merge pull request #3248 from ERGO-Code/fix-dual-flip-after-postsolve

Fix dual flip after postsolve

Author: jajhall

check/TestPresolve.cpp

@@ -1154,6 +1154,65 @@ TEST_CASE("add-to-matrix", "[highs_test_presolve]") {
   h.resetGlobalScheduler(true);
 }
 
+TEST_CASE("max-lp-dual-postsolve", "[highs_test_presolve]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+
+  // Start with blending problem
+  HighsLp lp;
+  lp.model_name_ = "blending";
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {8, 10};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {inf, inf};
+  lp.row_lower_ = {-inf, -inf};
+  lp.row_upper_ = {120, 210};
+  lp.a_matrix_.format_ = MatrixFormat::kColwise;
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {0, 1, 0, 1};
+  lp.a_matrix_.value_ = {0.3, 0.7, 0.5, 0.5};
+  lp.sense_ = ObjSense::kMaximize;
+  lp.offset_ = 0;
+  h.passModel(lp);
+
+  REQUIRE(h.setOptionValue("solver", kIpxString) == HighsStatus::kOk);
+  REQUIRE(h.setOptionValue("run_crossover", kHighsOffString) ==
+          HighsStatus::kOk);
+
+  h.run();
+  if (dev_run) h.writeSolution("", kSolutionStylePretty);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+
+  double optimal_objective = h.getObjectiveValue();
+  HighsSolution solution = h.getSolution();
+
+  // Add a redundant constraint that's removed in presolve so that
+  // postsolve must be performed
+  std::vector<HighsInt> index = {0, 1};
+  std::vector<double> value = {1, 1};
+  REQUIRE(h.addRow(-inf, 1000, 2, index.data(), value.data()) ==
+          HighsStatus::kOk);
+  h.clearSolver();
+
+  h.run();
+  if (dev_run) h.writeSolution("", kSolutionStylePretty);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+
+  REQUIRE(doubleEqual(optimal_objective, h.getObjectiveValue()));
+  // Check that the dual values are unchanged for the original problem
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+    REQUIRE(
+        doubleEqual(solution.col_dual[iCol], h.getSolution().col_dual[iCol]));
+  for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++)
+    REQUIRE(
+        doubleEqual(solution.row_dual[iRow], h.getSolution().row_dual[iRow]));
+
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+
+  h.resetGlobalScheduler(true);
+}
+
 TEST_CASE("issue-3140", "[highs_test_presolve]") {
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);

highs/ipm/hipo/ipm/Solver.cpp

@@ -1213,7 +1213,7 @@ isSuccess Solver::checkTerminationKkt() {
                       highs_solution, Hoptions_, "During HiPO solve");
 
     if (model_status == HighsModelStatus::kOptimal) {
-      logger_.printInfo("Check successfull\n");
+      logger_.printInfo("Check successful\n");
       return true;
     } else
       logger_.printInfo("Check failed\n");
@@ -1282,7 +1282,7 @@ isFailure Solver::initialiseLinearSolver() {
 }
 
 isSuccess Solver::switchToMultifrontal() {
-  bool switch_successfull = false;
+  bool switch_successful = false;
 
   if (LS_->type() == kUpLookingType && options_.factor == kHighsChooseString) {
     LS_.reset(new FactorHighsSolver(*kkt_, options_, model_, regul_, info_,
@@ -1297,11 +1297,11 @@ isSuccess Solver::switchToMultifrontal() {
       it_->largest_dx_x_ = 0;
       it_->largest_dy_y_ = 0;
       info_.error = kOk;
-      switch_successfull = true;
+      switch_successful = true;
     }
   }
 
-  return switch_successfull;
+  return switch_successful;
 }
 
 void Solver::printHeader() const {

highs/lp_data/Highs.cpp

@@ -4009,8 +4009,13 @@ HighsPostsolveStatus Highs::runPostsolve() {
   assert(model_.lp_.a_matrix_.isColwise());
   calculateRowValuesQuad(model_.lp_, presolve_.data_.recovered_solution_);
 
-  if (have_dual_solution && model_.lp_.sense_ == ObjSense::kMaximize)
-    presolve_.negateReducedLpColDuals();
+  if (have_dual_solution && model_.lp_.sense_ == ObjSense::kMaximize) {
+    // Negate the dual values since the incumbent LP is a maximization
+    for (HighsInt iCol = 0; iCol < model_.lp_.num_col_; iCol++)
+      presolve_.data_.recovered_solution_.col_dual[iCol] *= -1;
+    for (HighsInt iRow = 0; iRow < model_.lp_.num_row_; iRow++)
+      presolve_.data_.recovered_solution_.row_dual[iRow] *= -1;
+  }
 
   // Ensure that the postsolve status is used to set
   // presolve_.postsolve_status_, as well as being returned

highs/presolve/PresolveComponent.cpp

@@ -21,13 +21,6 @@ HighsStatus PresolveComponent::init(const HighsLp& lp, HighsTimer& timer,
   return HighsStatus::kOk;
 }
 
-void PresolveComponent::negateReducedLpColDuals() {
-  for (HighsInt col = 0; col < data_.reduced_lp_.num_col_; col++)
-    data_.recovered_solution_.col_dual[col] =
-        -data_.recovered_solution_.col_dual[col];
-  return;
-}
-
 HighsPresolveStatus PresolveComponent::run() {
   presolve::HPresolve presolve;
   presolve.setInput(data_.reduced_lp_, *options_,

highs/presolve/PresolveComponent.h

@@ -75,8 +75,6 @@ class PresolveComponent : public HighsComponent {
   HighsLp& getReducedProblem() { return data_.reduced_lp_; }
   HighsPresolveLog& getPresolveLog() { return data_.presolve_log_; }
 
-  void negateReducedLpColDuals();
-
   PresolveComponentInfo info_;
   PresolveComponentData data_;
   const HighsOptions* options_;