Merge pull request #1984 from ERGO-Code/fix-1956-1977-1981

Forces calculation of primal/dual ray and computed dual objective
ERGO-Code · Oct 26, 2024 · ace2579 · ace2579
2 parents 537423d + 0f478a9
commit ace2579
Show file tree

Hide file tree

Showing 19 changed files with 1,120 additions and 437 deletions.
diff --git a/check/CMakeLists.txt b/check/CMakeLists.txt
@@ -267,7 +267,7 @@ if (NOT FAST_BUILD OR ALL_TESTS)
                   DEPENDS unit_tests_all)
           set_tests_properties (${name}${setting} PROPERTIES
                   PASS_REGULAR_EXPRESSION
-                  "Model   status      : ${solutionstatus}")
+                  "Model status        : ${solutionstatus}")
 
           if(${solutionstatus} STREQUAL "Optimal")
               if(${setting} STREQUAL "--presolve=off")

diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp
@@ -2,14 +2,33 @@
 #include "Highs.h"
 #include "catch.hpp"
 
-const bool dev_run = true;  // false;
+const bool dev_run = false;
 
 struct IterationCount {
   HighsInt simplex;
   HighsInt ipm;
   HighsInt crossover;
 };
 
+void testDualObjective(const std::string model) {
+  HighsStatus return_status;
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  highs.readModel(model_file);
+  return_status = highs.run();
+  REQUIRE(return_status == HighsStatus::kOk);
+  double dual_objective;
+  return_status = highs.getDualObjectiveValue(dual_objective);
+  REQUIRE(return_status == HighsStatus::kOk);
+  double primal_objective = highs.getInfo().objective_function_value;
+  double relative_primal_dual_gap =
+      std::fabs(primal_objective - dual_objective) /
+      std::max(1.0, std::fabs(primal_objective));
+  REQUIRE(relative_primal_dual_gap < 1e-12);
+}
 void testSolver(Highs& highs, const std::string solver,
                 IterationCount& default_iteration_count,
                 const HighsInt int_simplex_strategy = 0) {
@@ -469,6 +488,42 @@ TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") {
     printf("Max   dual infeasibility   = %g\n", info.max_dual_infeasibility);
     printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
   }
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+}
+
+TEST_CASE("dual-objective-max", "[highs_lp_solver]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.sense_ = ObjSense::kMaximize;
+  lp.offset_ = 10;
+  lp.col_cost_ = {8, 10};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf};
+  lp.row_lower_ = {-kHighsInf, -kHighsInf};
+  lp.row_upper_ = {80, 120};
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {0, 1, 0, 1};
+  lp.a_matrix_.value_ = {1, 1, 2, 4};
+  highs.passModel(lp);
+  highs.run();
+  double dual_objective;
+  HighsStatus return_status = highs.getDualObjectiveValue(dual_objective);
+  REQUIRE(return_status == HighsStatus::kOk);
+  double primal_objective = highs.getInfo().objective_function_value;
+  double relative_primal_dual_gap =
+      std::fabs(primal_objective - dual_objective) /
+      std::max(1.0, std::fabs(primal_objective));
+  REQUIRE(relative_primal_dual_gap < 1e-12);
+}
+
+TEST_CASE("dual-objective", "[highs_lp_solver]") {
+  testDualObjective("avgas");
+  testDualObjective("adlittle");
+  testDualObjective("etamacro");
+  testDualObjective("stair");
 }
 
 void testStandardForm(const HighsLp& lp) {

diff --git a/check/TestMipSolver.cpp b/check/TestMipSolver.cpp
@@ -653,6 +653,37 @@ TEST_CASE("MIP-get-saved-solutions-presolve", "[highs_test_mip_solver]") {
   std::remove(solution_file.c_str());
 }
 
+TEST_CASE("IP-infeasible-unbounded", "[highs_test_mip_solver]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  double delta = 0.2;
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 0;
+  lp.col_cost_ = {-1, 0};
+  lp.integrality_ = {HighsVarType::kInteger, HighsVarType::kInteger};
+  highs.setOptionValue("presolve", kHighsOffString);
+  for (HighsInt k = 0; k < 2; k++) {
+    for (HighsInt l = 0; l < 2; l++) {
+      if (l == 0) {
+        // Infeasible
+        lp.col_lower_ = {0, delta};
+        lp.col_upper_ = {kHighsInf, 1 - delta};
+      } else {
+        // Unbounded
+        lp.col_lower_ = {0, -delta};
+        lp.col_upper_ = {kHighsInf, 1 + delta};
+      }
+      // Solve
+      highs.passModel(lp);
+      highs.run();
+      REQUIRE(highs.getModelStatus() ==
+              HighsModelStatus::kUnboundedOrInfeasible);
+    }
+    highs.setOptionValue("presolve", kHighsOnString);
+  }
+}
+
 TEST_CASE("IP-with-fract-bounds-no-presolve", "[highs_test_mip_solver]") {
   Highs highs;
   // No presolve