From abf551b569e0405842d1c28d3cdb9500ff867451 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Mon, 15 Jan 2024 16:22:59 +1100 Subject: [PATCH] SOCP: move QP terms into aux SOC #229 #192 When any cones found, move QP objective terms into an aux SOC, in the easy case Specially for Mosek --- include/mp/flat/expr_quadratic.h | 7 ++ include/mp/flat/obj_std.h | 8 ++- include/mp/flat/redef/conic/cones.h | 64 ++++++++++++++++--- .../categorized/fast/conic/modellist.json | 12 +++- .../fast/conic/socp_12__issue229_2.mod | 7 ++ test/end2end/scripts/python/Solver.py | 1 + 6 files changed, 87 insertions(+), 12 deletions(-) create mode 100644 test/end2end/cases/categorized/fast/conic/socp_12__issue229_2.mod diff --git a/include/mp/flat/expr_quadratic.h b/include/mp/flat/expr_quadratic.h index da7abd882..014fc0718 100644 --- a/include/mp/flat/expr_quadratic.h +++ b/include/mp/flat/expr_quadratic.h @@ -107,6 +107,13 @@ class QuadTerms { /// Sort and eliminate duplicates void sort_terms(); + /// Clear + void clear() { + coefs_.clear(); + vars1_.clear(); + vars2_.clear(); + } + /// Test equality bool equals(const QuadTerms& qt) const { return *this == qt; diff --git a/include/mp/flat/obj_std.h b/include/mp/flat/obj_std.h index 6f83e9d42..80286efee 100644 --- a/include/mp/flat/obj_std.h +++ b/include/mp/flat/obj_std.h @@ -26,8 +26,10 @@ class LinearObjective { name_(std::move(nm)){ } /// Get sense obj::Type obj_sense() const { return sense_; } - /// Get lin terms + /// Get lin terms, const const LinTerms& GetLinTerms() const { return lt_; } + /// Get lin terms + LinTerms& GetLinTerms() { return lt_; } /// Get N terms int num_terms() const { assert(check()); return (int)lt_.size(); } /// Validate @@ -55,8 +57,10 @@ class QuadraticObjective : public LinearObjective { QuadraticObjective(LinearObjective&& lc, QuadTerms&& qt) : LinearObjective(std::move(lc)), qt_(std::move(qt)) { sort_qp_terms(); } - /// Get QP terms + /// Get QP terms, const const QuadTerms& GetQPTerms() const { return qt_; } + /// Get QP terms + QuadTerms& GetQPTerms() { return qt_; } /// Sort QP terms void sort_qp_terms() { diff --git a/include/mp/flat/redef/conic/cones.h b/include/mp/flat/redef/conic/cones.h index 9801e6937..38414f8b4 100644 --- a/include/mp/flat/redef/conic/cones.h +++ b/include/mp/flat/redef/conic/cones.h @@ -85,6 +85,7 @@ class ConicConverter : public MCKeeper { RunQConesFromNonQC(); if (MC().IfPassExpCones()) RunExpCones(); + TryReformulateQPObjective(); if (IfNeedSOCP2QC()) SetupSOCP2QC(); WarnOnMix(); @@ -119,11 +120,62 @@ class ConicConverter : public MCKeeper { } } + /// Ideally also for exp/log terms + void TryReformulateQPObjective() { + if (HasAnyCones() + && MC().HasQPObjective() && 1==MC().num_objs()) { + // See if the QP objective is SOCP-easy + auto& obj = MC().get_objectives().at(0); + auto& qp = obj.GetQPTerms(); + auto objsns = (obj::MAX==obj.obj_sense()) ? 1.0 : -1.0; + for (auto i=qp.size(); i--; ) { + if (qp.coef(i) * objsns > 0.0 + || qp.var1(i) != qp.var2(i)) // not x1[i]==x2[i] + return; + } + // Set up AutoLink + auto obj_src = // source value node for this obj + MC().GetValuePresolver().GetSourceNodes().GetObjValues()().Select(0); + pre::AutoLinkScope auto_link_scope{ MC(), obj_src }; + // Aux vars + int z = (int)MC().AddVar(0.0, MC().Infty()); + int z1 = (int)MC().MakeFixedVar(1.0); + // Add cone + std::vector c = {{1.0, 0.5}}; + c.insert(c.end(), qp.coefs().begin(), qp.coefs().end()); + if (objsns > 0) // negate coefs if maximizing + for (size_t i=2; i x = {{z, z1}}; + x.insert(x.end(), qp.vars1().begin(), qp.vars1().end()); + MC().AddConstraint( + RotatedQuadraticConeConstraint( + std::move(x), std::move(c))); + // Add linear term, remove QP terms + qp.clear(); + obj.GetLinTerms().add_term(-objsns, z); + } + } + + /// Any cones at all + bool HasAnyCones() const { + return MC().GetNumberOfAddable((QuadraticConeConstraint*)0)>0 || + MC().GetNumberOfAddable((RotatedQuadraticConeConstraint*)0)>0 || + HasAnyNonSOCPCones(); + } + + /// any non-SOCP cones? + bool HasAnyNonSOCPCones() const { + return MC().GetNumberOfAddable((ExponentialConeConstraint*)0)>0 || + MC().GetNumberOfAddable((PowerConeConstraint*)0)>0 || + MC().GetNumberOfAddable((GeometricConeConstraint*)0)>0; + } + bool IfNeedSOCP2QC() { return MC().SOCP2QCMode() >= 2 // compulsory || (1 == MC().SOCP2QCMode() - && 0 == MC().NumExpConesRecognized() + && !HasAnyNonSOCPCones() // Might also have SOCP from sqrt(), abs(). // Some QC -> SOCP but not all, even if 0 succeeded. && (MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded() @@ -137,14 +189,10 @@ class ConicConverter : public MCKeeper { void WarnOnMix() { if ( !MC().ModelAPICanMixConicQCAndQC()) { // cannot mix - if ((MC().NumExpConesRecognized() // exp cones + if ((HasAnyCones() // exp cones && (MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded() || MC().HasQPObjective())) // and quadratics left in - || // Some QC -> SOCP but not all - ((MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded() - || MC().HasQPObjective()) // or a quadratic obj - && MC().NumQC2SOCPSucceeded() // Warn only if some succeeded - && !MC().IfConvertSOCP2QC()) // and not decided to convert + && !MC().IfConvertSOCP2QC() // and not decided to convert ) { MC().AddWarning("Mix QC+cones", "Not all quadratic constraints could " @@ -153,7 +201,7 @@ class ConicConverter : public MCKeeper { "additionally, further convertion back to QC\n" "not desired (option cvt:socp2qc) or other cone types present;\n" "solver might not accept the model.\n" - "Try to express all SOCP cones in standard forms,\n" + "Try to express all cones in standard forms,\n" "not in the objective.\n" "See mp.ampl.com/model-guide.html."); } diff --git a/test/end2end/cases/categorized/fast/conic/modellist.json b/test/end2end/cases/categorized/fast/conic/modellist.json index 9d32d38db..e4d338ffe 100644 --- a/test/end2end/cases/categorized/fast/conic/modellist.json +++ b/test/end2end/cases/categorized/fast/conic/modellist.json @@ -28,13 +28,13 @@ "name" : "TotalVariation2D", "objective" : 9.569949819774543, "files": ["totalvariation2d.mod", "totalvar2d_4e_4.dat"], - "tags" : ["socp", "socp_hard_to_recognize"] + "tags" : ["socp"] }, { "name" : "TotalVariation2D_Quadr", "objective" : 9.569949819774543, "files": ["totalvariation2d_quad.mod", "totalvar2d_4e_4.dat"], - "tags" : ["socp", "socp_hard_to_recognize"] + "tags" : ["socp"] }, { "name" : "socp_06_rsoc_constants", @@ -77,6 +77,14 @@ "solve_result_num": 0 } }, + { + "name" : "socp_12__issue229_2", + "tags" : ["socp"], + "objective" : -1607.1428, + "values": { + "solve_result_num": 0 + } + }, { "name" : "expcones_01__plain", "objective" : 0.7821882953, diff --git a/test/end2end/cases/categorized/fast/conic/socp_12__issue229_2.mod b/test/end2end/cases/categorized/fast/conic/socp_12__issue229_2.mod new file mode 100644 index 000000000..cd5b1eb10 --- /dev/null +++ b/test/end2end/cases/categorized/fast/conic/socp_12__issue229_2.mod @@ -0,0 +1,7 @@ +var x {1..7} >= 0; +var y >= 0; + +maximize obj: sum {j in 1..7} -x[j]^2 - y^2; + +subj to conQ: sum {j in 1..7} x[j]^2 <= y^2; +subj to conL: sum {j in 1..7} x[j] = 75; diff --git a/test/end2end/scripts/python/Solver.py b/test/end2end/scripts/python/Solver.py index 8d8f2ebbf..379f6c207 100644 --- a/test/end2end/scripts/python/Solver.py +++ b/test/end2end/scripts/python/Solver.py @@ -729,6 +729,7 @@ def __init__(self, exeName, timeout=None, nthreads=None, ModelTags.continuous, ModelTags.integer, ModelTags.binary, ModelTags.sos, ModelTags.quadratic, ModelTags.quadratic_obj, + ModelTags.socp, ModelTags.return_mipgap, ModelTags.multiobj,