Exp cones #192

CHANGES.mp.md

@@ -1,6 +1,13 @@
 Summary of recent updates to the AMPL MP Library
 ================================================
 
+
+## 20230515
+- *Recognize exponential conic constraints*.
+  Exponential cones are recognized and passed to the
+  solver, if supported.
+
+
 ## 20230424
 - *Pass variable names* if read from a `col` file with the 
   same name of the `nl` file being read.

CMakeLists.txt

@@ -214,7 +214,7 @@ endmacro()
 include_directories(include)
 
 set(CMAKE_POSITION_INDEPENDENT_CODE TRUE)
-set(MP_DATE 20230426)
+set(MP_DATE 20230515)
 
 set(MP_SYSINFO "${CMAKE_SYSTEM_NAME} ${CMAKE_SYSTEM_PROCESSOR}")
 

include/mp/flat/constr_general.h

@@ -190,16 +190,16 @@ using ComplementarityQuadratic = ComplementarityConstraint<QuadraticExpr>;
 
 /// Quadratic cone
 DEF_STATIC_CONSTR_WITH_PRM( QuadraticConeConstraint, VarArray, DblParamArray,
-														"Quadratic cone p1*x1 >= sqrt((p2*x2)^2 + ...)),"
+                            "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 p1*x1*p2*x2 >= sqrt((p3*x3)^2 + ...)),"
+                            "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 p1*x1 >= p2*x2*exp(p3*x3 / (p2*x2)),"
-                            " x1, x2 >= 0, with factors p1..p3");
+                            "Exponential cone ax >= by exp(cz / (by)),"
+                            " where ax, by >= 0, with factors a,b,c");
 /// Power cone
 DEF_STATIC_CONSTR_WITH_PRM( PowerConeConstraint, VarArray, DblParamArray,
                             "Power cone with factors ");

include/mp/flat/redef/conic/cones.h

@@ -722,7 +722,7 @@ class Convert1ExpC : public MCKeeper<MCType> {
   /// The RHS is just b * (v=exp(z))?
   RhsTraits ClassifyRhsLin(double b, int v) {
     assert(v>=0);
-    RhsTraits result {{b, b}, {-1} /*constant*/};
+    RhsTraits result {{b, b}, {-1} /*y = constant 1*/};
     if (b>=0.0) {
       if (auto pConExp = MC().template
           GetInitExpressionOfType<ExpConstraint>(v)) {
@@ -774,6 +774,20 @@ class Convert1ExpC : public MCKeeper<MCType> {
   }
 
   bool AddExpCone(LhsTraits l, RhsTraits r) {
+    assert(l.vars2del_.empty());
+    for (auto v2d: r.vars2del_)       // unuse result vars
+      MC().DecrementVarUsage(v2d);
+    std::array<int, 3> args
+    {l.vars_[0], r.vars_[0], r.vars_[1]};
+    for (int i=0; i<3; ++i) {
+      if (args[i]<0) {                // marked as -1?
+        args[i] = MC().MakeFixedVar(1.0);
+      }
+    }
+    MC().AddConstraint(
+          ExponentialConeConstraint(
+            args,
+            {l.coefs_[0], r.coefs_[0], r.coefs_[1]}));
     return true;
   }
 

solvers/mosek/CHANGES.mosek.md

@@ -1,9 +1,16 @@
 Summary of recent updates to MOSEK for AMPL
 ===========================================
 
+
+## 20230515
+- *Exponential cones*. MP driver recognizes exponential
+  cones from their algebraic representation and passes
+  them to Mosek.
+
+
 ## 20230505
 - *Updated Mosek library* to version 10.0.43. It includes a 
-  number of bug fixes, among which a numerical problem that 
+  number of bug fixes, including a numerical problem that 
   could occur with disjunctive constraints
 
 

solvers/mosek/mosekmodelapi.cc

@@ -261,6 +261,26 @@ void MosekModelAPI::AddConstraint(
 	MOSEK_CCALL( MSK_appendaccseq(lp(), domidx, nnz, numafe_prev, NULL) );
 }
 
+void MosekModelAPI::AddConstraint(
+    const ExponentialConeConstraint& ec) {
+  MSKint64t numafe_prev=0;
+  MOSEK_CCALL( MSK_getnumafe(lp(), &numafe_prev) );
+  auto nnz = ec.GetArguments().size();
+  // Is it too slow to add cones 1 by 1?
+  MOSEK_CCALL( MSK_appendafes(lp(), nnz) );
+  std::vector<MSKint64t> afeidx(nnz);
+  for (auto idx=nnz; idx--; )        // Fill new AFE indexes
+    afeidx[idx] = numafe_prev+idx;
+  MOSEK_CCALL(
+        MSK_putafefentrylist(lp(), nnz,
+                             afeidx.data(),
+                             ec.GetArguments().data(),  // assumes it's int32
+                             ec.GetParameters().data()) );
+  MSKint64t domidx=-1;
+  MOSEK_CCALL( MSK_appendprimalexpconedomain(lp(), &domidx) );
+  MOSEK_CCALL( MSK_appendaccseq(lp(), domidx, nnz, numafe_prev, NULL) );
+}
+
 /// Add indicator as disjunction
 /// (binvar!=binval) \/ (lincon).
 /// 1st afe: binvar+binval-1=0

solvers/mosek/mosekmodelapi.h

@@ -52,38 +52,40 @@ class MosekModelAPI :
 
 
   /// The linear range constraint, if fully supported with basis info etc.
-	ACCEPT_CONSTRAINT(LinConRange, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(LinConRange, Recommended, CG_Algebraic)
   void AddConstraint(const LinConRange& lc);
 
   /// LinCon(LE/EQ/GE) should have 'Recommended' for all backends
   /// and have an implementation,
   /// or a conversion rule is needed in a derived FlatConverter
-	ACCEPT_CONSTRAINT(LinConLE, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(LinConLE, Recommended, CG_Algebraic)
   void AddConstraint(const LinConLE& lc);
-	ACCEPT_CONSTRAINT(LinConEQ, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(LinConEQ, Recommended, CG_Algebraic)
   void AddConstraint(const LinConEQ& lc);
-	ACCEPT_CONSTRAINT(LinConGE, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(LinConGE, Recommended, CG_Algebraic)
   void AddConstraint(const LinConGE& lc);
 
   /// QuadConRange is optional.
-	ACCEPT_CONSTRAINT(QuadConRange, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(QuadConRange, Recommended, CG_Algebraic)
   void AddConstraint(const QuadConRange& qc);
 
   /// If using quadratics,
   /// QuadCon(LE/EQ/GE) should have 'Recommended'
   /// and have an implementation.
-	ACCEPT_CONSTRAINT(QuadConLE, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(QuadConLE, Recommended, CG_Algebraic)
   void AddConstraint(const QuadConLE& qc);
-	ACCEPT_CONSTRAINT(QuadConEQ, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(QuadConEQ, Recommended, CG_Algebraic)
   void AddConstraint(const QuadConEQ& qc);
-	ACCEPT_CONSTRAINT(QuadConGE, Recommended, CG_Algebraic)
+  ACCEPT_CONSTRAINT(QuadConGE, Recommended, CG_Algebraic)
   void AddConstraint(const QuadConGE& qc);
 
-	/// Cones
-	ACCEPT_CONSTRAINT(QuadraticConeConstraint, Recommended, CG_Conic)
-	void AddConstraint(const QuadraticConeConstraint& qc);
-	ACCEPT_CONSTRAINT(RotatedQuadraticConeConstraint, Recommended, CG_Conic)
-	void AddConstraint(const RotatedQuadraticConeConstraint& qc);
+  /// Cones
+  ACCEPT_CONSTRAINT(QuadraticConeConstraint, Recommended, CG_Conic)
+  void AddConstraint(const QuadraticConeConstraint& qc);
+  ACCEPT_CONSTRAINT(RotatedQuadraticConeConstraint, Recommended, CG_Conic)
+  void AddConstraint(const RotatedQuadraticConeConstraint& qc);
+  ACCEPT_CONSTRAINT(ExponentialConeConstraint, Recommended, CG_Conic)
+  void AddConstraint(const ExponentialConeConstraint& ec);
 
   /// Linear indicator constraints can be used as
   /// auxiliary constraints for logical conditions.
@@ -100,9 +102,9 @@ class MosekModelAPI :
   /// Set ``option pl_linearize 0;`` in AMPL if the solver
   /// supports PL natively.
   /// MOSEK 10 has no SOS
-	ACCEPT_CONSTRAINT(SOS1Constraint, NotAccepted, CG_General)
+  ACCEPT_CONSTRAINT(SOS1Constraint, NotAccepted, CG_General)
   void AddConstraint(const SOS1Constraint& cc);
-	ACCEPT_CONSTRAINT(SOS2Constraint, NotAccepted, CG_General)
+  ACCEPT_CONSTRAINT(SOS2Constraint, NotAccepted, CG_General)
   void AddConstraint(const SOS2Constraint& cc);
 
 

test/end2end/cases/categorized/fast/conic/expcones_01__plain.mod

@@ -0,0 +1,15 @@
+# expcones_01__plain.mod
+# From https://docs.mosek.com/latest/capi/tutorial-ceo-shared.html
+
+var x {1..3};
+
+s.t. Nonneg {i in 1..2}: x[i] >= 0;
+
+minimize Obj:
+   x[1] + x[2];
+
+s.t. ExpCone:
+   x[1] >= x[2] * exp( x[3] / x[2] );
+
+s.t. LinCon:
+   sum {i in 1..3} x[i] == 1;

test/end2end/cases/categorized/fast/conic/expcones_02.mod

@@ -0,0 +1,14 @@
+# expcones_02.mod
+
+var x {1..3};
+
+s.t. Nonneg {i in 1..2}: x[i] >= 0;
+
+minimize Obj:
+   x[1] + x[2];
+
+s.t. ExpCone:
+   2*x[1] >= 3*x[2] * exp( x[3] / (3*x[2]) );
+
+s.t. LinCon:
+   sum {i in 1..3} x[i] == 1;

test/end2end/cases/categorized/fast/conic/expcones_03.mod

@@ -0,0 +1,14 @@
+# expcones_03.mod
+
+var x {1..3};
+
+s.t. Nonneg {i in 1..2}: x[i] >= 0;
+
+minimize Obj:
+   x[1] + x[2];
+
+s.t. ExpCone:
+   2*x[1] >= 3*x[2] * exp( 4*x[3] / (3*x[2]) );
+
+s.t. LinCon:
+   sum {i in 1..3} x[i] == 1;

test/end2end/cases/categorized/fast/conic/expcones_04__negvar.mod

@@ -0,0 +1,14 @@
+# expcones_04__negvar.mod
+
+var x {1..3};
+
+s.t. Neg {i in 1..2}: x[i] <= 0;
+
+maximize Obj:
+   x[1] + x[2];
+
+s.t. ExpCone:
+   -2*x[1] >= -3*x[2] * exp( 4*x[3] / (-3*x[2]) );
+
+s.t. LinCon:
+   -x[1]-x[2]+x[3] == 1;

test/end2end/cases/categorized/fast/conic/expcones_05__const.mod

@@ -0,0 +1,20 @@
+# expcones_05__const.mod
+
+var x {1..3};
+
+s.t. Neg {i in 1..2}: x[i] <= 0;
+
+maximize Obj:
+   x[1] + x[2];
+
+s.t. ExpCone01:
+   -1.1 <= 3*x[2] * exp( 4*x[3] / (-3*x[2]) );
+
+s.t. ExpCone02:
+   -2*x[1] >= 0.3 * exp( 4*x[3] / (0.3) );
+
+s.t. ExpCone03:
+   1.1 >= -3*x[2] * exp( 0.4 / (-3*x[2]) );
+
+s.t. LinCon:
+   -x[1]-x[2]+x[3] == 1;

test/end2end/cases/categorized/fast/conic/modellist.json

@@ -50,5 +50,30 @@
     "name" : "socp_08_sqrt_constants",
     "objective" : 11.42915114,
     "tags" : ["socp"]
+  },
+  {
+    "name" : "expcones_01__plain",
+    "objective" : 0.7821882953,
+    "tags" : ["expcones"]
+  },
+  {
+    "name" : "expcones_02",
+    "objective" : 0.6242239555,
+    "tags" : ["expcones"]
+  },
+  {
+    "name" : "expcones_03",
+    "objective" : 0.8691893611,
+    "tags" : ["expcones"]
+  },
+  {
+    "name" : "expcones_04__negvar",
+    "objective" : -0.8691893611,
+    "tags" : ["expcones"]
+  },
+  {
+    "name" : "expcones_05__const",
+    "objective" : -0.8988331541,
+    "tags" : ["expcones"]
   }
 ]

test/end2end/scripts/python/Model.py

@@ -13,6 +13,7 @@ class ModelTags(enum.Enum):
     quadraticnonconvex = 3
     socp = 4
     socp_hard_to_recognize = 4.1  ## For solvers recognizing from quadratics
+    expcones = 4.4
     nonlinear = 5
     complementarity = 6
     arc = 7

test/end2end/scripts/python/Solver.py

@@ -808,6 +808,7 @@ def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
 
                  ModelTags.quadratic, ModelTags.quadratic_obj,
                  ModelTags.socp, ModelTags.socp_hard_to_recognize,
+                 ModelTags.expcones,
 
                  ModelTags.warmstart, ModelTags.mipstart,
                  ModelTags.return_mipgap, ModelTags.sens, ModelTags.sstatus}