From 73805c4e18547d78463a46d672e9e41d75908eb6 Mon Sep 17 00:00:00 2001 From: Gleb Belov Date: Tue, 21 Mar 2023 16:47:59 +1100 Subject: [PATCH] Rotated SOC from quadratics: more constants #192 Allowing constants in place of variables in RSOC-equivalent quadratics --- include/mp/flat/redef/conic/cones.h | 84 ++++++++++++++----- include/mp/flat/redef/conic/qcones2qc.h | 4 +- .../categorized/fast/conic/modellist.json | 5 ++ .../fast/conic/socp_06_rsoc_constants.mod | 21 +++++ 4 files changed, 91 insertions(+), 23 deletions(-) create mode 100644 test/end2end/cases/categorized/fast/conic/socp_06_rsoc_constants.mod diff --git a/include/mp/flat/redef/conic/cones.h b/include/mp/flat/redef/conic/cones.h index b405d8314..1ff9e965e 100644 --- a/include/mp/flat/redef/conic/cones.h +++ b/include/mp/flat/redef/conic/cones.h @@ -97,8 +97,9 @@ class ConicConverter : public MCKeeper { }; -/** Converter for a single quadratic inequality constraint. - * Generic, receives only body (assumes rhs=0) and sense. +/** Converter for a single linear / quadratic + * inequality constraint. + * Generic, receives only body and sense. */ template class Convert1QC : public MCKeeper { @@ -120,7 +121,6 @@ class Convert1QC : public MCKeeper { /// DoRun. Body: quadratic. /// - /// Currently considering rhs=0. /// Accept non-(+-1) coefficients. /// /// @param body: quadratic constraint body. @@ -131,7 +131,7 @@ class Convert1QC : public MCKeeper { bool DoRun(const QuadAndLinTerms& body, int sens, double rhs) { assert((sens==1 || sens==-1) && "sens 1 or -1 only"); - if (body.GetLinTerms().empty()) { + if (body.GetLinTerms().size()<=1) { QPTermsTraits qptt; if (Characterize(body.GetQPTerms(), qptt)) return ProceedQCWithTraits(body, sens, rhs, qptt); @@ -165,21 +165,44 @@ class Convert1QC : public MCKeeper { return true; } - /// Proceed with a quadratic constraint provided traits - /// @return true iff cone added + /// Proceed with a quadratic constraint provided traits. + /// + /// Cases of StdSOC: + /// 1. x^2 >= ||y||^2 + /// 2. c^2 >= ||y||^2 + /// 3. x^2 >= ||y||^2 + c^2 + /// 4. 0 >= ||y||^2 + c^2 -- fail? + /// + /// Cases of RotSOC: + /// 1. 2xy >= ||z||^2 [ + c^2 ] + /// 2. 2 c1 y >= ||z||^2 [ + c2^2 ] + /// + /// @return true iff cone added. bool ProceedQCWithTraits(const QuadAndLinTerms& body, int sens, double rhs, const QPTermsTraits& qptt) { - assert(body.GetLinTerms().empty()); + const auto& lint = body.GetLinTerms(); + assert(lint.size() <= 1); + assert(qptt.nDiffVars <= 1); const auto& qpterms = body.GetQPTerms(); - if (qptt.nDiffVars) { // Rotated cone? - if (0.0 == std::fabs(rhs)) { // rhs==0 - if ((sens>0 && qptt.coef12>0.0 && + if (qptt.nDiffVars || lint.size()) { // Rotated cone? + const auto fLT1 = lint.size()==1; // Constant in the lhs? + const auto fND1 = qptt.nDiffVars==1; + if (fLT1 + fND1 == 1 && + rhs * sens >= 0.0) { // ... >= ... + rhs + if (fLT1 && // Check lin term + (MC().lb( lint.var(0) ) < 0.0 || + lint.coef(0)*sens < 0.0)) + return false; + if (fND1 && qptt.coef12*sens<0.0) + return false; + if ((sens>0 && qptt.nSamePos==0 && qptt.nSameNeg) || - (sens<0 && qptt.coef12<0.0 && + (sens<0 && qptt.nSameNeg==0 && qptt.nSamePos)) - return AddRotatedQC(qpterms, qptt.iDiffVars); + return AddRotatedQC(qpterms, lint, rhs, + qptt.iDiffVars); } } else if (0.0 >= rhs*sens) { // Standard cone? if (0.0 == std::fabs(rhs)) { @@ -200,9 +223,10 @@ class Convert1QC : public MCKeeper { return false; } - /// DoRun. Body: linear. + /// DoRun. Body: linear, i.e., var1 (or const) >= var2 + /// (then checking if var2 := ||.||). /// - /// Currently considering rhs=0 or rhs>0. + /// Considering const>=0. /// Accept non-(+-1) coefficients. /// /// @param body: linear constraint body. @@ -236,8 +260,8 @@ class Convert1QC : public MCKeeper { // it's -k*sqrt(...) <= rhs(>0), k>0, // which is always true TODO if (auto rhs_args = CheckNorm2(lint.var(0))) { - return ContinueStdSOC(1.0, - MC().MakeFixedVar(std::fabs(rhs)), + return ContinueStdSOC(std::fabs(rhs), + MC().MakeFixedVar( 1.0 ), lint.coef(0), rhs_args); } } @@ -417,22 +441,38 @@ class Convert1QC : public MCKeeper { /// Add rotated cone from pure-quadratic constraint - bool AddRotatedQC(const QuadTerms& qpterms, int iDiffVars) { - std::vector x(qpterms.size()+1); - std::vector c(qpterms.size()+1); + bool AddRotatedQC(const QuadTerms& qpterms, + const LinTerms& lint, + double rhs, int iDiffVars) { + assert(lint.size() + (iDiffVars>=0) == 1); + const auto rhsNon0 = (std::fabs(rhs)!=0.0); + std::vector x(qpterms.size()+1+lint.size()+rhsNon0); + std::vector c(qpterms.size()+1+lint.size()+rhsNon0); size_t iPush=1; + if (lint.size()) { + x[0] = lint.var(0); + x[1] = MC().MakeFixedVar( 1.0 ); + c[0] = std::fabs(lint.coef(0)) / 2.0; + c[1] = 1.0; // it's 2xy >= z^2 + ... + } for (int i=0; i<(int)qpterms.size(); ++i) { if (i==iDiffVars) { x[0] = qpterms.var1(i); x[1] = qpterms.var2(i); c[0] = std::fabs(qpterms.coef(i)); - c[1] = 0.5; // it's 2xy >= z^2 + ... + c[1] = 0.5; // it's 2xy >= z^2 + ... } else { ++iPush; x.at(iPush) = qpterms.var1(i); c.at(iPush) = std::sqrt(std::fabs(qpterms.coef(i))); } } + if (rhsNon0) { // Add ... >= ... + |rhs| + ++iPush; + x.at(iPush) = MC().MakeFixedVar( 1.0 ); + c.at(iPush) = std::sqrt( std::fabs(rhs) ); + } + assert(iPush == x.size()-1); MC().AddConstraint( RotatedQuadraticConeConstraint( std::move(x), std::move(c))); @@ -447,8 +487,8 @@ class Convert1QC : public MCKeeper { std::vector x(qpterms.size() + fNewVar); std::vector c(qpterms.size() + fNewVar); if (fNewVar) { - c[0] = 1.0; - x[0] = MC().MakeFixedVar(std::sqrt(std::fabs(rhs))); + c[0] = std::sqrt(std::fabs(rhs)); + x[0] = MC().MakeFixedVar( 1.0 ); } size_t iPush=0; for (int i=0; i<(int)qpterms.size(); ++i) { diff --git a/include/mp/flat/redef/conic/qcones2qc.h b/include/mp/flat/redef/conic/qcones2qc.h index fa1f7065d..3ef1540e4 100644 --- a/include/mp/flat/redef/conic/qcones2qc.h +++ b/include/mp/flat/redef/conic/qcones2qc.h @@ -27,7 +27,9 @@ class QConeConverter : 0!=GetMC().IfPassQuadCon(); } - /// Convert to (c[0]*x[0])^2 >= sum(i)((c[i]*x[i])^2). + /// Convert to + /// (c[0]*x[0])^2 >= sum(i)((c[i]*x[i])^2) + /// with x[0]>=0. void Convert(const ItemType& ac, int ) { const auto& x = ac.GetArguments(); auto c = ac.GetParameters(); diff --git a/test/end2end/cases/categorized/fast/conic/modellist.json b/test/end2end/cases/categorized/fast/conic/modellist.json index d6b7b3b40..aa5bb40cd 100644 --- a/test/end2end/cases/categorized/fast/conic/modellist.json +++ b/test/end2end/cases/categorized/fast/conic/modellist.json @@ -35,5 +35,10 @@ "objective" : 9.569949819774543, "files": ["totalvariation2d_quad.mod", "totalvar2d_4e_4.dat"], "tags" : ["socp", "socp_hard_to_recognize"] + }, + { + "name" : "socp_06_rsoc_constants", + "objective" : 7.993662815, + "tags" : ["socp"] } ] diff --git a/test/end2end/cases/categorized/fast/conic/socp_06_rsoc_constants.mod b/test/end2end/cases/categorized/fast/conic/socp_06_rsoc_constants.mod new file mode 100644 index 000000000..b6d9564f0 --- /dev/null +++ b/test/end2end/cases/categorized/fast/conic/socp_06_rsoc_constants.mod @@ -0,0 +1,21 @@ +# socp_04_2cones.mod + +var x {1..3} >=0, <=20; + +minimize Obj: + 2*x[1] + 3*x[2] - 0.5*x[3]; + +s.t. StdCone: + sqrt(19.03*x[1]^2 + x[3]^2) <= sqrt(3.72)*x[2]; + +s.t. RotatedCone1: + -25*x[1]*x[2] <= -16*x[3]^2 - 7; + +s.t. RotatedCone2: + -25*x[2] <= -16*x[3]^2; + +s.t. RotatedCone3: + -25*x[2] <= -16*x[3]^2 - 7; + +s.t. LinCon: + sum {i in 1..3} x[i] == 5;