IIS "light"

FEATURES.md

@@ -8,4 +8,6 @@ Now handling correctly the case where an infeasible MIP has a feasible relaxatio
 
 Fixed minor bug exposed by [#2441](https://github.com/ERGO-Code/HiGHS/issues/2441) in Highs::setSolution() for a sparse user solution when the moidel is empty, and only clearing the dual data before solving with modified objective in Highs::multiobjectiveSolve() so that user-supplied solution is not cleared.
 
+The irreducible infeasibility system (IIS) facility now detects infeasibility due to bounds on constraint activity values being incompatible with constraint bounds. A kIisStrategyLight mode for the iis_strategy option has been introduced so that only infeasibility due to incompatible variable/constraint bounds and constraint activity values is checked for. The LP corresponding to any known IIS is now formed and held as a data member of the HighsIis class. 
+
 

check/TestIis.cpp

@@ -5,7 +5,7 @@
 #include "Highs.h"
 #include "catch.hpp"
 
-const bool dev_run = false;
+const bool dev_run = false;  // true;  //
 const double inf = kHighsInf;
 
 void testIis(const std::string& model, const HighsIis& iis);
@@ -19,6 +19,8 @@
     const double rhs_penalty,
     const double require_feasibility_objective_function_value);
 
+void checkIisLp(HighsLp& lp, const HighsIis& iis);
+
 TEST_CASE("lp-incompatible-bounds", "[iis]") {
   // LP has row0 and col2 with inconsistent bounds.
   //
@@ -42,10 +44,21 @@
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
   highs.passModel(lp);
+  // Perform the light IIS check
+  highs.setOptionValue("iis_strategy", kIisStrategyLight);
+  HighsIis iis;
+  REQUIRE(highs.getIis(iis) == HighsStatus::kOk);
+  checkIisLp(lp, iis);
+
+  highs.passModel(lp);
+  highs.getIis(iis);
+  REQUIRE(highs.checkIis() == HighsStatus::kOk);
+
+  // Perform full IIS
   REQUIRE(highs.run() == HighsStatus::kOk);
   REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible);
   highs.setOptionValue("iis_strategy", kIisStrategyFromLpRowPriority);
-  HighsIis iis;
+
   REQUIRE(highs.getIis(iis) == HighsStatus::kOk);
   REQUIRE(iis.col_index_.size() == 0);
   REQUIRE(iis.row_index_.size() == 1);
@@ -87,7 +100,13 @@
   REQUIRE(iis.row_index_.size() == 1);
   REQUIRE(iis.row_index_[0] == empty_row);
   REQUIRE(iis.row_bound_[0] == kIisBoundStatusLower);
-  REQUIRE(highs.changeRowBounds(empty_row, -2, -1) == HighsStatus::kOk);
+  double new_lower = -2;
+  double new_upper = -1;
+  assert(new_upper < 0);
+  REQUIRE(highs.changeRowBounds(empty_row, new_lower, new_upper) ==
+          HighsStatus::kOk);
+  lp.row_lower_[empty_row] = new_lower;
+  lp.row_upper_[empty_row] = new_upper;
   REQUIRE(highs.run() == HighsStatus::kOk);
   REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible);
   REQUIRE(highs.getIis(iis) == HighsStatus::kOk);
@@ -96,11 +115,80 @@
   REQUIRE(iis.row_index_[0] == empty_row);
   REQUIRE(iis.row_bound_[0] == kIisBoundStatusUpper);
 
+  checkIisLp(lp, iis);
+
+  highs.passModel(lp);
+  highs.getIis(iis);
+  REQUIRE(highs.checkIis() == HighsStatus::kOk);
+
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("lp-get-iis-light", "[iis]") {
+  HighsLp lp;
+  lp.model_name_ = "lp-get-iis-light";
+  lp.num_col_ = 4;
+  lp.num_row_ = 3;
+  lp.col_cost_ = {0, 0, 0, 0};
+  lp.col_lower_ = {10, 10, 10, 10};
+  lp.col_upper_ = {20, 20, 20, 20};
+  lp.row_lower_ = {-inf, -10, -34};
+  lp.row_upper_ = {30, 15, inf};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.num_col_ = lp.num_col_;
+  lp.a_matrix_.num_row_ = lp.num_row_;
+  lp.a_matrix_.start_ = {0, 3, 7, 10};
+  lp.a_matrix_.index_ = {0, 1, 2, 0, 1, 2, 3, 1, 2, 3};
+  lp.a_matrix_.value_ = {1.5, 2, 1, 4, -2, 1, 2, -2, -1.5, -1};
+  //
+  // 1.5w + 2x + y <= 30
+  //
+  // -10 <= 4w -2x + y + 2z <= 15
+  //
+  // -2x -1.5y -z >= -34
+  //
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.passModel(lp);
+  highs.setOptionValue("iis_strategy", kIisStrategyLight);
+  HighsIis iis;
+
+  for (int l = 0; l < 3; l++) {
+    for (int k = 0; k < 2; k++) {
+      REQUIRE(highs.getIis(iis) == HighsStatus::kOk);
+      REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible);
+      const bool write_model = false;
+      if (dev_run && write_model) {
+        highs.writeModel("");
+        printf("\nNow pass IIS LP to write it out\n");
+        highs.passModel(iis.lp_);
+        highs.writeModel("");
+      }
+      checkIisLp(lp, iis);
+
+      highs.passModel(lp);
+      highs.getIis(iis);
+      REQUIRE(highs.checkIis() == HighsStatus::kOk);
+
+      if (k == 0) {
+        // Now flip to column-wise for code coverage
+        lp.a_matrix_.ensureColwise();
+        highs.passModel(lp);
+      }
+    }
+    if (l == 0) {
+      lp.row_upper_[0] = inf;
+    } else if (l == 1) {
+      lp.row_upper_[1] = inf;
+    }
+    highs.passModel(lp);
+  }
   highs.resetGlobalScheduler(true);
 }
 
 TEST_CASE("lp-get-iis", "[iis]") {
   HighsLp lp;
+  lp.model_name_ = "lp-get-iis";
   lp.num_col_ = 2;
   lp.num_row_ = 3;
   lp.col_cost_ = {0, 0};
@@ -109,9 +197,20 @@
   lp.row_lower_ = {-inf, -inf, -inf};
   lp.row_upper_ = {8, 9, -2};
   lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.num_col_ = lp.num_col_;
+  lp.a_matrix_.num_row_ = lp.num_row_;
   lp.a_matrix_.start_ = {0, 2, 4, 6};
   lp.a_matrix_.index_ = {0, 1, 0, 1, 0, 1};
   lp.a_matrix_.value_ = {2, 1, 1, 3, 1, 1};
+  //
+  // 2x + y <= 8
+  //
+  // x + 3y <= 9
+  //
+  // x +  y <= -2
+  //
+  // x, y \in [0, inf)
+  //
   //  lp.col_name_ = {"Col0", "Col1"};
   //  lp.row_name_ = {"Row0", "Row1", "Row2"};
   Highs highs;
@@ -128,6 +227,12 @@
   REQUIRE(iis.col_index_[1] == 1);
   REQUIRE(iis.row_index_[0] == 2);
 
+  checkIisLp(lp, iis);
+
+  highs.passModel(lp);
+  highs.getIis(iis);
+  REQUIRE(highs.checkIis() == HighsStatus::kOk);
+
   highs.resetGlobalScheduler(true);
 }
 
@@ -197,6 +302,7 @@
 TEST_CASE("lp-feasibility-relaxation", "[iis]") {
   // Using infeasible LP from AMPL documentation
   HighsLp lp;
+  lp.model_name_ = "ampl_infeas";
   lp.num_col_ = 2;
   lp.num_row_ = 3;
   lp.col_cost_ = {1, -2};
@@ -453,6 +559,23 @@
       printf("Model %s has IIS with %d columns and %d rows\n", model.c_str(),
              int(num_iis_col), int(num_iis_row));
     testIis(model, iis);
+
+    const bool write_model = false;
+    if (dev_run && write_model) highs.writeModel("");
+    HighsLp lp = highs.getLp();
+
+    if (dev_run && write_model) {
+      printf("\nNow pass IIS LP to write it out\n");
+      highs.passModel(iis.lp_);
+      highs.writeModel("");
+    }
+
+    checkIisLp(lp, iis);
+
+    highs.passModel(lp);
+    highs.getIis(iis);
+    REQUIRE(highs.checkIis() == HighsStatus::kOk);
+
   } else {
     REQUIRE(num_iis_col == 0);
     REQUIRE(num_iis_row == 0);
@@ -473,3 +596,127 @@
   REQUIRE(h.getInfo().objective_function_value ==
           require_feasibility_objective_function_value);
 }
+
+void checkIisLp(HighsLp& lp, const HighsIis& iis) {
+  const HighsLp& iis_lp = iis.lp_;
+  HighsInt iis_num_col = iis.col_index_.size();
+  HighsInt iis_num_row = iis.row_index_.size();
+  REQUIRE(iis_lp.num_col_ == iis_num_col);
+  REQUIRE(iis_lp.num_row_ == iis_num_row);
+
+  lp.a_matrix_.num_col_ = lp.num_col_;
+  lp.a_matrix_.num_row_ = lp.num_row_;
+  lp.a_matrix_.ensureColwise();
+  std::vector<HighsInt> iis_row;
+  iis_row.assign(lp.num_row_, -1);
+  double bound;
+  for (HighsInt iisRow = 0; iisRow < iis_num_row; iisRow++) {
+    HighsInt iRow = iis.row_index_[iisRow];
+    if (iRow < 0 || iRow >= lp.num_row_) {
+      printf("iRow out of range\n");
+    }
+    iis_row[iRow] = iisRow;
+    HighsInt row_bound = iis.row_bound_[iisRow];
+    bound =
+        row_bound == kIisBoundStatusLower || row_bound == kIisBoundStatusBoxed
+            ? lp.row_lower_[iRow]
+            : -kHighsInf;
+    const bool lower_ok = iis_lp.row_lower_[iisRow] == bound;
+    if (!lower_ok)
+      printf(
+          "!lower_ok iis_lp.row_lower_[iisRow] = %g; bound = %g; "
+          "iis_lp.row_lower_[iisRow] - bound = %g\n",
+          iis_lp.row_lower_[iisRow], bound, iis_lp.row_lower_[iisRow] - bound);
+    REQUIRE(iis_lp.row_lower_[iisRow] == bound);
+    bound =
+        row_bound == kIisBoundStatusUpper || row_bound == kIisBoundStatusBoxed
+            ? lp.row_upper_[iRow]
+            : kHighsInf;
+    REQUIRE(iis_lp.row_upper_[iisRow] == bound);
+  }
+
+  // Work through the LP columns and matrix, checking the costs,
+  // bounds and matrix index/value
+  const HighsInt illegal_index = -1;
+  const double illegal_value = kHighsInf;
+  std::vector<HighsInt> index;
+  std::vector<double> value;
+  for (HighsInt iisCol = 0; iisCol < iis_num_col; iisCol++) {
+    HighsInt iCol = iis.col_index_[iisCol];
+    REQUIRE(iis_lp.col_cost_[iisCol] == lp.col_cost_[iCol]);
+    HighsInt col_bound = iis.col_bound_[iisCol];
+    bound =
+        col_bound == kIisBoundStatusLower || col_bound == kIisBoundStatusBoxed
+            ? lp.col_lower_[iCol]
+            : -kHighsInf;
+    REQUIRE(iis_lp.col_lower_[iisCol] == bound);
+    bound =
+        col_bound == kIisBoundStatusUpper || col_bound == kIisBoundStatusBoxed
+            ? lp.col_upper_[iCol]
+            : kHighsInf;
+    REQUIRE(iis_lp.col_upper_[iisCol] == bound);
+    // Use index/value to scatter the IIS matrix column
+    index.assign(iis_num_row, illegal_index);
+    value.assign(iis_num_row, illegal_value);
+    for (HighsInt iEl = iis_lp.a_matrix_.start_[iisCol];
+         iEl < iis_lp.a_matrix_.start_[iisCol + 1]; iEl++) {
+      HighsInt iisRow = iis_lp.a_matrix_.index_[iEl];
+      HighsInt iRow = iis.row_index_[iisRow];
+      index[iisRow] = iRow;
+      value[iisRow] = iis_lp.a_matrix_.value_[iEl];
+    }
+    for (HighsInt iEl = lp.a_matrix_.start_[iCol];
+         iEl < lp.a_matrix_.start_[iCol + 1]; iEl++) {
+      HighsInt iRow = lp.a_matrix_.index_[iEl];
+      HighsInt iisRow = iis_row[iRow];
+      if (iisRow >= 0) {
+        const bool index_ok = index[iisRow] == iRow;
+        if (!index_ok && dev_run) {
+          printf("checkIisLp: IIS LP matrix index incorrect %d != %d\n",
+                 int(index[iisRow]), int(iRow));
+        }
+        REQUIRE(index_ok);
+        index[iisRow] = illegal_index;
+        const bool value_ok = value[iisRow] == lp.a_matrix_.value_[iEl];
+        if (!value_ok && dev_run) {
+          printf("checkIisLp: IIS LP matrix value incorrect %g != %g\n",
+                 value[iisRow], lp.a_matrix_.value_[iEl]);
+        }
+        REQUIRE(value_ok);
+        value[iisRow] = illegal_value;
+      }
+    }
+  }
+  // Work through the IIS LP matrix, making sure that the index/value
+  // are correct
+  for (HighsInt iisCol = 0; iisCol < iis_num_col; iisCol++) {
+    HighsInt iCol = iis.col_index_[iisCol];
+    // Use index/value to scatter the LP matrix column
+    index.assign(lp.num_row_, illegal_index);
+    value.assign(lp.num_row_, illegal_value);
+    for (HighsInt iEl = lp.a_matrix_.start_[iCol];
+         iEl < lp.a_matrix_.start_[iCol + 1]; iEl++) {
+      HighsInt iRow = lp.a_matrix_.index_[iEl];
+      HighsInt iisRow = iis_row[iRow];
+      index[iRow] = iisRow;
+      value[iRow] = lp.a_matrix_.value_[iEl];
+    }
+    for (HighsInt iEl = iis_lp.a_matrix_.start_[iisCol];
+         iEl < iis_lp.a_matrix_.start_[iisCol + 1]; iEl++) {
+      HighsInt iisRow = iis_lp.a_matrix_.index_[iEl];
+      HighsInt iRow = iis.row_index_[iisRow];
+      const bool index_ok = index[iRow] == iisRow;
+      if (!index_ok && dev_run) {
+        printf("checkIisLp: IIS LP matrix index incorrect %d != %d\n",
+               int(index[iRow]), int(iisRow));
+      }
+      REQUIRE(index_ok);
+      const bool value_ok = value[iRow] == iis_lp.a_matrix_.value_[iEl];
+      if (!value_ok && dev_run) {
+        printf("checkIisLp: IIS LP matrix value incorrect %g != %g\n",
+               value[iRow], iis_lp.a_matrix_.value_[iEl]);
+      }
+      REQUIRE(value_ok);
+    }
+  }
+}

highs/Highs.h

@@ -573,6 +573,11 @@ class Highs {
    */
   HighsStatus getIis(HighsIis& iis);
 
+  /**
+   * @brief Checks the IIS information for the incumbent model
+   */
+  HighsStatus checkIis() const;
+
   /**
    * @brief Get the current model objective function value
    */
@@ -1660,7 +1665,8 @@ class Highs {
 
   HighsStatus elasticityFilterReturn(
       const HighsStatus return_status, const bool feasible_model,
-      const HighsInt original_num_col, const HighsInt original_num_row,
+      const std::string& original_model_name, const HighsInt original_num_col,
+      const HighsInt original_num_row,
       const std::vector<double>& original_col_cost,
       const std::vector<double>& original_col_lower,
       const std::vector<double> original_col_upper,

highs/lp_data/HConst.h

@@ -270,10 +270,11 @@ enum PresolveRuleType : int {
 
 enum IisStrategy {
   kIisStrategyMin = 0,
-  kIisStrategyFromLpRowPriority = kIisStrategyMin,  // 0
-  kIisStrategyFromLpColPriority,                    // 1
-  //  kIisStrategyFromRayRowPriority,                     // 2
-  //  kIisStrategyFromRayColPriority,                     // 3
+  kIisStrategyLight = kIisStrategyMin,  // 0
+  kIisStrategyFromLpRowPriority,        // 1
+  kIisStrategyFromLpColPriority,        // 2
+  //  kIisStrategyFromRayRowPriority,                     // 3
+  //  kIisStrategyFromRayColPriority,                     // 4
   kIisStrategyMax = kIisStrategyFromLpColPriority
 };