diff --git a/cmake/sources-python.cmake b/cmake/sources-python.cmake index ac1c94c7e4..33a3081c72 100644 --- a/cmake/sources-python.cmake +++ b/cmake/sources-python.cmake @@ -215,6 +215,7 @@ set(highs_sources_python highs/mip/HighsImplications.cpp highs/mip/HighsLpAggregator.cpp highs/mip/HighsLpRelaxation.cpp + highs/mip/HighsMachineSchedSeparator.cpp highs/mip/HighsMipAnalysis.cpp highs/mip/HighsMipSolver.cpp highs/mip/HighsMipSolverData.cpp @@ -339,6 +340,7 @@ set(highs_headers_python highs/mip/HighsImplications.h highs/mip/HighsLpAggregator.h highs/mip/HighsLpRelaxation.h + highs/mip/HighsMachineSchedSeparator.h highs/mip/HighsMipAnalysis.h highs/mip/HighsMipSolver.h highs/mip/HighsMipSolverData.h diff --git a/cmake/sources.cmake b/cmake/sources.cmake index 89897224c9..1b54006ef4 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -371,6 +371,7 @@ set(highs_sources mip/HighsImplications.cpp mip/HighsLpAggregator.cpp mip/HighsLpRelaxation.cpp + mip/HighsMachineSchedSeparator.cpp mip/HighsMipAnalysis.cpp mip/HighsMipSolver.cpp mip/HighsMipSolverData.cpp @@ -498,6 +499,7 @@ set(highs_headers mip/HighsImplications.h mip/HighsLpAggregator.h mip/HighsLpRelaxation.h + mip/HighsMachineSchedSeparator.h mip/HighsMipAnalysis.h mip/HighsMipSolver.h mip/HighsMipSolverData.h diff --git a/highs/meson.build b/highs/meson.build index 93b4a6f0c8..c164f195ef 100644 --- a/highs/meson.build +++ b/highs/meson.build @@ -232,6 +232,7 @@ _srcs = [ 'mip/HighsImplications.cpp', 'mip/HighsLpAggregator.cpp', 'mip/HighsLpRelaxation.cpp', + 'mip/HighsMachineSchedSeparator.cpp', 'mip/HighsMipAnalysis.cpp', 'mip/HighsMipSolver.cpp', 'mip/HighsMipSolverData.cpp', diff --git a/highs/mip/HighsMachineSchedSeparator.cpp b/highs/mip/HighsMachineSchedSeparator.cpp new file mode 100644 index 0000000000..d91b2f01e6 --- /dev/null +++ b/highs/mip/HighsMachineSchedSeparator.cpp @@ -0,0 +1,259 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* */ +/* This file is part of the HiGHS linear optimization suite */ +/* */ +/* Available as open-source under the MIT License */ +/* */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/**@file mip/HighsMachineSchedSeparator.cpp + */ + +#include "mip/HighsMachineSchedSeparator.h" + +#include + +#include "../extern/pdqsort/pdqsort.h" +#include "mip/HighsCutGeneration.h" +#include "mip/HighsLpRelaxation.h" +#include "mip/HighsMipSolverData.h" +#include "mip/HighsTransformedLp.h" + +bool HighsMachineSchedSeparator::findSingleMachineScheduleClique( + std::vector>& vals, + std::vector>& inds, std::vector& rhss, + double& releasedate, const HighsMipSolver& mipsolver) { + struct pair_hash { + size_t operator()(const std::pair& p) const { + return HighsHashHelpers::hash(p); + } + }; + enum class ArcType { + kIfBinOne, // 0 + kIfBinZero, // 1 + }; + HighsInt largestDegree = 0; + HighsInt largestDegreeCol = -1; + std::vector degrees(mipsolver.numCol()); + std::unordered_map, + std::tuple, pair_hash> + adjacency; + + auto addEntry = [&](HighsInt posCol, HighsInt negCol, HighsInt binCol, + double p, ArcType t) { + // My_ji + s_i - s_j >= p_ji + // y_ji = val -> s_i >= s_j + p_ji - M * val + // Make an arc from negCol (j) to posCol (i) + if (p < 0) return; + auto it = adjacency.find({negCol, posCol}); + if (it != adjacency.end()) { + if (std::get<0>(it->second) < p) { + std::get<0>(it->second) = p; + std::get<1>(it->second) = binCol; + std::get<2>(it->second) = t; + } + } else { + degrees[posCol]++; + if (degrees[posCol] > largestDegree || + (degrees[posCol] == largestDegree && posCol < largestDegreeCol)) { + largestDegreeCol = posCol; + largestDegree = degrees[posCol]; + } + adjacency.emplace(std::make_pair(negCol, posCol), + std::make_tuple(p, binCol, t)); + } + }; + + HighsInt numRows = 0; + const HighsInt maxRows = std::min(HighsInt{1000}, 2 * mipsolver.numRow()); + adjacency.reserve(maxRows + 2); + for (HighsInt row = 0; row != mipsolver.numRow(); row++) { + double rowLower = mipsolver.model_->row_lower_[row]; + double rowUpper = mipsolver.model_->row_upper_[row]; + if (rowLower == rowUpper) continue; + HighsInt start = mipsolver.mipdata_->ARstart_[row]; + HighsInt end = mipsolver.mipdata_->ARstart_[row + 1]; + if (end - start != 3) continue; + bool machineSchedRow = true; + HighsInt posContCol = -1; + HighsInt negContCol = -1; + HighsInt binCol = -1; + double binCoef = 0; + for (HighsInt i = start; i != end; i++) { + HighsInt col = mipsolver.mipdata_->ARindex_[i]; + if (mipsolver.mipdata_->domain.col_lower_[col] == -kHighsInf) { + // TODO: Could some jobs be modelled in reverse? + machineSchedRow = false; + break; + } + if (mipsolver.mipdata_->domain.isBinary(col)) { + if (binCol != -1) { + machineSchedRow = false; + break; + } + binCol = col; + binCoef = mipsolver.mipdata_->ARvalue_[i]; + } else if (mipsolver.mipdata_->ARvalue_[i] == -1) { + negContCol = col; + } else if (mipsolver.mipdata_->ARvalue_[i] == 1) { + posContCol = col; + } else { + machineSchedRow = false; + break; + } + } + if (!machineSchedRow || binCol == -1 || negContCol == -1 || + posContCol == -1 || posContCol == negContCol) + continue; + // We want to put the row into form: + // Mx_ji + s_i - s_j >= p_ji, p_ji >= 0 + // x_ji = 1 -> s_i >= s_j + p_ij - M + if (rowUpper != kHighsInf) { + // Given My_ij + s_i - s_j <= d + // The row becomes (after multiplying by -1): + // -My_ij + s_j - s_i >= -d + // Add implication s_j >= s_i + p_ij + M, p_ij + M > 0 when binCol = 1 + // Add implication s_j >= s_i + p_ij, p, p_ij > 0 when binCol = 0 + double rhs_0 = -rowUpper; + double rhs_1 = -rowUpper + binCoef; + // TODO: If both rhs_0 > 0 && rhs_1 > 0 then strengthen the claim, i.e., + // TODO: min{rhs_0, rhs_1} + y_ji * (max{rhs_0, rhs_1} - min{rhs_0, rhs_1} + if (rhs_0 > 0 || rhs_1 > 0) { + if (rhs_0 > 0) { + addEntry(negContCol, posContCol, binCol, rhs_0, ArcType::kIfBinZero); + } else { + addEntry(negContCol, posContCol, binCol, rhs_1, ArcType::kIfBinOne); + } + numRows++; + } + } + if (rowLower != -kHighsInf) { + // Given Mx_ij + s_i - s_j >= d + // Add implication s_i >= s_j + pji - M, p_ji - M > 0 when binCol = 1 + // Add implication s_i >= s_j + pji, p_ji > 0 when binCol = 0 + double rhs_0 = rowLower; + double rhs_1 = rowLower - binCoef; + if (rhs_0 > 0 || rhs_1 > 0) { + if (rhs_0 > 0) { + addEntry(posContCol, negContCol, binCol, rhs_0, ArcType::kIfBinZero); + } else if (rhs_1 > 0) { + addEntry(posContCol, negContCol, binCol, rhs_1, ArcType::kIfBinOne); + } + numRows++; + } + } + if (numRows >= maxRows) break; + } + + // A clique of size 3 needs at least 6 arcs + if (numRows <= 5) return false; + + // Greedily search neighbours of largest degree column for a double-sided + // clique (corresponds to a single machine schedule) + std::vector potentialNeighbours; + potentialNeighbours.reserve(largestDegree); + for (const auto& arc : adjacency) { + HighsInt col = std::get<1>(arc.first); + if (col == largestDegreeCol) { + potentialNeighbours.emplace_back(std::get<0>(arc.first)); + } + } + pdqsort(potentialNeighbours.begin(), potentialNeighbours.end(), + [&](const HighsInt c1, const HighsInt c2) { + return degrees[c1] > degrees[c2]; + }); + + std::vector neighbours; + neighbours.reserve(largestDegree + 1); + neighbours.emplace_back(largestDegreeCol); + releasedate = mipsolver.mipdata_->domain.col_lower_[largestDegreeCol]; + std::vector processingTimes; + processingTimes.resize(largestDegree + 1, kHighsInf); + // Iterate over potential neighbours and check validity + for (HighsInt col : potentialNeighbours) { + bool valid_neighbour = true; + for (HighsInt neighbour : neighbours) { + auto fromArc = adjacency.find({col, neighbour}); + auto toArc = adjacency.find({neighbour, col}); + // Need to verify that the to and fromArc exist and are opposites + if (toArc == adjacency.end() || fromArc == adjacency.end() || + std::get<1>(fromArc->second) != std::get<1>(toArc->second) || + std::get<2>(fromArc->second) == std::get<2>(toArc->second)) { + valid_neighbour = false; + break; + } + } + if (!valid_neighbour) continue; + size_t newNeighbourIndex = neighbours.size(); + // Extract the processing times from the arcs + for (size_t i = 0; i != neighbours.size(); ++i) { + HighsInt neighbour = neighbours[i]; + const auto fromArc = adjacency.find({col, neighbour}); + const auto toArc = adjacency.find({neighbour, col}); + processingTimes[newNeighbourIndex] = std::min( + processingTimes[newNeighbourIndex], std::get<0>(fromArc->second)); + processingTimes[i] = + std::min(processingTimes[i], std::get<0>(toArc->second)); + } + releasedate = + std::min(mipsolver.mipdata_->domain.col_lower_[col], releasedate); + neighbours.emplace_back(col); + } + if (neighbours.size() < 3) return false; + + // Now populate the actual inequalities + vals.resize(neighbours.size(), std::vector(neighbours.size())); + inds.resize(neighbours.size(), std::vector(neighbours.size())); + rhss.resize(neighbours.size()); + for (size_t i = 0; i != neighbours.size(); ++i) { + rhss[i] -= releasedate; + HighsInt col = neighbours[i]; + for (size_t j = 0; j != neighbours.size(); ++j) { + size_t jj = j >= i ? j + 1 : j; + if (jj >= neighbours.size()) continue; + HighsInt neighbour = neighbours[jj]; + const auto toArc = adjacency.find({neighbour, col}); + assert(toArc != adjacency.end()); + inds[i][j] = std::get<1>(toArc->second); + vals[i][j] = processingTimes[jj]; + if (std::get<2>(toArc->second) == ArcType::kIfBinZero) { + rhss[i] -= vals[i][j]; + vals[i][j] *= -1; + } + } + // Put the job start time on the LHS + inds[i].back() = col; + vals[i].back() = -1; + } + + return true; +} + +void HighsMachineSchedSeparator::separateLpSolution( + HighsLpRelaxation& lpRelaxation, HighsLpAggregator& lpAggregator, + HighsTransformedLp& transLp, HighsCutPool& cutpool) { + if (separated) return; + const HighsMipSolver& mip = lpRelaxation.getMipSolver(); + std::vector> vals; + std::vector> inds; + std::vector rhss; + double releasedate; + has_single_machine_schedule = + findSingleMachineScheduleClique(vals, inds, rhss, releasedate, mip); + if (!has_single_machine_schedule) { + separated = true; + return; + } + + // Load the cuts!!! + const std::vector& lpSolution = lpRelaxation.getSolution().col_value; + for (size_t i = 0; i != inds.size(); ++i) { + double viol = -rhss[i]; + for (size_t j = 0; j != inds[i].size(); j++) { + viol += lpSolution[inds[i][j]] * vals[i][j]; + } + if (viol >= 10 * mip.mipdata_->feastol) + cutpool.addCut(mip, inds[i].data(), vals[i].data(), inds[i].size(), + rhss[i]); + } + separated = true; +} diff --git a/highs/mip/HighsMachineSchedSeparator.h b/highs/mip/HighsMachineSchedSeparator.h new file mode 100644 index 0000000000..ae59c59259 --- /dev/null +++ b/highs/mip/HighsMachineSchedSeparator.h @@ -0,0 +1,60 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* */ +/* This file is part of the HiGHS linear optimization suite */ +/* */ +/* Available as open-source under the MIT License */ +/* */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/**@file mip/HighsMachineSchedSeparator.h + * @brief Class for separating release-date cuts, see Section 4.1 in + * Recognition and Exploitation of Single-Machine Scheduling Subproblems in + * Mixed Integer Programs, Reinout Lambertus Henricus Wijfjes, Msc Thesis, 2022, + * or Cutting Plane Algorithm for the Single Machine Scheduling Problem + * with Release Times, G. L. Nemhauser and M. W. P. Savelsbergh, 1992. + * + * Given a set of jobs N with start times s_j, processing times p_ij, + * release times r_j, and binary dependency variables y_ij, + * where y = 1 -> s_i - s_j < 0. + * To simplify the problem: We set p_j = min{p_ji : i \in N / j} + * A valid inequality is then: + * s_j >= min{r_{i} : i \in N} + \sum_{i \in N / j} p_i y_ij + * + */ + +#ifndef MIP_HIGHS_MACHINE_SCHED_SEPARATOR_H_ +#define MIP_HIGHS_MACHINE_SCHED_SEPARATOR_H_ + +#include "mip/HighsMipSolver.h" +#include "mip/HighsMipSolverData.h" +#include "mip/HighsSeparator.h" +#include "util/HighsRandom.h" + +/// Helper class to compute single-row relaxations from the current LP +/// relaxation by substituting bounds and aggregating rows +class HighsMachineSchedSeparator : public HighsSeparator { + private: + HighsRandom randgen; + bool has_single_machine_schedule = false; + bool separated = false; + + public: + void separateLpSolution(HighsLpRelaxation& lpRelaxation, + HighsLpAggregator& lpAggregator, + HighsTransformedLp& transLp, + HighsCutPool& cutpool) override; + + bool findSingleMachineScheduleClique(std::vector>& vals, + std::vector>& inds, + std::vector& rhss, + double& releasedate, + const HighsMipSolver& mipsolver); + + HighsMachineSchedSeparator(const HighsMipSolver& mipsolver) + : HighsSeparator(mipsolver, kMachineSchedSepaString) { + randgen.initialise(mipsolver.options_mip_->random_seed); + has_single_machine_schedule = false; + separated = false; + } +}; + +#endif diff --git a/highs/mip/HighsMipAnalysis.cpp b/highs/mip/HighsMipAnalysis.cpp index 2f36e01ec7..2a69c69bc8 100644 --- a/highs/mip/HighsMipAnalysis.cpp +++ b/highs/mip/HighsMipAnalysis.cpp @@ -42,6 +42,8 @@ void HighsMipAnalysis::setupMipTime(const HighsOptions& options) { std::make_pair(kPathAggrSepaString, kMipClockPathAggrSepa)); sepa_name_clock.push_back( std::make_pair(kModKSepaString, kMipClockModKSepa)); + sepa_name_clock.push_back( + std::make_pair(kMachineSchedSepaString, kMipClockMachineSchedSepa)); } } diff --git a/highs/mip/HighsSeparation.cpp b/highs/mip/HighsSeparation.cpp index a85691a1d0..fd316d4609 100644 --- a/highs/mip/HighsSeparation.cpp +++ b/highs/mip/HighsSeparation.cpp @@ -16,6 +16,7 @@ #include "mip/HighsImplications.h" #include "mip/HighsLpAggregator.h" #include "mip/HighsLpRelaxation.h" +#include "mip/HighsMachineSchedSeparator.h" #include "mip/HighsMipSolverData.h" #include "mip/HighsModkSeparator.h" #include "mip/HighsPathSeparator.h" @@ -31,6 +32,7 @@ HighsSeparation::HighsSeparation(const HighsMipSolver& mipsolver) { separators.emplace_back(new HighsTableauSeparator(mipsolver)); separators.emplace_back(new HighsPathSeparator(mipsolver)); separators.emplace_back(new HighsModkSeparator(mipsolver)); + separators.emplace_back(new HighsMachineSchedSeparator(mipsolver)); } HighsInt HighsSeparation::separationRound(HighsDomain& propdomain, diff --git a/highs/mip/HighsSeparator.h b/highs/mip/HighsSeparator.h index 75fb1a41e0..2eb5d60cdc 100644 --- a/highs/mip/HighsSeparator.h +++ b/highs/mip/HighsSeparator.h @@ -22,6 +22,8 @@ const std::string kCliqueSepaString = "Separation: Clique"; const std::string kTableauSepaString = "Separation: Tableau"; const std::string kPathAggrSepaString = "Separation: Path aggregation"; const std::string kModKSepaString = "Separation: Mod-k"; +const std::string kMachineSchedSepaString = + "Separation: Single Machine Scheduling"; class HighsLpRelaxation; class HighsTransformedLp; diff --git a/highs/mip/MipTimer.h b/highs/mip/MipTimer.h index f1044de065..4f82373190 100644 --- a/highs/mip/MipTimer.h +++ b/highs/mip/MipTimer.h @@ -95,6 +95,7 @@ enum iClockMip { kMipClockTableauSepa, kMipClockPathAggrSepa, kMipClockModKSepa, + kMipClockMachineSchedSepa, // LP solves kMipClockDuSimplexBasisSolveLp, @@ -233,6 +234,8 @@ class MipTimer { timer_pointer->clock_def(kPathAggrSepaString.c_str()); clock[kMipClockModKSepa] = timer_pointer->clock_def(kModKSepaString.c_str()); + clock[kMipClockMachineSchedSepa] = + timer_pointer->clock_def(kMachineSchedSepaString.c_str()); // Presolve - Should correspond to kMipClockRunPresolve clock[kMipClockProbingPresolve] = @@ -488,7 +491,7 @@ class MipTimer { void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) { const std::vector mip_clock_list{ kMipClockImplboundSepa, kMipClockCliqueSepa, kMipClockTableauSepa, - kMipClockPathAggrSepa, kMipClockModKSepa}; + kMipClockPathAggrSepa, kMipClockModKSepa, kMipClockMachineSchedSepa}; reportMipClockList("MipSeparation", mip_clock_list, mip_timer_clock, kMipClockTotal); //, tolerance_percent_report); }; diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 1ebf6d9ca5..107bd95ff1 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -4466,12 +4466,13 @@ HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack, // check if variable is implied integer HPRESOLVE_CHECKED_CALL(static_cast(convertImpliedInteger(col))); - // shift integral variables to have a lower bound of zero + // shift "binary" variables to have a lower bound of zero if (model->integrality_[col] != HighsVarType::kContinuous && model->col_lower_[col] != 0.0 && (model->col_lower_[col] != -kHighsInf || model->col_upper_[col] != kHighsInf) && - model->col_upper_[col] - model->col_lower_[col] > 0.5) { + model->col_upper_[col] - model->col_lower_[col] > 0.5 && + model->col_upper_[col] - model->col_lower_[col] < 1.5) { // substitute with the bound that is smaller in magnitude and only // substitute if bound is not large for an integer if (std::abs(model->col_upper_[col]) > std::abs(model->col_lower_[col])) {