diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 5ca948690..824849ce9 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -190,15 +190,16 @@ using ComplementarityQuadratic = ComplementarityConstraint; /// Quadratic cone DEF_STATIC_CONSTR_WITH_PRM( QuadraticConeConstraint, VarArray, DblParamArray, - "Quadratic cone x1 >= sqrt(x2^2 + ...)) with factors " - "applied to the arguments, aka MOSEK 10 affine cones"); + "Quadratic cone p1*x1 >= sqrt((p2*x2)^2 + ...))," + " with factors p1..pn"); /// Rotated quadratic cone DEF_STATIC_CONSTR_WITH_PRM( RotatedQuadraticConeConstraint, VarArray, DblParamArray, - "Rotated quadratic cone x1*x2 >= sqrt(x3^2 + ...)) with factors " - "applied to the arguments, aka MOSEK 10 affine cones"); + "Rotated quadratic cone p1*x1*p2*x2 >= sqrt((p3*x3)^2 + ...))," + " x1, x2 >= 0, with factors p1..pn"); /// Exponential cone DEF_STATIC_CONSTR_WITH_PRM( ExponentialConeConstraint, VarArray3, DblParamArray3, - "Exponential cone with factors"); + "Exponential cone p1*x1 >= p2*x2*exp(p3*x3 / (p2*x2))," + " x1, x2 >= 0, with factors p1..p3"); /// Power cone DEF_STATIC_CONSTR_WITH_PRM( PowerConeConstraint, VarArray, DblParamArray, "Power cone with factors "); diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h index 7cfc1dd2a..0bb84cbea 100644 --- a/include/mp/flat/converter.h +++ b/include/mp/flat/converter.h @@ -826,6 +826,12 @@ class FlatConverter : (int)GetConstraintAcceptance((RotatedQuadraticConeConstraint*)nullptr)); } + /// Whether the ModelAPI accepts exp cones + int ModelAPIAcceptsExponentialCones() { + return + (int)GetConstraintAcceptance((ExponentialConeConstraint*)nullptr); + } + private: struct Options { @@ -837,6 +843,7 @@ class FlatConverter : int passQuadObj_ = ModelAPIAcceptsQuadObj(); int passQuadCon_ = ModelAPIAcceptsQC(); int passSOCPCones_ = 0; + int passExpCones_ = 0; int relax_ = 0; }; @@ -893,6 +900,13 @@ class FlatConverter : "0*/1: Multiply out and pass quadratic constraint terms to the solver, " "vs. linear approximation.", options_.passQuadCon_, 0, 1); + if (ModelAPIAcceptsExponentialCones()) + GetEnv().AddOption("cvt:expcones expcones", + ModelAPIAcceptsExponentialCones()>1 ? + "0/1*: Recognize exponential cones." : + "0*/1: Recognize exponential cones.", + options_.passExpCones_, 0, 1); + options_.passExpCones_ = ModelAPIAcceptsExponentialCones()>1; if (ModelAPIAcceptsQuadraticCones()) GetEnv().AddOption("cvt:socp passsocp socp", ModelAPIAcceptsQuadraticCones()>1 ? @@ -939,6 +953,9 @@ class FlatConverter : /// Whether we pass SOCP cones bool IfPassSOCPCones() const { return options_.passSOCPCones_; } + /// Whether we pass exp cones + bool IfPassExpCones() const { return options_.passExpCones_; } + public: /// Typedef ModelAPIType. For tests diff --git a/include/mp/flat/redef/conic/cones.h b/include/mp/flat/redef/conic/cones.h index 40d53e67e..3ac467718 100644 --- a/include/mp/flat/redef/conic/cones.h +++ b/include/mp/flat/redef/conic/cones.h @@ -12,7 +12,6 @@ #include #include -#include "mp/flat/constr_std.h" #include "mp/flat/constr_keeper.h" #include "mp/valcvt-link.h" @@ -20,9 +19,22 @@ namespace mp { /// Abstract converter for a single constraint. -template +/// +/// @param MCType: ModelConverter, e.g., FlatConverter. +/// @param Con: the source constraint type. +/// @param ConvertBase: another class +/// providing conversion methods DoRun +/// for quadratic and linear constraints. +template class CvtBase> class Convert1 { }; +/// Declare a base converter into quadratic cones +template +class Convert1QC; + +/// Declare base converter into exp cones +template +class Convert1ExpC; /** A class to store a ref to ModelConverter. * @@ -62,34 +74,58 @@ class ConicConverter : public MCKeeper { MCKeeper(mc) { } /// Run conversions - void Run() { - if (MC().IfPassSOCPCones()) { - Walk(); - Walk(); - Walk(); - - Walk(); - Walk(); - Walk(); - } else - if (MC().IfPassQuadCon() && - (MC().GetNumberOfAddable((PowConstraint*)0)>0 || - MC().GetNumberOfAddable((AbsConstraint*)0)>0)) { - // Still collect QCones expressed by 2-norms. - // They are to be converted to quadratics. - Walk(); - Walk(); - Walk(); - } + void Run() { + // We _could_ walk everything just once + // and see which cones are there. + // Or, even walk expression trees from exponents, etc. + RunQCones(); + RunExpCones(); } protected: - /// Walk a single constraint type - template + void RunQCones() { + if (MC().IfPassSOCPCones()) { // convert everything to QuadraticCones. + Walk(); + Walk(); + Walk(); + + if (MC().GetNumberOfAddable((PowConstraint*)0)>0 || + MC().GetNumberOfAddable((AbsConstraint*)0)>0) { + Walk(); + Walk(); + Walk(); + } + } else + if (MC().IfPassQuadCon() && + (MC().GetNumberOfAddable((PowConstraint*)0)>0 || + MC().GetNumberOfAddable((AbsConstraint*)0)>0)) { + // Still collect QCones expressed by 2-norms. + // They are to be converted to quadratics. + Walk(); + Walk(); + Walk(); + } + } + + void RunExpCones() { + if (MC().IfPassExpCones() && // convert everything to ExpCones. + MC().GetNumberOfAddable((ExpConstraint*)0)>0) { + Walk(); + Walk(); // also ExpA ?? + Walk(); + + Walk(); + Walk(); + Walk(); + } + } + + /// Walk a single constraint type + template class CvtBase> void Walk() { auto& ck = MC().GetConstraintKeeper((Con*)nullptr); - ck.ForEachActive(Convert1 { MC() }); + ck.ForEachActive(Convert1 { MC() }); } /// Retrieve the MC @@ -98,7 +134,7 @@ class ConicConverter : public MCKeeper { /** Converter for a single linear / quadratic - * inequality constraint. + * inequality constraint into a Quadratic cone. * Generic, receives only body and sense. */ template @@ -561,17 +597,203 @@ class Convert1QC : public MCKeeper { }; +/** Converter for a single linear / quadratic + * inequality constraint into an exponential cone. + * Generic, receives only body and sense. + */ +template +class Convert1ExpC : public MCKeeper { +public: + /// Constructor + Convert1ExpC(MCType& mc) : MCKeeper(mc) { } + + +protected: + /// Characterize target subexpression. + template + struct SubexprTraits { + /// Construct + SubexprTraits( + std::array c={}, + std::array v={}, + std::vector v2d={}) + : coefs_(c), vars_(v), vars2del_(v2d) { } + /// Valid? + operator bool() const { return valid_; } + /// Variables to use. + /// Coefs + std::array coefs_; + /// -1 means 'use constant 1'. + std::array vars_; + /// Vars to delete + std::vector vars2del_; + /// Valid? + bool valid_ {false}; + }; + + /// LHS of ax >= by * exp( cz / by ) + using LhsTraits = SubexprTraits<1>; + /// RHS + using RhsTraits = SubexprTraits<2>; + + /// DoRun. Body: quadratic. + /// + /// Accept non-(+-1) coefficients. + /// + /// @param body: quadratic constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const QuadAndLinTerms& body, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + if (1 == body.GetQPTerms().size()) { + const auto& lt = body.GetLinTerms(); // But we have c0*x + c1*y*zz <=> 0 + const auto& qt = body.GetQPTerms(); + LhsTraits lhst; + if (0.0 == std::fabs(rhs) && // rhs==0: + 1 == lt.size()) // ax >= by exp(cz / by) ? + lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0)); + else if (0 == lt.size()) // const >= by exp(cz / by) ? + lhst = ClassifyLhs(-rhs*sens, -1); + + if (lhst) { + if (auto rhst = ClassifyRhsQuadr( + -qt.coef(0)*sens, qt.var1(0), qt.var2(0))) + return AddExpCone(lhst, rhst); + else if (auto rhst = ClassifyRhsQuadr( + -qt.coef(0)*sens, qt.var2(0), qt.var1(0))) + return AddExpCone(lhst, rhst); + } + } + return false; + } + + /// DoRun. + /// Body: linear, i.e., var1 (or const) >= var2 (or const). + /// In particular, the following cases. + /// /// + /// Considering const>=0. + /// Accept non-(+-1) coefficients. + /// + /// @param body: linear constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const LinTerms& lt, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + if (0.0 == std::fabs(rhs) && // rhs==0: y==const + 2 == lt.size()) { // ax >= b exp(cz / b) + if (auto rhst = ClassifyRhsLin(-lt.coef(0)*sens, lt.var(0))) { + // c2*x >= -c1*exp() + if (auto lhst = ClassifyLhs(lt.coef(1)*sens, lt.var(1))) { + return AddExpCone(lhst, rhst); + } + } else if (auto rhst = ClassifyRhsLin(-lt.coef(1)*sens, lt.var(1))) { + // c1*x >= -c2*exp() + if (auto lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0))) { + return AddExpCone(lhst, rhst); + } + } + } else if (1 == lt.size()) { // const >= b exp(cz / b) + if (auto lhst = ClassifyLhs(-rhs*sens, -1)) { + if (auto rhst = ClassifyRhsLin(-lt.coef(0)*sens, lt.var(0))) { + return AddExpCone(lhst, rhst); + } + } + } + return false; + } + + +protected: + LhsTraits ClassifyLhs(double a, int x) { + LhsTraits result {{a}, {x}}; + if ((x<0 && a>=0.0) || // Lhs is const, >=0 + (x>=0 && a>=0.0 && MC().lb(x)>=0.0) || // c*v >= 0 + (x>=0 && a<=0.0 && MC().ub(x)>=0.0)) + result.valid_ = true; + return result; + } + + /// The RHS is just b * (v=exp(z))? + RhsTraits ClassifyRhsLin(double b, int v) { + assert(v>=0); + RhsTraits result {{b, b}, {-1} /*constant*/}; + if (b>=0.0) { + if (auto pConExp = MC().template + GetInitExpressionOfType(v)) { + result.vars_[1] = pConExp->GetArguments()[0]; // the z + result.vars2del_ = {v}; // delete v + result.valid_ = true; + } + } + return result; + } + + /// The RHS is by * (v1 = exp( v2 = cz / by )) + /// or by * (v1 = exp( v2 = c / by ))? + /// In our case, c/by is stored as [constvar=c/b]/y + RhsTraits ClassifyRhsQuadr(double b, int y, int v1) { + assert(y>=0 && v1>=0); + RhsTraits result {{b, b}}; + if ((b>=0.0 && MC().lb(y)>=0.0) || // by >= 0 + (b<=0.0 && MC().ub(y)<=0.0)) { + if (auto pConExp = MC().template + GetInitExpressionOfType(v1)) { + int v2 = pConExp->GetArguments()[0]; // the z + if (auto pConDiv = MC().template + GetInitExpressionOfType(v2)) { + int z = pConDiv->GetArguments()[0]; // the z0 + int y0 = pConDiv->GetArguments()[1]; + if (y0==y) { // v2 = z/y + result.vars_ = {y, z}; + result.vars2del_ = {v1, v2}; // delete v1, v2 + result.valid_ = true; + } else if (auto pConLin = MC().template + GetInitExpressionOfType(y0)) { + const auto& ae = pConLin->GetAffineExpr(); + if (0.0 == std::fabs(ae.constant_term() && + 1==ae.GetBody().size())) { + const auto& body = ae.GetBody(); // can we ever get here? + if (y == body.var(0)) { // v2 = z / (c1*y) + result.coefs_ = {b, b/body.coef(0)}; + result.vars_ = {y, z}; + result.vars2del_ = {v1, v2, y0}; // delete + result.valid_ = true; + } + } + } + } + } + } + return result; + } + + bool AddExpCone(LhsTraits l, RhsTraits r) { + return true; + } + + /// Retrieve the MC + using MCKeeper::MC; +}; + + /** Converter for a single algebraic inequality constraint. * Receives range constraints and proceeds * if they are inequalities. */ -template +template class CvtBase> class Convert1 > : - public Convert1QC { + AlgebraicConstraint, CvtBase > : + public CvtBase { public: /// Constructor - Convert1(MCType& mc) : Convert1QC(mc) { } + Convert1(MCType& mc) : CvtBase(mc) { } /// Constraint type using ConType = AlgebraicConstraint; @@ -587,7 +809,7 @@ class Convert1(i) }; return - Convert1QC::DoRun(con.GetBody(), + CvtBase::DoRun(con.GetBody(), lbf ? 1 : -1, lbf ? con.lb() : con.ub()); } @@ -602,13 +824,14 @@ class Convert1 +template class CvtBase> class Convert1 > > : - public Convert1QC { + AlgebraicConstraint< Body, AlgConRhs >, CvtBase > : + public CvtBase { public: /// Constructor - Convert1(MCType& mc) : Convert1QC(mc) { } + Convert1(MCType& mc) : CvtBase(mc) { } /// The constraint type using ConType = @@ -621,7 +844,7 @@ class Convert1(i) }; return - Convert1QC::DoRun(con.GetBody(), + CvtBase::DoRun(con.GetBody(), sens, con.rhs()); }