diff --git a/cpp/cuopt_cli.cpp b/cpp/cuopt_cli.cpp index ac568e07cf..4c82c26154 100644 --- a/cpp/cuopt_cli.cpp +++ b/cpp/cuopt_cli.cpp @@ -135,7 +135,6 @@ int run_single_file(const std::string& file_path, std::make_unique>(); } - // Populate the problem from MPS data model cuopt::linear_programming::populate_from_mps_data_model(problem_interface.get(), mps_data_model); const bool is_mip = (problem_interface->get_problem_category() == diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 06eacb3408..aee5e4c3cb 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -81,6 +81,7 @@ #define CUOPT_PRESOLVE_FILE "presolve_file" #define CUOPT_RANDOM_SEED "random_seed" #define CUOPT_PDLP_PRECISION "pdlp_precision" +#define CUOPT_SC_BIG_M "sc_big_m" #define CUOPT_MIP_HYPER_HEURISTIC_POPULATION_SIZE "mip_hyper_heuristic_population_size" #define CUOPT_MIP_HYPER_HEURISTIC_NUM_CPUFJ_THREADS "mip_hyper_heuristic_num_cpufj_threads" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 14c4d227bc..18ad86da34 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -31,6 +31,14 @@ struct benchmark_info_t { template class solver_settings_t; +template +class mip_solver_settings_t; + +namespace detail { +template +struct mip_solver_settings_accessor; +} // namespace detail + template class mip_solver_settings_t { public: @@ -86,6 +94,7 @@ class mip_solver_settings_t { f_t time_limit = std::numeric_limits::infinity(); f_t work_limit = std::numeric_limits::infinity(); + f_t sc_big_m = f_t(1e5); i_t node_limit = std::numeric_limits::max(); bool heuristics_only = false; i_t reliability_branching = -1; @@ -147,6 +156,19 @@ class mip_solver_settings_t { std::vector mip_callbacks_; friend class solver_settings_t; + friend struct detail::mip_solver_settings_accessor; }; +namespace detail { + +template +struct mip_solver_settings_accessor { + static void clear_mip_callbacks(mip_solver_settings_t& settings) + { + settings.mip_callbacks_.clear(); + } +}; + +} // namespace detail + } // namespace cuopt::linear_programming diff --git a/cpp/include/cuopt/linear_programming/optimization_problem_interface.hpp b/cpp/include/cuopt/linear_programming/optimization_problem_interface.hpp index 767e62e746..afadbcaf7d 100644 --- a/cpp/include/cuopt/linear_programming/optimization_problem_interface.hpp +++ b/cpp/include/cuopt/linear_programming/optimization_problem_interface.hpp @@ -20,7 +20,7 @@ namespace cuopt::linear_programming { -enum class var_t { CONTINUOUS = 0, INTEGER }; +enum class var_t { CONTINUOUS = 0, INTEGER, SEMI_CONTINUOUS }; enum class problem_category_t : int8_t { LP = 0, MIP = 1, IP = 2 }; template diff --git a/cpp/include/cuopt/linear_programming/optimization_problem_utils.hpp b/cpp/include/cuopt/linear_programming/optimization_problem_utils.hpp index 90e853f530..741824485f 100644 --- a/cpp/include/cuopt/linear_programming/optimization_problem_utils.hpp +++ b/cpp/include/cuopt/linear_programming/optimization_problem_utils.hpp @@ -87,9 +87,14 @@ void populate_from_mps_data_model(optimization_problem_interface_t* pr if (!char_variable_types.empty()) { std::vector enum_variable_types(char_variable_types.size()); for (size_t i = 0; i < char_variable_types.size(); ++i) { - enum_variable_types[i] = (char_variable_types[i] == 'I' || char_variable_types[i] == 'B') - ? var_t::INTEGER - : var_t::CONTINUOUS; + const char c = char_variable_types[i]; + if (c == 'I' || c == 'B') { + enum_variable_types[i] = var_t::INTEGER; + } else if (c == 'S') { + enum_variable_types[i] = var_t::SEMI_CONTINUOUS; + } else { + enum_variable_types[i] = var_t::CONTINUOUS; + } } problem->set_variable_types(enum_variable_types.data(), enum_variable_types.size()); // Problem category (LP/MIP/IP) is auto-detected by set_variable_types @@ -253,7 +258,9 @@ void populate_from_data_model_view(optimization_problem_interface_t* p data_model->get_variable_types().data() + data_model->get_variable_types().size(), enum_variable_types.begin(), [](const auto val) -> var_t { - return (val == 'I' || val == 'B') ? var_t::INTEGER : var_t::CONTINUOUS; + if (val == 'I' || val == 'B') return var_t::INTEGER; + if (val == 'S') return var_t::SEMI_CONTINUOUS; + return var_t::CONTINUOUS; }); problem->set_variable_types(enum_variable_types.data(), enum_variable_types.size()); // Problem category (LP/MIP/IP) is auto-detected by set_variable_types diff --git a/cpp/libmps_parser/src/mps_parser.cpp b/cpp/libmps_parser/src/mps_parser.cpp index 586544331f..1eda61d1c4 100644 --- a/cpp/libmps_parser/src/mps_parser.cpp +++ b/cpp/libmps_parser/src/mps_parser.cpp @@ -243,14 +243,14 @@ BoundType convert(std::string_view str) return LowerBoundIntegerVariable; } else if (str == "UI") { return UpperBoundIntegerVariable; - } else if (str == "LC") { - return SemiContiniousVariable; + } else if (str == "SC" || str == "LC") { + return SemiContinuousVariable; } else { mps_parser_expects(false, error_type_t::ValidationError, "Invalid variable bound type found in BOUNDS section! Bound type=%s", std::string(str).c_str()); - return SemiContiniousVariable; + return SemiContinuousVariable; } } @@ -1378,6 +1378,7 @@ void mps_parser_t::read_bound_and_value(std::string_view line, switch (bound_type) { case LowerBound: { variable_lower_bounds[var_id] = get_numerical_bound(line, start); + lower_bounds_defined_for_var_id.insert(var_id); break; } case UpperBound: { @@ -1394,15 +1395,18 @@ void mps_parser_t::read_bound_and_value(std::string_view line, const f_t val = get_numerical_bound(line, start); variable_lower_bounds[var_id] = val; variable_upper_bounds[var_id] = val; + lower_bounds_defined_for_var_id.insert(var_id); break; } case Free: { variable_lower_bounds[var_id] = -std::numeric_limits::infinity(); variable_upper_bounds[var_id] = +std::numeric_limits::infinity(); + lower_bounds_defined_for_var_id.insert(var_id); break; } case LowerBoundNegInf: variable_lower_bounds[var_id] = -std::numeric_limits::infinity(); + lower_bounds_defined_for_var_id.insert(var_id); break; case UpperBoundInf: variable_upper_bounds[var_id] = +std::numeric_limits::infinity(); @@ -1411,6 +1415,7 @@ void mps_parser_t::read_bound_and_value(std::string_view line, variable_lower_bounds[var_id] = 0; variable_upper_bounds[var_id] = 1; var_types[var_id] = 'I'; + lower_bounds_defined_for_var_id.insert(var_id); break; case LowerBoundIntegerVariable: // CPLEX MPS file references seems to imply that integer variables default to an upper bound @@ -1420,6 +1425,7 @@ void mps_parser_t::read_bound_and_value(std::string_view line, } variable_lower_bounds[var_id] = get_numerical_bound(line, start); var_types[var_id] = 'I'; + lower_bounds_defined_for_var_id.insert(var_id); break; case UpperBoundIntegerVariable: variable_upper_bounds[var_id] = get_numerical_bound(line, start); @@ -1431,11 +1437,24 @@ void mps_parser_t::read_bound_and_value(std::string_view line, } var_types[var_id] = 'I'; break; - case SemiContiniousVariable: - mps_parser_expects(false, - error_type_t::ValidationError, - "Unsupported semi continous bound type found! Line=%s", - std::string(line).c_str()); + case SemiContinuousVariable: + // SC bound type: value is the upper bound U. + if (fixed_mps_format) { + const auto maybe_value = + start == std::string_view::npos ? std::string_view{} : trim(line.substr(start, 12)); + variable_upper_bounds[var_id] = maybe_value.empty() ? +std::numeric_limits::infinity() + : get_numerical_bound(line, start); + } else { + const auto maybe_value = + start == std::string_view::npos ? std::string_view{} : trim(line.substr(start)); + if (!maybe_value.empty()) { + std::stringstream ss{std::string(maybe_value)}; + ss >> variable_upper_bounds[var_id]; + } else { + variable_upper_bounds[var_id] = +std::numeric_limits::infinity(); + } + } + var_types[var_id] = 'S'; break; default: mps_parser_expects(false, diff --git a/cpp/libmps_parser/src/mps_parser.hpp b/cpp/libmps_parser/src/mps_parser.hpp index facad14c66..9db0ff4605 100644 --- a/cpp/libmps_parser/src/mps_parser.hpp +++ b/cpp/libmps_parser/src/mps_parser.hpp @@ -41,7 +41,7 @@ enum BoundType { BinaryVariable, LowerBoundIntegerVariable, UpperBoundIntegerVariable, - SemiContiniousVariable, + SemiContinuousVariable, }; // enum BoundType /** @@ -135,6 +135,7 @@ class mps_parser_t { std::unordered_map var_names_map{}; std::unordered_set ignored_objective_names{}; std::unordered_set bounds_defined_for_var_id{}; + std::unordered_set lower_bounds_defined_for_var_id{}; static constexpr f_t unset_range_value = std::numeric_limits::infinity(); /* Reads an MPS input file into a buffer. diff --git a/cpp/libmps_parser/tests/mps_parser_test.cpp b/cpp/libmps_parser/tests/mps_parser_test.cpp index f915fb2df5..b0ec17d149 100644 --- a/cpp/libmps_parser/tests/mps_parser_test.cpp +++ b/cpp/libmps_parser/tests/mps_parser_test.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -422,6 +423,33 @@ TEST(mps_bounds, upper_inf_var_bound) EXPECT_EQ(std::numeric_limits::infinity(), mps.variable_upper_bounds[1]); } +TEST(mps_bounds, semi_continuous_var_bounds_from_dataset) +{ + struct Case { + const char* file; + int n_vars; + double lower; + double upper; + }; + const std::vector cases = { + {"mip/sc_standard.mps", 2, 2.0, 10.0}, + {"mip/sc_lb_zero.mps", 2, 0.0, 10.0}, + {"mip/sc_no_ub.mps", 2, 2.0, 1e30}, + }; + + for (const auto& c : cases) { + SCOPED_TRACE(c.file); + auto mps = read_from_mps(c.file, false); + + ASSERT_EQ(c.n_vars, static_cast(mps.var_types.size())); + EXPECT_EQ('S', mps.var_types[0]); + ASSERT_EQ(c.n_vars, static_cast(mps.variable_lower_bounds.size())); + ASSERT_EQ(c.n_vars, static_cast(mps.variable_upper_bounds.size())); + EXPECT_DOUBLE_EQ(c.lower, mps.variable_lower_bounds[0]); + EXPECT_DOUBLE_EQ(c.upper, mps.variable_upper_bounds[0]); + } +} + TEST(mps_ranges, fixed_ranges) { std::string file = "linear_programming/good-mps-fixed-ranges.mps"; diff --git a/cpp/src/grpc/client/solve_remote.cpp b/cpp/src/grpc/client/solve_remote.cpp index 19908557e8..fb39a6d184 100644 --- a/cpp/src/grpc/client/solve_remote.cpp +++ b/cpp/src/grpc/client/solve_remote.cpp @@ -20,6 +20,8 @@ #include #include +#include + namespace cuopt::linear_programming { // Buffer added to the solver's time_limit to account for worker startup, @@ -209,6 +211,15 @@ std::unique_ptr> solve_mip_remote( // Check if user has set incumbent callbacks auto mip_callbacks = settings.get_mip_callbacks(); + const auto var_types = cpu_problem.get_variable_types_host(); + const bool has_sc_variables = + thrust::count(var_types.begin(), var_types.end(), var_t::SEMI_CONTINUOUS) > 0; + if (has_sc_variables && !mip_callbacks.empty()) { + CUOPT_LOG_WARN( + "Disabling remote MIP get/set callbacks: semi-continuous models are not " + "supported with callbacks"); + mip_callbacks.clear(); + } bool has_incumbents = !mip_callbacks.empty(); bool enable_tracking = has_incumbents; diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index c23b1d27ca..148eba603c 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -113,6 +113,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_HEURISTIC_INITIAL_INFEASIBILITY_WEIGHT, &mip_settings.heuristic_params.initial_infeasibility_weight, f_t(1e-9), std::numeric_limits::infinity(), f_t(1000.0), "constraint violation penalty seed"}, {CUOPT_MIP_HYPER_HEURISTIC_RELAXED_LP_TIME_LIMIT, &mip_settings.heuristic_params.relaxed_lp_time_limit, f_t(1e-9), std::numeric_limits::infinity(), f_t(1.0), "base relaxed LP time cap in heuristics"}, {CUOPT_MIP_HYPER_HEURISTIC_RELATED_VARS_TIME_LIMIT, &mip_settings.heuristic_params.related_vars_time_limit, f_t(1e-9), std::numeric_limits::infinity(), f_t(30.0), "time for related-variable structure build"}, + {CUOPT_SC_BIG_M, &mip_settings.sc_big_m, f_t(1.0), std::numeric_limits::infinity(), f_t(1e5), "big-M value for semi-continuous variables with no finite upper bound"}, }; // Int parameters diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 13649682a6..9d5ef320f2 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -36,6 +36,7 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/local_search/line_segment_search/line_segment_search.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_update_data.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/semi_continuous.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conditional_bound_strengthening.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/multi_probe.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/probing_cache.cu diff --git a/cpp/src/mip_heuristics/presolve/semi_continuous.cu b/cpp/src/mip_heuristics/presolve/semi_continuous.cu new file mode 100644 index 0000000000..7b710e7a63 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/semi_continuous.cu @@ -0,0 +1,264 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "semi_continuous.cuh" + +#include "bounds_presolve.cuh" + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace cuopt::linear_programming::detail { + +namespace { + +constexpr double sc_infinity_threshold = 1e30; + +template +bool is_effectively_infinite_sc_upper_bound(f_t ub) +{ + return !std::isfinite(ub) || ub >= static_cast(sc_infinity_threshold); +} + +} // namespace + +template +bool reformulate_semi_continuous(optimization_problem_t& op_problem, + const mip_solver_settings_t& settings, + std::vector* used_fallback_big_m) +{ + // 1. Identify semi-continuous variables + auto var_types = op_problem.get_variable_types_host(); + auto var_ub = op_problem.get_variable_upper_bounds_host(); + std::vector sc_indices; + bool normalized_large_sc_ub = false; + for (i_t i = 0; i < static_cast(var_types.size()); ++i) { + if (var_types[i] != var_t::SEMI_CONTINUOUS) { continue; } + sc_indices.push_back(i); + if (is_effectively_infinite_sc_upper_bound(var_ub[i])) { + CUOPT_LOG_WARN( + "Semi-continuous var %d upper bound %.6g exceeds semi-continuous infinity " + "threshold %.6g; treating it as +inf", + i, + static_cast(var_ub[i]), + sc_infinity_threshold); + var_ub[i] = std::numeric_limits::infinity(); + normalized_large_sc_ub = true; + } + } + if (sc_indices.empty()) { return false; } + if (normalized_large_sc_ub) { + op_problem.set_variable_upper_bounds(var_ub.data(), var_ub.size()); + } + + const i_t n_orig = op_problem.get_n_variables(); + const i_t n_sc = static_cast(sc_indices.size()); + const auto* handle_ptr = op_problem.get_handle_ptr(); + const f_t big_m = settings.sc_big_m; + if (used_fallback_big_m != nullptr) { used_fallback_big_m->assign(n_orig, uint8_t{0}); } + + CUOPT_LOG_INFO("Reformulating %d semi-continuous variable(s) before presolve", n_sc); + + // 2. Build a relaxed copy where SC vars become continuous [0, original_ub]. + // This lets GPU bounds propagation derive tight upper bounds from the + // constraint structure without the binary domain {0} ∪ [L, U]. + optimization_problem_t op_relaxed(op_problem); + { + auto relaxed_types = var_types; + auto relaxed_ub = var_ub; + auto relaxed_lb = op_problem.get_variable_lower_bounds_host(); + for (i_t idx : sc_indices) { + relaxed_types[idx] = var_t::CONTINUOUS; + // Relax to the convex hull of {0} U [L, U] before running GPU bound propagation. + relaxed_lb[idx] = std::min(f_t(0), relaxed_lb[idx]); + if (std::isfinite(relaxed_ub[idx])) { relaxed_ub[idx] = std::max(f_t(0), relaxed_ub[idx]); } + } + op_relaxed.set_variable_types(relaxed_types.data(), n_orig); + op_relaxed.set_variable_lower_bounds(relaxed_lb.data(), n_orig); + op_relaxed.set_variable_upper_bounds(relaxed_ub.data(), n_orig); + } + + // 3. Run GPU bounds propagation on the relaxed problem to tighten UBs. + // Skip propagation when there are no constraints (nothing to propagate). + auto tight_ub = var_ub; // fallback: normalized original UBs + + if (op_relaxed.get_n_constraints() > 0) { + problem_t temp_pb(op_relaxed, settings.get_tolerances()); + mip_solver_context_t ctx(handle_ptr, &temp_pb, settings); + + typename bound_presolve_t::settings_t bp_settings; + bp_settings.time_limit = 5.0; + bound_presolve_t bps(ctx, bp_settings); + bps.resize(temp_pb); + bps.solve(temp_pb); + + // Copy tightened upper bounds from GPU to host + raft::copy(tight_ub.data(), bps.upd.ub.data(), n_orig, handle_ptr->get_stream()); + handle_ptr->sync_stream(); + } + + // 4. Fetch all host arrays we need to extend with the new binary variables + // and linking constraints. + auto var_lb = op_problem.get_variable_lower_bounds_host(); + auto obj_c = op_problem.get_objective_coefficients_host(); + auto A_vals = op_problem.get_constraint_matrix_values_host(); + auto A_idx = op_problem.get_constraint_matrix_indices_host(); + auto A_off = op_problem.get_constraint_matrix_offsets_host(); + auto clb = op_problem.get_constraint_lower_bounds_host(); + auto cub = op_problem.get_constraint_upper_bounds_host(); + + // Optional arrays — only extend if they were originally set + auto b_rhs = op_problem.get_constraint_bounds_host(); + auto row_types_h = op_problem.get_row_types_host(); + + // Ensure objective and variable arrays are sized to n_orig + if (obj_c.empty()) { obj_c.assign(n_orig, f_t(0)); } + + // 5. Count how many SC vars truly need the binary-variable reformulation. + // If 0 is already inside [L, U], then "x=0 OR L<=x<=U" simplifies to + // plain continuous [L, U] — no binary needed. + std::vector needs_binary(n_sc, true); + i_t n_binary_needed = 0; + for (i_t s = 0; s < n_sc; ++s) { + const i_t idx = sc_indices[s]; + needs_binary[s] = + !(var_lb[idx] <= f_t(0) && std::isfinite(var_ub[idx]) && var_ub[idx] >= f_t(0)) && + !(var_lb[idx] <= f_t(0) && !std::isfinite(var_ub[idx])); + if (needs_binary[s]) { ++n_binary_needed; } + } + + // Extend variable arrays (one binary per SC var that actually needs it) + var_types.resize(n_orig + n_binary_needed, var_t::INTEGER); + var_lb.resize(n_orig + n_binary_needed, f_t(0)); + var_ub.resize(n_orig + n_binary_needed, f_t(1)); + obj_c.resize(n_orig + n_binary_needed, f_t(0)); + + // 6. For each SC variable: derive U when needed, then either add binary + 2 + // linking constraints or simply relax to continuous if 0 is already in + // the interval [L, U]. + i_t binary_count = 0; + for (i_t s = 0; s < n_sc; ++s) { + const i_t idx = sc_indices[s]; + const f_t L = var_lb[idx]; + const f_t orig_u = var_ub[idx]; + + if (!needs_binary[s]) { + // 0 already lies in [L, U], so the SC disjunction is just the interval itself. + CUOPT_LOG_DEBUG( + "SC var %d interval [%.6g, %.6g] already contains 0; treating it as continuous", + idx, + L, + orig_u); + var_types[idx] = var_t::CONTINUOUS; + continue; + } + + // Use GPU-propagated upper bound for positive-side SC variables when available. + // For negative-side intervals, keep the original upper bound because the relaxed + // convex hull includes 0 and is not useful for tightening the negative upper edge. + f_t U = orig_u; + if (orig_u >= f_t(0) || !std::isfinite(orig_u)) { U = tight_ub[idx]; } + if (!std::isfinite(orig_u) && std::isfinite(U)) { + CUOPT_LOG_DEBUG( + "Semi-continuous var %d upper bound was tightened from %.6g to %.6g by " + "bounds strengthening", + idx, + static_cast(orig_u), + static_cast(U)); + } + if (!std::isfinite(U)) { U = orig_u; } + if (!std::isfinite(U)) { + cuopt_assert(std::isfinite(big_m) && big_m >= L, + "Semi-continuous fallback sc_big_m must be finite and >= lower bound"); + U = big_m; + CUOPT_LOG_DEBUG( + "Semi-continuous var %d has no finite upper bound after bounds " + "strengthening; using fallback sc_big_m %.6g", + idx, + static_cast(big_m)); + if (used_fallback_big_m != nullptr) { (*used_fallback_big_m)[idx] = uint8_t{1}; } + } + + CUOPT_LOG_DEBUG("SC var %d: L=%.6g, U=%.6g (after propagation)", idx, L, U); + + const i_t b_idx = n_orig + binary_count; + ++binary_count; + + // Convert SC var to the convex hull of {0} U [L, U]. + var_types[idx] = var_t::CONTINUOUS; + var_lb[idx] = std::min(f_t(0), L); + var_ub[idx] = std::max(f_t(0), U); + + // Constraint 1: x_i - L * b_i >= 0 (clb=0, cub=+inf) + A_vals.push_back(f_t(1)); + A_idx.push_back(idx); + A_vals.push_back(-L); + A_idx.push_back(b_idx); + A_off.push_back(A_off.back() + 2); + clb.push_back(f_t(0)); + cub.push_back(std::numeric_limits::infinity()); + if (!b_rhs.empty()) { b_rhs.push_back(f_t(0)); } + if (!row_types_h.empty()) { row_types_h.push_back('G'); } + + // Constraint 2: x_i - U * b_i <= 0 (clb=-inf, cub=0) + A_vals.push_back(f_t(1)); + A_idx.push_back(idx); + A_vals.push_back(-U); + A_idx.push_back(b_idx); + A_off.push_back(A_off.back() + 2); + clb.push_back(-std::numeric_limits::infinity()); + cub.push_back(f_t(0)); + if (!b_rhs.empty()) { b_rhs.push_back(f_t(0)); } + if (!row_types_h.empty()) { row_types_h.push_back('L'); } + } + + // 7. Rebuild op_problem with the extended data. + const i_t new_n_vars = n_orig + n_binary_needed; + const i_t new_n_cons = static_cast(clb.size()); + const i_t new_nnz = static_cast(A_vals.size()); + const i_t added_constraints = 2 * n_binary_needed; + + CUOPT_LOG_INFO("SC reformulation added %d variable(s) and %d constraint(s)", + n_binary_needed, + added_constraints); + + op_problem.set_objective_coefficients(obj_c.data(), new_n_vars); + op_problem.set_variable_lower_bounds(var_lb.data(), new_n_vars); + op_problem.set_variable_upper_bounds(var_ub.data(), new_n_vars); + op_problem.set_variable_types(var_types.data(), new_n_vars); + op_problem.set_csr_constraint_matrix( + A_vals.data(), new_nnz, A_idx.data(), new_nnz, A_off.data(), new_n_cons + 1); + op_problem.set_constraint_lower_bounds(clb.data(), new_n_cons); + op_problem.set_constraint_upper_bounds(cub.data(), new_n_cons); + if (!b_rhs.empty()) { op_problem.set_constraint_bounds(b_rhs.data(), new_n_cons); } + if (!row_types_h.empty()) { op_problem.set_row_types(row_types_h.data(), new_n_cons); } + + return true; +} + +#if MIP_INSTANTIATE_FLOAT +template bool reformulate_semi_continuous(optimization_problem_t&, + const mip_solver_settings_t&, + std::vector*); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template bool reformulate_semi_continuous(optimization_problem_t&, + const mip_solver_settings_t&, + std::vector*); +#endif + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/presolve/semi_continuous.cuh b/cpp/src/mip_heuristics/presolve/semi_continuous.cuh new file mode 100644 index 0000000000..f57ef114f3 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/semi_continuous.cuh @@ -0,0 +1,41 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include + +namespace cuopt::linear_programming::detail { + +/** + * @brief Reformulate semi-continuous variables in-place inside the MIP solver. + * + * A semi-continuous variable x satisfies: x = 0 OR L <= x <= U (0 < L < U). + * Reformulation introduces a binary variable b and two linking constraints: + * x - L * b >= 0 (forces x >= L when b=1; allows x=0 when b=0) + * x - U * b <= 0 (forces x <= U when b=1; forces x=0 when b=0) + * b in {0, 1}, x in [0, U] + * + * GPU bounds propagation (bound_presolve_t) is used to derive tight upper bounds + * for SC variables that have infinite original upper bounds. If propagation cannot + * derive a finite bound, settings.sc_big_m is used as a fallback. + * + * This must be called before problem_t construction and Papilo presolve. + * + * @tparam i_t Integer index type + * @tparam f_t Floating-point value type + * @param[in,out] op_problem The optimization problem (modified in-place) + * @param[in] settings MIP solver settings (provides sc_big_m and tolerances) + * @returns true if any semi-continuous variables were found and reformulated. + */ +template +bool reformulate_semi_continuous(optimization_problem_t& op_problem, + const mip_solver_settings_t& settings, + std::vector* used_fallback_big_m = nullptr); + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index be01516657..d061db6fbd 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -23,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +50,8 @@ #include +#include + namespace cuopt::linear_programming { // This serves as both a warm up but also a mandatory initial call to setup cuSparse and cuBLAS @@ -60,6 +64,38 @@ static void init_handler(const raft::handle_t* handle_ptr) handle_ptr->get_cusparse_handle(), CUSPARSE_POINTER_MODE_DEVICE, handle_ptr->get_stream())); } +template +static void normalize_zero_lb_semi_continuous(optimization_problem_t& op_problem) +{ + if (op_problem.get_variable_types().is_empty()) { return; } + + auto var_types = op_problem.get_variable_types_host(); + auto var_lb = op_problem.get_variable_lower_bounds_host(); + bool modified = false; + + for (i_t i = 0; i < static_cast(var_types.size()); ++i) { + if (var_types[i] == var_t::SEMI_CONTINUOUS && var_lb[i] == f_t(0)) { + CUOPT_LOG_WARN( + "SC var %d has zero lower bound; treating it as continuous to match common solver behavior", + i); + var_types[i] = var_t::CONTINUOUS; + modified = true; + } + } + + if (modified) { op_problem.set_variable_types(var_types.data(), var_types.size()); } +} + +template +static bool has_semi_continuous_variables(const optimization_problem_t& op_problem) +{ + if (op_problem.get_variable_types().is_empty()) { return false; } + + auto var_types = op_problem.get_variable_types_host(); + return std::any_of( + var_types.begin(), var_types.end(), [](var_t type) { return type == var_t::SEMI_CONTINUOUS; }); +} + template static void invoke_solution_callbacks( const std::vector& mip_callbacks, @@ -281,9 +317,31 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, raft::common::nvtx::range fun_scope("Running solver"); // This is required as user might forget to set some fields + normalize_zero_lb_semi_continuous(op_problem); + if (has_semi_continuous_variables(op_problem) && !settings.initial_solutions.empty()) { + CUOPT_LOG_WARN( + "Ignoring %zu user initial solution(s): semi-continuous warm starts are not supported yet", + settings.initial_solutions.size()); + settings.initial_solutions.clear(); + } problem_checking_t::check_problem_representation(op_problem); problem_checking_t::check_initial_solution_representation(op_problem, settings); + // Reformulate semi-continuous variables (x = 0 OR L <= x <= U) before Papilo presolve. + // Uses GPU bounds propagation to derive tight upper bounds for SC vars with infinite UB. + // Track n_orig so that auxiliary binary variables added by reformulation can be stripped + // from the solution before returning it to the caller. + const i_t n_orig_before_sc = op_problem.get_n_variables(); + std::vector sc_used_fallback_big_m; + const bool had_sc = + detail::reformulate_semi_continuous(op_problem, settings, &sc_used_fallback_big_m); + if (had_sc && !settings.get_mip_callbacks().empty()) { + CUOPT_LOG_WARN( + "Disabling MIP get/set callbacks: semi-continuous models are not supported " + "with callbacks"); + detail::mip_solver_settings_accessor::clear_mip_callbacks(settings); + } + CUOPT_LOG_INFO( "Solving a problem with %d constraints, %d variables (%d integers), and %d nonzeros", op_problem.get_n_constraints(), @@ -325,7 +383,7 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, } } if (run_presolve && has_set_solution_callback) { - CUOPT_LOG_WARN("Presolve is disabled because set_solution callbacks are provided."); + CUOPT_LOG_INFO("Presolve is disabled because set_solution callbacks are provided."); run_presolve = false; } @@ -544,6 +602,31 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, } } + // Strip auxiliary binary variables that were injected by SC reformulation. + // The caller only knows about the original n_orig_before_sc variables. + if (had_sc && sol.get_solution().size() > static_cast(n_orig_before_sc)) { + sol.get_solution().resize(n_orig_before_sc, op_problem.get_handle_ptr()->get_stream()); + } + + if (had_sc && (sol.get_termination_status() == mip_termination_status_t::FeasibleFound || + sol.get_termination_status() == mip_termination_status_t::Optimal)) { + auto host_solution = + cuopt::host_copy(sol.get_solution(), op_problem.get_handle_ptr()->get_stream()); + const f_t active_tol = settings.tolerances.integrality_tolerance; + i_t num_active_fallback_big_m = 0; + for (i_t i = 0; i < static_cast(sc_used_fallback_big_m.size()); ++i) { + if (!sc_used_fallback_big_m[i]) { continue; } + if (host_solution[i] >= settings.sc_big_m - active_tol) { ++num_active_fallback_big_m; } + } + if (num_active_fallback_big_m > 0) { + return mip_solution_t{ + cuopt::logic_error("Semi-continuous solution is active at fallback sc_big_m; result may " + "depend on an artificial upper bound", + cuopt::error_type_t::RuntimeError), + op_problem.get_handle_ptr()->get_stream()}; + } + } + if (sol.get_termination_status() == mip_termination_status_t::FeasibleFound || sol.get_termination_status() == mip_termination_status_t::Optimal) { sol.log_detailed_summary(); @@ -553,6 +636,7 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, CUOPT_LOG_INFO("Writing solution to file %s", settings.sol_file.c_str()); sol.write_to_sol_file(settings.sol_file, op_problem.get_handle_ptr()->get_stream()); } + return sol; } catch (const cuopt::logic_error& e) { CUOPT_LOG_ERROR("Error in solve_mip: %s", e.what()); diff --git a/cpp/src/pdlp/optimization_problem.cu b/cpp/src/pdlp/optimization_problem.cu index 87ff9dab08..c0ecdcb41d 100644 --- a/cpp/src/pdlp/optimization_problem.cu +++ b/cpp/src/pdlp/optimization_problem.cu @@ -233,14 +233,19 @@ void optimization_problem_t::set_variable_types(const var_t* variable_ variable_types_.resize(size, stream_view_); raft::copy(variable_types_.data(), variable_types, size, stream_view_); - // Auto-detect problem category based on variable types - i_t n_integer = thrust::count_if(handle_ptr_->get_thrust_policy(), - variable_types_.begin(), - variable_types_.end(), - [] __device__(auto val) { return val == var_t::INTEGER; }); - if (n_integer == size) { + // Auto-detect problem category based on variable types. + // SEMI_CONTINUOUS vars will be reformulated into binary + continuous before solving, + // so a problem with only SC vars is treated as MIP. + i_t n_discrete = thrust::count_if( + handle_ptr_->get_thrust_policy(), + variable_types_.begin(), + variable_types_.end(), + [] __device__(auto val) { + return val == var_t::INTEGER || val == var_t::SEMI_CONTINUOUS; + }); + if (n_discrete == size) { problem_category_ = problem_category_t::IP; - } else if (n_integer > 0) { + } else if (n_discrete > 0) { problem_category_ = problem_category_t::MIP; } else { problem_category_ = problem_category_t::LP; diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 29a7f32db6..46cc474b76 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -1596,7 +1596,8 @@ optimization_problem_solution_t solve_lp( template cuopt::linear_programming::optimization_problem_t mps_data_model_to_optimization_problem( - raft::handle_t const* handle_ptr, const cuopt::mps_parser::mps_data_model_t& data_model) + raft::handle_t const* handle_ptr, + const cuopt::mps_parser::mps_data_model_t& data_model) { cuopt_expects(handle_ptr != nullptr, error_type_t::ValidationError, @@ -1635,7 +1636,11 @@ cuopt::linear_programming::optimization_problem_t mps_data_model_to_op data_model.get_variable_types().cbegin(), data_model.get_variable_types().cend(), enum_variable_types.begin(), - [](const auto val) -> var_t { return val == 'I' ? var_t::INTEGER : var_t::CONTINUOUS; }); + [](const auto val) -> var_t { + if (val == 'I' || val == 'B') return var_t::INTEGER; + if (val == 'S') return var_t::SEMI_CONTINUOUS; + return var_t::CONTINUOUS; + }); op_problem.set_variable_types(enum_variable_types.data(), enum_variable_types.size()); } @@ -1814,7 +1819,7 @@ std::unique_ptr> solve_lp( \ template optimization_problem_t mps_data_model_to_optimization_problem( \ raft::handle_t const* handle_ptr, \ - const cuopt::mps_parser::mps_data_model_t& data_model); \ + const cuopt::mps_parser::mps_data_model_t& data_model); \ template void set_pdlp_solver_mode(pdlp_solver_settings_t& settings); #if MIP_INSTANTIATE_FLOAT diff --git a/cpp/src/pdlp/utilities/problem_checking.cu b/cpp/src/pdlp/utilities/problem_checking.cu index b10850de27..88dd356cfb 100644 --- a/cpp/src/pdlp/utilities/problem_checking.cu +++ b/cpp/src/pdlp/utilities/problem_checking.cu @@ -217,6 +217,38 @@ void problem_checking_t::check_problem_representation( op_problem.get_objective_coefficients().size(), op_problem.get_variable_upper_bounds().size()); } + if (!op_problem.get_variable_types().is_empty()) { + cuopt_expects( + op_problem.get_variable_types().size() == op_problem.get_objective_coefficients().size(), + error_type_t::ValidationError, + "Sizes for vectors related to the variables are not the same. The objective " + "vector has size %zu and the variable types vector has size %zu.", + op_problem.get_objective_coefficients().size(), + op_problem.get_variable_types().size()); + + auto var_types = op_problem.get_variable_types_host(); + auto var_lb = op_problem.get_variable_lower_bounds_host(); + auto var_ub = op_problem.get_variable_upper_bounds_host(); + for (i_t i = 0; i < static_cast(var_types.size()); ++i) { + if (var_types[i] != var_t::SEMI_CONTINUOUS) { continue; } + cuopt_expects(var_lb[i] > f_t(0), + error_type_t::ValidationError, + "Semi-continuous variable must have a strictly positive lower bound, but has " + "lower bound %g.", + static_cast(var_lb[i])); + cuopt_expects(var_ub[i] > f_t(0), + error_type_t::ValidationError, + "Semi-continuous variable must have a strictly positive upper bound, but has " + "upper bound %g.", + static_cast(var_ub[i])); + cuopt_expects(var_lb[i] < var_ub[i], + error_type_t::ValidationError, + "Semi-continuous variable must satisfy lower bound < upper bound, but has " + "bounds [%g, %g].", + static_cast(var_lb[i]), + static_cast(var_ub[i])); + } + } // Check constraints sizes cuopt_expects( diff --git a/cpp/tests/mip/CMakeLists.txt b/cpp/tests/mip/CMakeLists.txt index f2cf53ff6c..13eec1772a 100644 --- a/cpp/tests/mip/CMakeLists.txt +++ b/cpp/tests/mip/CMakeLists.txt @@ -8,6 +8,9 @@ ConfigureTest(MIP_TEST ${CMAKE_CURRENT_SOURCE_DIR}/miplib_test.cu ) +ConfigureTest(SEMI_CONTINUOUS_TEST + ${CMAKE_CURRENT_SOURCE_DIR}/semi_continuous_test.cu +) ConfigureTest(PROBLEM_TEST ${CMAKE_CURRENT_SOURCE_DIR}/problem_test.cu ) diff --git a/cpp/tests/mip/semi_continuous_test.cu b/cpp/tests/mip/semi_continuous_test.cu new file mode 100644 index 0000000000..192ebf18cb --- /dev/null +++ b/cpp/tests/mip/semi_continuous_test.cu @@ -0,0 +1,121 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2024-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../linear_programming/utilities/pdlp_test_utilities.cuh" +#include "cuopt/linear_programming/mip/solver_settings.hpp" + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +namespace cuopt::linear_programming::test { + +struct sc_result_t { + std::string file; + double objective; + double sc_value; +}; + +optimization_problem_t make_sc_problem(raft::handle_t const* handle, + double sc_lb, + double sc_ub) +{ + optimization_problem_t problem(handle); + + const std::vector coefficients = {1.0, 1.0}; + const std::vector indices = {0, 1}; + const std::vector offsets = {0, 2}; + const std::vector row_lower = {1.0}; + const std::vector row_upper = {1.0}; + const std::vector obj = {1.0, 0.0}; + const std::vector var_lower = {sc_lb, 0.0}; + const std::vector var_upper = {sc_ub, 1.0}; + const std::vector var_types = {var_t::SEMI_CONTINUOUS, var_t::CONTINUOUS}; + + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + problem.set_constraint_lower_bounds(row_lower.data(), row_lower.size()); + problem.set_constraint_upper_bounds(row_upper.data(), row_upper.size()); + problem.set_objective_coefficients(obj.data(), obj.size()); + problem.set_variable_lower_bounds(var_lower.data(), var_lower.size()); + problem.set_variable_upper_bounds(var_upper.data(), var_upper.size()); + problem.set_variable_types(var_types.data(), var_types.size()); + + return problem; +} + +TEST(mip_solve, semi_continuous_regressions) +{ + const raft::handle_t handle_{}; + mip_solver_settings_t settings; + settings.time_limit = 10.; + + const std::vector valid_test_instances = { + {"mip/sc_standard.mps", 8., 0.}, + {"mip/sc_no_ub.mps", 8., 0.}, + {"mip/sc_lb_zero.mps", 8., 0.}, + {"mip/sc_inferred_ub.mps", -4., 4.}, + }; + + for (const auto& test_instance : valid_test_instances) { + auto path = make_path_absolute(test_instance.file); + auto problem = cuopt::mps_parser::parse_mps(path, false); + auto solution = solve_mip(&handle_, problem, settings); + + EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal) + << test_instance.file; + ASSERT_EQ(solution.get_solution().size(), static_cast(problem.get_n_variables())) + << test_instance.file; + + auto host_solution = + cuopt::host_copy(solution.get_solution(), solution.get_solution().stream()); + EXPECT_NEAR(solution.get_objective_value(), test_instance.objective, 1e-6) + << test_instance.file; + EXPECT_NEAR(host_solution[0], test_instance.sc_value, 1e-6) << test_instance.file; + } +} + +TEST(mip_solve, semi_continuous_invalid_bounds_rejected) +{ + const raft::handle_t handle_{}; + mip_solver_settings_t settings; + settings.time_limit = 10.; + + const std::vector> invalid_bounds = { + {-3.0, 5.0}, + {-5.0, -1.0}, + {-4.0, 0.0}, + {5.0, 5.0}, + {6.0, 5.0}, + }; + + for (const auto& [lb, ub] : invalid_bounds) { + SCOPED_TRACE(::testing::Message() << "bounds=[" << lb << ", " << ub << "]"); + auto problem = make_sc_problem(&handle_, lb, ub); + + auto solution = solve_mip(problem, settings); + const auto& error = solution.get_error_status(); + EXPECT_EQ(error.get_error_type(), cuopt::error_type_t::ValidationError); + EXPECT_NE(std::string(error.what()).find("Semi-continuous variable"), std::string::npos); + } +} + +} // namespace cuopt::linear_programming::test diff --git a/datasets/mip/sc_inferred_ub.mps b/datasets/mip/sc_inferred_ub.mps new file mode 100644 index 0000000000..a9b4415cf5 --- /dev/null +++ b/datasets/mip/sc_inferred_ub.mps @@ -0,0 +1,15 @@ +NAME sc_inferred_ub +ROWS + N OBJ + L cap +COLUMNS + x OBJ -1 + x cap 1 + y cap 1 +RHS + RHS1 cap 4 +BOUNDS + LO BND1 x 2 + SC BND1 x 1e+30 + UP BND1 y 10 +ENDATA diff --git a/datasets/mip/sc_lb_zero.mps b/datasets/mip/sc_lb_zero.mps new file mode 100644 index 0000000000..a39f434ed9 --- /dev/null +++ b/datasets/mip/sc_lb_zero.mps @@ -0,0 +1,16 @@ +NAME sc_lb_zero +ROWS + N OBJ + G cover +COLUMNS + x OBJ 3 + x cover 1 + y OBJ 2 + y cover 1 +RHS + RHS1 cover 4 +BOUNDS + SC BND1 x 10 + LO BND1 y -10 + UP BND1 y 10 +ENDATA diff --git a/datasets/mip/sc_no_ub.mps b/datasets/mip/sc_no_ub.mps new file mode 100644 index 0000000000..565a475d11 --- /dev/null +++ b/datasets/mip/sc_no_ub.mps @@ -0,0 +1,17 @@ +NAME sc_no_ub +ROWS + N OBJ + G cover +COLUMNS + x OBJ 3 + x cover 1 + y OBJ 2 + y cover 1 +RHS + RHS1 cover 4 +BOUNDS + LO BND1 x 2 + SC BND1 x 1e+30 + LO BND1 y -10 + UP BND1 y 10 +ENDATA diff --git a/datasets/mip/sc_standard.mps b/datasets/mip/sc_standard.mps new file mode 100644 index 0000000000..8eaa3691fa --- /dev/null +++ b/datasets/mip/sc_standard.mps @@ -0,0 +1,17 @@ +NAME sc_standard +ROWS + N OBJ + G cover +COLUMNS + x OBJ 3 + x cover 1 + y OBJ 2 + y cover 1 +RHS + RHS1 cover 4 +BOUNDS + LO BND1 x 2 + SC BND1 x 10 + LO BND1 y -10 + UP BND1 y 10 +ENDATA