diff --git a/check/TestIis.cpp b/check/TestIis.cpp index 456d99f177..70dcc73520 100644 --- a/check/TestIis.cpp +++ b/check/TestIis.cpp @@ -207,14 +207,23 @@ TEST_CASE("lp-get-iis", "[iis]") { highs.setOptionValue("output_flag", dev_run); highs.passModel(lp); HighsIis iis; - REQUIRE(highs.getIis(iis) == HighsStatus::kOk); - REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible); - REQUIRE(iis.col_index_.size() == 2); - REQUIRE(iis.row_index_.size() == 1); - REQUIRE(iis.col_index_[0] == 0); - REQUIRE(iis.col_index_[1] == 1); - REQUIRE(iis.row_index_[0] == 2); - + const HighsLp& highs_lp = highs.getLp(); + // First pass with incumbent matrix colwise; second with it + // rowwise + highs.ensureColwise(); + REQUIRE(highs_lp.a_matrix_.isColwise()); + for (HighsInt k = 0; k < 2; k++) { + REQUIRE(highs.getIis(iis) == HighsStatus::kOk); + REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible); + REQUIRE(iis.col_index_.size() == 2); + REQUIRE(iis.row_index_.size() == 1); + REQUIRE(iis.col_index_[0] == 0); + REQUIRE(iis.col_index_[1] == 1); + REQUIRE(iis.row_index_[0] == 2); + highs.clearSolver(); + highs.ensureRowwise(); + REQUIRE(highs_lp.a_matrix_.isRowwise()); + } highs.resetGlobalScheduler(true); } @@ -397,7 +406,7 @@ void testMps(std::string& model, const HighsInt iis_strategy, std::string model_file = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps"; Highs highs; - highs.setOptionValue("output_flag", dev_run); + // highs.setOptionValue("output_flag", dev_run); REQUIRE(highs.readModel(model_file) == HighsStatus::kOk); // if (iis_strategy == kIisStrategyFromRayRowPriority || diff --git a/docs/make.jl b/docs/make.jl index 4287467b79..724620e8f4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,9 +2,6 @@ * * * This file is part of the HiGHS linear optimization suite * * * -* Written and engineered 2008-2024 by Julian Hall, Ivet Galabova, * -* Leona Gottwald and Michael Feldmeier * -* * * Available as open-source under the MIT License * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *=# @@ -68,7 +65,8 @@ Documenter.makedocs( "structures/classes/HighsSparseMatrix.md", "structures/classes/HighsLp.md", "structures/classes/HighsHessian.md", - "structures/classes/HighsModel.md" + "structures/classes/HighsModel.md", + "structures/classes/HighsIis.md" ], "Structures" => Any[ "structures/structs/index.md", diff --git a/docs/src/guide/advanced.md b/docs/src/guide/advanced.md index 8a39f41241..037344a4d4 100644 --- a/docs/src/guide/advanced.md +++ b/docs/src/guide/advanced.md @@ -29,3 +29,34 @@ and yield a specific row or column of ``B^{-1}A``. In all cases, HiGHS can return the number and indices of the nonzeros in the result. +## [Irreducible infeasibility system (IIS) detection](@id highs-iis) + +An Irreducible infeasibility system (IIS) consists of a set of +variables and a set of constraints in a model, together with +variable/constraint bound information, that cannot be satisfied (so is +infeasible). It is irreducible in that if any constraint or variable +bound is removed, then the system can be satisfied (so is feasible). + +HiGHS has an IIS facility that is under development. Currently it can +only be used for LPs. The full IIS calculation is expensive, since it +requires the solution of multiple LPs. Although there is a prototype +implementation, it is not as robust or efficient as it will +be. Otherwise, there is a simple, cheap test that looks for +infeasibility due to incompatible variable oe constraint bounds, or +constraint bounds that cannot be satisfied given the range of values +on the constraint activity implied by bounds on variables. + +The choice of IIS strategy is defined by the [iis_strategy](@id option-iis-strategy) option, which can take the value + +- `kIisStrategyLight` = 0: The cheap test +- `kIisStrategyFromLpRowPriority` = 1: The full IIS calculation, aiming to have a minimal number of rows in the IIS +- `kIisStrategyFromLpColPriority` = 2: The full IIS calculation, aiming to have a minimal number of columns in the IIS + +### IIS-related methods in the `Highs` class + +- `const HighsLp& getIisLp()`: Return a const reference to the internal IIS LP instance +- `HighsStatus getIis(HighsIis& iis)`: Try to find an IIS for the incumbent model. Gets the internal [HighsIis](@ref) instance, returning `HighsStatus::kError` if the calculation failed. Note that if the incumbent model is found to be feasible, this is a "success", and `HighsStatus::kOk` is returned. +- ` HighsStatus writeIisModel(const std::string& filename = "")`: Write out the internal IIS LP instance to a file. + + + diff --git a/docs/src/guide/further.md b/docs/src/guide/further.md index ed111f6682..34b0de5851 100644 --- a/docs/src/guide/further.md +++ b/docs/src/guide/further.md @@ -125,8 +125,8 @@ linear objective is represented by the following data, held in the Multi-objective optimization in HiGHS is defined by the following methods -- [addLinearObjective](@ref Multi-objective-optimization] - Add a single `HighsLinearObjective` instance to any already stored in HiGHS -- [clearLinearObjectives](@ref Multi-objective-optimization] - Clears any linear objectives stored in HiGHS +- [addLinearObjective](@ref Multi-objective-optimization) - Add a single `HighsLinearObjective` instance to any already stored in HiGHS +- [clearLinearObjectives](@ref Multi-objective-optimization) - Clears any linear objectives stored in HiGHS When there is at least one `HighsLinearObjective` instance in HiGHS, the `col_cost_` data in the incumbent model is ignored. diff --git a/docs/src/options/definitions.md b/docs/src/options/definitions.md index 8bc8b039e0..df0e2cd8d9 100644 --- a/docs/src/options/definitions.md +++ b/docs/src/options/definitions.md @@ -486,7 +486,7 @@ - Range: [0, inf] - Default: 1e-07 -## iis\_strategy +## [iis\_strategy](@id option-iis-strategy) - Strategy for IIS calculation: Light test / Full and prioritise rows / Full and prioritise columns (0/1/2) - Type: integer - Range: {0, 2} diff --git a/docs/src/structures/classes/HighsHessian.md b/docs/src/structures/classes/HighsHessian.md index 5e165ff8ee..c1e1342947 100644 --- a/docs/src/structures/classes/HighsHessian.md +++ b/docs/src/structures/classes/HighsHessian.md @@ -2,8 +2,8 @@ A Hessian matrix is communicated via an instance of the HighsHessian class. -- dim_: Scalar of type integer - Dimension of the Hessian -- format\_: Scalar of [HessianFormat](@ref) type - Format of the Hessian -- start\_: Vector of integer type - Start of each compressed column in the Hessian -- index\_: Vector of integer type - Indices of the nonzeros in the Hessian -- value\_: Vector of double type - Values of the nonzeros in the Hessian +- `dim_`: Scalar of type integer - Dimension of the Hessian +- `format_`: Scalar of [HessianFormat](@ref) type - Format of the Hessian +- `start_`: Vector of integer type - Start of each compressed column in the Hessian +- `index_`: Vector of integer type - Indices of the nonzeros in the Hessian +- `value_`: Vector of double type - Values of the nonzeros in the Hessian diff --git a/docs/src/structures/classes/HighsIis.md b/docs/src/structures/classes/HighsIis.md new file mode 100644 index 0000000000..dbe6c4fa81 --- /dev/null +++ b/docs/src/structures/classes/HighsIis.md @@ -0,0 +1,16 @@ +# HighsIis + +Irreducible infeasibility system (IIS) data are communicated via an instance of the `HighsIis` class. + +- `valid_`: The data in the `HighsIis` instance is valid +- `strategy_`: The IIS strategy used +- `col_index_`: The indices of model columns in the IIS +- `row_index_`: The indices of model rows in the IIS +- `col_bound_`: The bounds on each column that define the IIS +- `row_bound_`: The bounds on each row that define the IIS +- `col_status_`: Indicates whether a column in the model is in an IIS, may be in an IIS, or is not in an IIS +- `row_status_`: Indicates whether a row in the model is in an IIS, may be in an IIS, or is not in an IIS +- `info_`: Data on the time and number of simplex iterations required to form the IIS +- `model_`: A [HighsModel](@ref) consisting of the variables, constraints and bounds in the IIS. Currently only its [HighsLp](@ref) instance is relevant + + diff --git a/docs/src/structures/classes/HighsLp.md b/docs/src/structures/classes/HighsLp.md index 44027869c7..0c73db7d3f 100644 --- a/docs/src/structures/classes/HighsLp.md +++ b/docs/src/structures/classes/HighsLp.md @@ -1,19 +1,19 @@ # HighsLp -An LP model is communicated via an instance of the HighsLp class +An LP model is communicated via an instance of the `HighsLp` class -- num\_col\_: Scalar of type integer - Number of columns in the model -- num\_row\_: Scalar of type integer - Number of rows in the model -- col\_cost\_: Vector of type double - Coefficients of the linear term in the objective function -- col\_lower\_: Vector of type double - Lower bounds on the variables -- col\_upper\_: Vector of type double - Upper bounds on the variables -- row\_lower\_: Vector of type double - Lower bounds on the constraints -- row\_upper\_: Vector of type double - Upper bounds on the constraints -- a\_matrix\_: Instance of [HighsSparseMatrix](@ref) class - Constraint matrix -- sense\_: Scalar of type [ObjSense](@ref) - Optimization sense of the model -- offset\_: Scalar of type double - Constant term in the objective function -- model\_name\_: Scalar of type string - Name of the model -- objective\_name\_: Scalar of type string - Name of the objective function -- col\_names\_: Vector of type string - Names of the variables -- row\_names\_: Vector of type string - Names of the constraints -- integrality\_: Vector of type [HighsVarType](@ref) - Type of each variable +- `num_col_`: Scalar of type integer - Number of columns in the model +- `num_row_`: Scalar of type integer - Number of rows in the model +- `col_cost_`: Vector of type double - Coefficients of the linear term in the objective function +- `col_lower_`: Vector of type double - Lower bounds on the variables +- `col_upper_`: Vector of type double - Upper bounds on the variables +- `row_lower_`: Vector of type double - Lower bounds on the constraints +- `row_upper_`: Vector of type double - Upper bounds on the constraints +- `a_matrix_`: Instance of [HighsSparseMatrix](@ref) class - Constraint matrix +- `sense_`: Scalar of type [ObjSense](@ref) - Optimization sense of the model +- `offset_`: Scalar of type double - Constant term in the objective function +- `model_name_`: Scalar of type string - Name of the model +- `objective_name_`: Scalar of type string - Name of the objective function +- `col_names_`: Vector of type string - Names of the variables +- `row_names_`: Vector of type string - Names of the constraints +- `integrality_`: Vector of type [HighsVarType](@ref) - Type of each variable diff --git a/docs/src/structures/classes/HighsModel.md b/docs/src/structures/classes/HighsModel.md index 21aa4c1e80..af6643304c 100644 --- a/docs/src/structures/classes/HighsModel.md +++ b/docs/src/structures/classes/HighsModel.md @@ -2,6 +2,5 @@ A QP model is communicated via an instance of the HighsModel class -- lp\_: Instance of [HighsLp](@ref) class - LP components of the model - -- hessian\_: Instance of [HighsHessian](@ref) class - Hessian matrix +- `lp_`: Instance of [HighsLp](@ref) class - LP components of the model +- `hessian_`: Instance of [HighsHessian](@ref) class - Hessian matrix diff --git a/docs/src/structures/classes/HighsSparseMatrix.md b/docs/src/structures/classes/HighsSparseMatrix.md index 411c9bcc1e..2838e50fc3 100644 --- a/docs/src/structures/classes/HighsSparseMatrix.md +++ b/docs/src/structures/classes/HighsSparseMatrix.md @@ -2,9 +2,9 @@ The constraint matrix of an LP model is communicated via an instance of the HighsSparseMatrix class -- format\_: Scalar of [MatrixFormat](@ref) type - Format of the matrix -- num\_col\_ : Scalar of integer type - Number of columns in the matrix -- num\_row\_: Scalar of integer type - Number of rows in the matrix -- start\_: Vector of integer type - Start of each compressed vector in the matrix -- index\_: Vector of integer type - Indices of the nonzeros in the matrix -- value\_: Vector of double type - Values of the nonzeros in the matrix +- `format_`: Scalar of [MatrixFormat](@ref) type - Format of the matrix +- `num_col_ `: Scalar of integer type - Number of columns in the matrix +- `num_row_`: Scalar of integer type - Number of rows in the matrix +- `start_`: Vector of integer type - Start of each compressed vector in the matrix +- `index_`: Vector of integer type - Indices of the nonzeros in the matrix +- `value_`: Vector of double type - Values of the nonzeros in the matrix diff --git a/docs/src/structures/classes/index.md b/docs/src/structures/classes/index.md index 9277bc5ef8..7f48a7030d 100644 --- a/docs/src/structures/classes/index.md +++ b/docs/src/structures/classes/index.md @@ -6,5 +6,6 @@ The data members of fundamental classes in HiGHS are defined in this section. * [HighsLp](@ref) * [HighsHessian](@ref) * [HighsModel](@ref) + * [HighsIis](@ref) Class data members for internal use only are not documented. diff --git a/highs/lp_data/Highs.cpp b/highs/lp_data/Highs.cpp index ace0ae07ae..efd6f64f14 100644 --- a/highs/lp_data/Highs.cpp +++ b/highs/lp_data/Highs.cpp @@ -1038,7 +1038,7 @@ HighsStatus Highs::optimizeModel() { // kHighsDebugLevelMax; // // if (model_.lp_.num_row_>0 && model_.lp_.num_col_>0) - // writeLpMatrixPicToFile(options_, "LpMatrix", model_.lp_); + // writeLpMatrixPicToFile(options_, "LpMatrix", model_.lp_); if (options_.highs_debug_level < min_highs_debug_level) options_.highs_debug_level = min_highs_debug_level; diff --git a/highs/lp_data/HighsIis.cpp b/highs/lp_data/HighsIis.cpp index 0bd8043606..4f7d731b15 100644 --- a/highs/lp_data/HighsIis.cpp +++ b/highs/lp_data/HighsIis.cpp @@ -124,11 +124,18 @@ bool HighsIis::trivial(const HighsLp& lp, const HighsOptions& options) { } // Now look for empty rows that cannot have zero activity std::vector count; - count.assign(lp.num_row_, 0); - for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { - for (HighsInt iEl = lp.a_matrix_.start_[iCol]; - iEl < lp.a_matrix_.start_[iCol + 1]; iEl++) - count[lp.a_matrix_.index_[iEl]]++; + // Get the row counts + if (lp.a_matrix_.isColwise()) { + count.assign(lp.num_row_, 0); + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + for (HighsInt iEl = lp.a_matrix_.start_[iCol]; + iEl < lp.a_matrix_.start_[iCol + 1]; iEl++) + count[lp.a_matrix_.index_[iEl]]++; + } + } else { + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) + count.push_back(lp.a_matrix_.start_[iRow + 1] - + lp.a_matrix_.start_[iRow]); } assert(this->row_index_.size() == 0); for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { @@ -289,6 +296,7 @@ HighsStatus HighsIis::getData(const HighsLp& lp, const HighsOptions& options, std::vector from_col; std::vector to_row; to_row.assign(lp.num_row_, -1); + // To get the IIS data needs the matrix to be column-wise assert(lp.a_matrix_.isColwise()); // Determine how to detect whether a row is in infeasible_row and // (then) gather information about it @@ -434,6 +442,8 @@ void HighsIis::getLp(const HighsLp& lp) { } iis_lp.num_col_ = iis_lp.col_cost_.size(); iis_lp.num_row_ = iis_lp.row_lower_.size(); + // The IIS LP matrix will have the same format as the incumbent LP + iis_lp.a_matrix_.format_ = lp.a_matrix_.format_; iis_lp.a_matrix_.num_col_ = iis_lp.num_col_; iis_lp.a_matrix_.num_row_ = iis_lp.num_row_; iis_lp.model_name_ = lp.model_name_ + "_IIS"; @@ -775,17 +785,25 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, return HighsStatus::kOk; } +bool lpDataOkReturn(const bool return_value) { return return_value; } + bool HighsIis::lpDataOk(const HighsLp& lp, const HighsOptions& options) const { const HighsLp& iis_lp = this->model_.lp_; HighsInt iis_num_col = this->col_index_.size(); HighsInt iis_num_row = this->row_index_.size(); - if (!(iis_lp.num_col_ == iis_num_col)) return false; - if (!(iis_lp.num_row_ == iis_num_row)) return false; + if (!(iis_lp.num_col_ == iis_num_col)) return lpDataOkReturn(false); + if (!(iis_lp.num_row_ == iis_num_row)) return lpDataOkReturn(false); const bool colwise = lp.a_matrix_.isColwise(); + const HighsInt illegal_index = -1; + const double illegal_value = kHighsInf; + // iis_row/col give the row/col in the IIS for each row/col in the + // LP, or an illegal index if the LP row/col isn't in the IIS std::vector iis_row; - iis_row.assign(lp.num_row_, -1); + iis_row.assign(lp.num_row_, illegal_index); + std::vector iis_col; + iis_col.assign(lp.num_col_, illegal_index); double bound; for (HighsInt iisRow = 0; iisRow < iis_num_row; iisRow++) { HighsInt iRow = this->row_index_[iisRow]; @@ -795,32 +813,31 @@ bool HighsIis::lpDataOk(const HighsLp& lp, const HighsOptions& options) const { row_bound == kIisBoundStatusLower || row_bound == kIisBoundStatusBoxed ? lp.row_lower_[iRow] : -kHighsInf; - if (iis_lp.row_lower_[iisRow] != bound) return false; + if (iis_lp.row_lower_[iisRow] != bound) return lpDataOkReturn(false); bound = row_bound == kIisBoundStatusUpper || row_bound == kIisBoundStatusBoxed ? lp.row_upper_[iRow] : kHighsInf; - if (iis_lp.row_upper_[iisRow] != bound) return false; + if (iis_lp.row_upper_[iisRow] != bound) return lpDataOkReturn(false); } // Work through the LP columns checking the zero costs and bounds for (HighsInt iisCol = 0; iisCol < iis_num_col; iisCol++) { HighsInt iCol = this->col_index_[iisCol]; - if (iis_lp.col_cost_[iisCol]) return false; + iis_col[iCol] = iisCol; + if (iis_lp.col_cost_[iisCol]) return lpDataOkReturn(false); HighsInt col_bound = this->col_bound_[iisCol]; bound = col_bound == kIisBoundStatusLower || col_bound == kIisBoundStatusBoxed ? lp.col_lower_[iCol] : -kHighsInf; - if (iis_lp.col_lower_[iisCol] != bound) return false; + if (iis_lp.col_lower_[iisCol] != bound) return lpDataOkReturn(false); bound = col_bound == kIisBoundStatusUpper || col_bound == kIisBoundStatusBoxed ? lp.col_upper_[iCol] : kHighsInf; - if (iis_lp.col_upper_[iisCol] != bound) return false; + if (iis_lp.col_upper_[iisCol] != bound) return lpDataOkReturn(false); } - const HighsInt illegal_index = -1; - const double illegal_value = kHighsInf; std::vector index; std::vector value; // Work through the LP matrix, checking the matrix index/value @@ -842,8 +859,9 @@ bool HighsIis::lpDataOk(const HighsLp& lp, const HighsOptions& options) const { HighsInt iRow = lp.a_matrix_.index_[iEl]; HighsInt iisRow = iis_row[iRow]; if (iisRow >= 0) { - if (index[iisRow] != iRow) return false; - if (value[iisRow] != lp.a_matrix_.value_[iEl]) return false; + if (index[iisRow] != iRow) return lpDataOkReturn(false); + if (value[iisRow] != lp.a_matrix_.value_[iEl]) + return lpDataOkReturn(false); index[iisRow] = illegal_index; value[iisRow] = illegal_value; } @@ -853,22 +871,23 @@ bool HighsIis::lpDataOk(const HighsLp& lp, const HighsOptions& options) const { for (HighsInt iisRow = 0; iisRow < iis_num_row; iisRow++) { HighsInt iRow = this->row_index_[iisRow]; // Use index/value to scatter the IIS matrix row - index.assign(iis_num_row, illegal_index); - value.assign(iis_num_row, illegal_value); + index.assign(iis_num_col, illegal_index); + value.assign(iis_num_col, illegal_value); for (HighsInt iEl = iis_lp.a_matrix_.start_[iisRow]; iEl < iis_lp.a_matrix_.start_[iisRow + 1]; iEl++) { HighsInt iisCol = iis_lp.a_matrix_.index_[iEl]; - HighsInt iCol = this->row_index_[iisCol]; + HighsInt iCol = this->col_index_[iisCol]; index[iisCol] = iCol; value[iisCol] = iis_lp.a_matrix_.value_[iEl]; } for (HighsInt iEl = lp.a_matrix_.start_[iRow]; iEl < lp.a_matrix_.start_[iRow + 1]; iEl++) { HighsInt iCol = lp.a_matrix_.index_[iEl]; - HighsInt iisCol = iis_row[iCol]; + HighsInt iisCol = iis_col[iCol]; if (iisCol >= 0) { - if (index[iisCol] != iCol) return false; - if (value[iisCol] != lp.a_matrix_.value_[iEl]) return false; + if (index[iisCol] != iCol) return lpDataOkReturn(false); + if (value[iisCol] != lp.a_matrix_.value_[iEl]) + return lpDataOkReturn(false); index[iisCol] = illegal_index; value[iisCol] = illegal_value; } @@ -894,35 +913,39 @@ bool HighsIis::lpDataOk(const HighsLp& lp, const HighsOptions& options) const { iEl < iis_lp.a_matrix_.start_[iisCol + 1]; iEl++) { HighsInt iisRow = iis_lp.a_matrix_.index_[iEl]; HighsInt iRow = this->row_index_[iisRow]; - if (index[iRow] != iisRow) return false; - if (value[iRow] != iis_lp.a_matrix_.value_[iEl]) return false; + if (index[iRow] != iisRow) return lpDataOkReturn(false); + if (value[iRow] != iis_lp.a_matrix_.value_[iEl]) + return lpDataOkReturn(false); } } } else { for (HighsInt iisRow = 0; iisRow < iis_num_row; iisRow++) { HighsInt iRow = this->row_index_[iisRow]; // Use index/value to scatter the LP matrix row - index.assign(lp.num_row_, illegal_index); - value.assign(lp.num_row_, illegal_value); + index.assign(lp.num_col_, illegal_index); + value.assign(lp.num_col_, illegal_value); for (HighsInt iEl = lp.a_matrix_.start_[iRow]; iEl < lp.a_matrix_.start_[iRow + 1]; iEl++) { HighsInt iCol = lp.a_matrix_.index_[iEl]; - HighsInt iisCol = iis_row[iCol]; + HighsInt iisCol = iis_col[iCol]; index[iCol] = iisCol; value[iCol] = lp.a_matrix_.value_[iEl]; } for (HighsInt iEl = iis_lp.a_matrix_.start_[iisRow]; iEl < iis_lp.a_matrix_.start_[iisRow + 1]; iEl++) { HighsInt iisCol = iis_lp.a_matrix_.index_[iEl]; - HighsInt iCol = this->row_index_[iisCol]; - if (index[iCol] != iisCol) return false; - if (value[iCol] != iis_lp.a_matrix_.value_[iEl]) return false; + HighsInt iCol = this->col_index_[iisCol]; + if (index[iCol] != iisCol) return lpDataOkReturn(false); + if (value[iCol] != iis_lp.a_matrix_.value_[iEl]) + return lpDataOkReturn(false); } } } - return true; + return lpDataOkReturn(true); } +bool lpOkReturn(const bool return_value) { return return_value; } + bool HighsIis::lpOk(const HighsOptions& options) const { // Check that the IIS LP is OK (infeasible and optimal if // any bound is relaxed) @@ -946,7 +969,7 @@ bool HighsIis::lpOk(const HighsOptions& options) const { if (h.getModelStatus() != HighsModelStatus::kInfeasible) { highsLogUser(log_options, HighsLogType::kError, "HighsIis: IIS LP is not infeasible\n"); - return false; + return lpOkReturn(false); } auto optimal = [&]() -> bool { if (options.log_dev_level > 0) h.writeModel(""); @@ -963,7 +986,7 @@ bool HighsIis::lpOk(const HighsOptions& options) const { "bound of %g yield IIS LP with status %s\n", int(iisCol), int(iCol), iis_lp.col_lower_[iisCol], h.modelStatusToString(h.getModelStatus()).c_str()); - return false; + return lpOkReturn(false); } h.changeColBounds(iisCol, iis_lp.col_lower_[iisCol], iis_lp.col_upper_[iisCol]); @@ -976,7 +999,7 @@ bool HighsIis::lpOk(const HighsOptions& options) const { "bound of %g yield IIS LP with status %s\n", int(iisCol), int(iCol), iis_lp.col_upper_[iisCol], h.modelStatusToString(h.getModelStatus()).c_str()); - return false; + return lpOkReturn(false); } h.changeColBounds(iisCol, iis_lp.col_lower_[iisCol], iis_lp.col_upper_[iisCol]); @@ -992,7 +1015,7 @@ bool HighsIis::lpOk(const HighsOptions& options) const { "of %g yield IIS LP with status %s\n", int(iisRow), int(iRow), iis_lp.row_lower_[iisRow], h.modelStatusToString(h.getModelStatus()).c_str()); - return false; + return lpOkReturn(false); } h.changeRowBounds(iisRow, iis_lp.row_lower_[iisRow], iis_lp.row_upper_[iisRow]); @@ -1001,15 +1024,15 @@ bool HighsIis::lpOk(const HighsOptions& options) const { h.changeRowBounds(iisRow, iis_lp.row_lower_[iisRow], kHighsInf); if (!optimal()) { highsLogUser(log_options, HighsLogType::kError, - "HighsIis: IIS rowumn %d (LP rowumn %d): relaxing upper " + "HighsIis: IIS row %d (LP row %d): relaxing upper " "bound of %g yield IIS LP with status %s\n", int(iisRow), int(iRow), iis_lp.row_upper_[iisRow], h.modelStatusToString(h.getModelStatus()).c_str()); - return false; + return lpOkReturn(false); } h.changeRowBounds(iisRow, iis_lp.row_lower_[iisRow], iis_lp.row_upper_[iisRow]); } } - return true; + return lpOkReturn(true); } diff --git a/highs/lp_data/HighsIis.h b/highs/lp_data/HighsIis.h index 92f68d8878..5e5ebad1af 100644 --- a/highs/lp_data/HighsIis.h +++ b/highs/lp_data/HighsIis.h @@ -13,7 +13,7 @@ #include "model/HighsModel.h" -const bool kIisDevReport = false; +const bool kIisDevReport = true; // false; enum IisBoundStatus { kIisBoundStatusDropped = -1, diff --git a/highs/lp_data/HighsInterface.cpp b/highs/lp_data/HighsInterface.cpp index 41c5c9614c..38f4d5b72b 100644 --- a/highs/lp_data/HighsInterface.cpp +++ b/highs/lp_data/HighsInterface.cpp @@ -1986,6 +1986,8 @@ HighsStatus Highs::getIisInterface() { this->iis_.col_status_.assign(lp.num_col_, kIisStatusNotInConflict); this->iis_.row_status_.assign(lp.num_row_, kIisStatusNotInConflict); } else { + // To get the IIS data needs the matrix to be column-wise + model_.lp_.a_matrix_.ensureColwise(); return_status = this->iis_.getData(lp, options_, basis_, infeasible_row_subset); if (return_status == HighsStatus::kOk) { @@ -2095,10 +2097,11 @@ HighsStatus Highs::elasticityFilter( const bool get_infeasible_row, std::vector& infeasible_row_subset) { // this->writeModel("infeasible.mps"); + // // Solve the feasibility relaxation problem for the given penalties, - // continuing to act as the elasticity filter get_infeasible_row is - // true, resulting in an infeasibility subset for further refinement - // as an IIS + // continuing to act as the elasticity filter if get_infeasible_row + // is true, resulting in an infeasibility subset for further + // refinement as an IIS // // Construct the e-LP: // diff --git a/highs/lp_data/Iis.md b/highs/lp_data/Iis.md new file mode 100644 index 0000000000..4e9cca5d48 --- /dev/null +++ b/highs/lp_data/Iis.md @@ -0,0 +1,60 @@ +# HiGHS irreducible infeasibility system (IIS) facility + +Further to documentation in https://ergo-code.github.io/HiGHS/stable/guide/advanced/ + +The IIS search is rooted in `Highs::getIisInterface()`, which first +checks whether the `Highs::model_status_` is +`HighsModelStatus::kOptimal` or `HighsModelStatus::kUnbounded`, in +which case the model is feasible so no IIS exists. Otherwise, the +trivial check (for inconsistent bounds or empty infeasible rows) - +performed by `HighsIis::trivial()` - and infeasible rows based on row +value bounds - performed by `HighsIis::rowValueBounds()` - are +performed. If `Highs::options_.iis_strategy` is `kIisStrategyLight` +then `Highs::getIisInterface()` returns. + +The "full" IIS calculation operates in two phases: after a set of +mutually infeasible rows has been identified, this is reduced to an +IIS. The set of mutually infeasible rows can be found in two ways. + +Firstly, if it is known that the model is infeasible, then the simplex +solver may have identified a dual ray. If there is a dual ray then its +nonzeros correspond to a set of mutually infeasible constraints. If +there is no dual ray - as might happen if the model's infeasibility +has been identified in presolve - then the incumbent model is solved +with `Highs::options_.presolve` = "off". Unfortunately the "ray route" +is not robust, so currently switched off. + +Secondly - and the only route at the moment - an elasticity filter +calculation is done to identify an infeasible subset of rows. This is +calculation performed in `HighsIis::elasticityFilter`. This method is +more general than is necessary for finding the set of mutually +infeasible rows for an IIS calculation, and can be called directly +using `Highs::feasibilityRelaxation`. + +The essence of the `HighsIis::elasticityFilter` is that it allows +lower bounds, upper bounds and RHS values to be violated. There are +penalties for doing so that can be global for each of these three +cases, or local to each column lower bound, upper bound and row +bound. The "elasticity LP" is constructed by adding elastic variables +to transform the constraints from +$$ +L <= Ax <= U;\qquad l <= x <= u +$$ +to +$$ +L <= Ax + e_L - e_U <= U;\qquad l <= x + e_l - e_u <= u, +$$ +where the elastic variables are not used if the corresponding bound is +infinite or the local/global penalty is negative. The original bounds +on the variables $x$ are removed, and the objective is the linear +function of the elastic variables given by the local/global +penalties. Note that the model modifications required to achieve this +formulation, are made by calls to methods in the `Highs` class so that +the value of any initial basis is maintained. + +For the purposes of IIS calculation, each elastic variable has a + +If the original constraints cannot be satisfied, then some + +After the elasticity LP has been solved, each elasticity variable whose optimal value is positive is fixed at zero is + diff --git a/highs/util/HFactorRefactor.cpp b/highs/util/HFactorRefactor.cpp index 723de0aef6..e404747b31 100644 --- a/highs/util/HFactorRefactor.cpp +++ b/highs/util/HFactorRefactor.cpp @@ -185,7 +185,7 @@ HighsInt HFactor::rebuild(HighsTimerClock* factor_timer_clock_pointer) { // // First of all complete the L factor with identity columns so // that FtranL counts the RHS entries in rows that don't yet have - // picots by running to completion. In the hyper-sparse code, + // pivots by running to completion. In the hyper-sparse code, // these will HOPEFULLY be skipped // // There are already l_start entries for the first stage rows, but diff --git a/highs/util/HighsSparseMatrix.cpp b/highs/util/HighsSparseMatrix.cpp index 96beced7ba..412bfab4c7 100644 --- a/highs/util/HighsSparseMatrix.cpp +++ b/highs/util/HighsSparseMatrix.cpp @@ -890,12 +890,12 @@ void HighsSparseMatrix::considerRowScaling( row_scale_value = min(max(min_allow_row_scale, row_scale_value), max_allow_row_scale); row_scale[iRow] = row_scale_value; - // Scale the rowumn + // Scale the row for (HighsInt iEl = this->start_[iRow]; iEl < this->start_[iRow + 1]; iEl++) this->value_[iEl] *= row_scale[iRow]; } else { - // Empty rowumn + // Empty row row_scale[iRow] = 1; } }