diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 4e9f827..16d3f04 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -87,6 +87,7 @@ extern "C" double eps_optimal_relative; double eps_feasible_relative; double eps_feas_polish_relative; + double eps_infeasible_relative; double time_sec_limit; int iteration_limit; } termination_criteria_t; diff --git a/python/README.md b/python/README.md index be188dd..0a84ed1 100644 --- a/python/README.md +++ b/python/README.md @@ -163,6 +163,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. | | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. | | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. | +| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-14` | Relative tolerance on the ray certificate for declaring primal/dual infeasibility. | | `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. | | `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations of L∞ Ruiz scaling. Improves numerical conditioning. | | `UsePCAlpha` | `has_pock_chambolle_alpha` | bool | `True` | Whether to use the Pock–Chambolle α step size adjustment. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index d6ca197..fb0cb2f 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -41,6 +41,7 @@ # tolerances "OptimalityTol": "eps_optimal_relative", "FeasibilityTol": "eps_feasible_relative", + "InfeasibleTol": "eps_infeasible_relative", # scaling / step size "GeoMeanIters": "geometric_mean_iterations", "RuizIters": "l_inf_ruiz_iterations", diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index fdc23b3..dcb7f6c 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -53,6 +53,7 @@ { "eps_optimal_relative", "eps_feasible_relative", + "eps_infeasible_relative", "time_sec_limit", "pock_chambolle_alpha", "artificial_restart_threshold", @@ -70,6 +71,7 @@ { "eps_optimal_relative", "eps_feasible_relative", + "eps_infeasible_relative", "eps_feas_polish_relative", "sv_tol", "infinite_bound", diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index aee848c..ffc6816 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -273,6 +273,7 @@ static py::dict get_default_params_py() // tolerances d["eps_optimal_relative"] = p.termination_criteria.eps_optimal_relative; d["eps_feasible_relative"] = p.termination_criteria.eps_feasible_relative; + d["eps_infeasible_relative"] = p.termination_criteria.eps_infeasible_relative; // limits d["time_sec_limit"] = p.termination_criteria.time_sec_limit; @@ -384,6 +385,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p // tolerances getf("eps_optimal_relative", p->termination_criteria.eps_optimal_relative); getf("eps_feasible_relative", p->termination_criteria.eps_feasible_relative); + getf("eps_infeasible_relative", p->termination_criteria.eps_infeasible_relative); // limits getf("time_sec_limit", p->termination_criteria.time_sec_limit); @@ -435,6 +437,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p throw std::invalid_argument("eps_optimal_relative must be positive."); if (p->termination_criteria.eps_feasible_relative <= 0.0) throw std::invalid_argument("eps_feasible_relative must be positive."); + if (p->termination_criteria.eps_infeasible_relative <= 0.0) + throw std::invalid_argument("eps_infeasible_relative must be positive."); if (p->termination_criteria.eps_feas_polish_relative <= 0.0) throw std::invalid_argument("eps_feas_polish_relative must be positive."); if (p->sv_tol <= 0.0) diff --git a/src/solver.cu b/src/solver.cu index 3e12d57..636c7b5 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -203,6 +203,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t * compute_fixed_point_error(state); compute_residual(state, params->optimality_norm); + compute_infeasibility_information(state); state->inner_count += params->termination_evaluation_frequency; state->total_count += params->termination_evaluation_frequency; diff --git a/src/utils.cu b/src/utils.cu index b4bad93..7ff69ce 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -254,6 +254,16 @@ void check_termination_criteria(pdhg_solver_state_t *solver_state, const termina solver_state->termination_reason = TERMINATION_REASON_OPTIMAL; return; } + if (primal_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) + { + solver_state->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE; + return; + } + if (dual_infeasibility_criteria_met(solver_state, criteria->eps_infeasible_relative)) + { + solver_state->termination_reason = TERMINATION_REASON_DUAL_INFEASIBLE; + return; + } if (solver_state->total_count >= criteria->iteration_limit) { solver_state->termination_reason = TERMINATION_REASON_ITERATION_LIMIT; @@ -318,6 +328,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->termination_criteria.time_sec_limit = 3600.0; params->termination_criteria.iteration_limit = INT32_MAX; params->termination_criteria.eps_feas_polish_relative = 1e-6; + params->termination_criteria.eps_infeasible_relative = 1e-14; params->restart_params.artificial_restart_threshold = 0.36; params->restart_params.sufficient_reduction_for_restart = 0.2; @@ -625,6 +636,9 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro PRINT_DIFF_DBL("eps_feas_polish_relative", params->termination_criteria.eps_feas_polish_relative, default_params.termination_criteria.eps_feas_polish_relative); + PRINT_DIFF_DBL("eps_infeasible_relative", + params->termination_criteria.eps_infeasible_relative, + default_params.termination_criteria.eps_infeasible_relative); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); PRINT_DIFF_DBL("infinite_bound", params->infinite_bound, default_params.infinite_bound); diff --git a/test/test_infeasible_unbounded.py b/test/test_infeasible_unbounded.py index 8fffaab..60eefe9 100644 --- a/test/test_infeasible_unbounded.py +++ b/test/test_infeasible_unbounded.py @@ -35,13 +35,15 @@ def test_infeasible_lp(base_lp_data, atol): model.setConstraintUpperBound(u) # turn off output model.setParams(OutputFlag=False, Presolve=False) + # set infeasibility tolerance + model.setParams(InfeasibleTol=1e-6) # optimize - #model.optimize() + model.optimize() # check status - #assert hasattr(model, "Status"), "Model.Status not exposed." - #assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" + assert hasattr(model, "Status"), "Model.Status not exposed." + assert model.Status == PDLP.PRIMAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check dual ray - #assert model.DualRayObj > atol, f"DualRayObj should be positive for dual infeasible, got {model.DualRayObj}" + assert model.DualRayObj > atol, f"DualRayObj should be positive for primal infeasible, got {model.DualRayObj}" def test_unbounded_lp(base_lp_data, atol): @@ -64,10 +66,12 @@ def test_unbounded_lp(base_lp_data, atol): model.setVariableLowerBound(lb) # turn off output model.setParams(OutputFlag=False, Presolve=False) + # set infeasibility tolerance + model.setParams(InfeasibleTol=1e-6) # optimize - #model.optimize() + model.optimize() # check status - #assert hasattr(model, "Status"), "Model.Status not exposed." - #assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" + assert hasattr(model, "Status"), "Model.Status not exposed." + assert model.Status == PDLP.DUAL_INFEASIBLE, f"Unexpected termination status: {model.Status}" # check primal ray - #assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}" \ No newline at end of file + assert model.PrimalRayLinObj < -atol, f"PrimalRayLinObj should be negative for dual infeasible, got {model.PrimalRayLinObj}" \ No newline at end of file