Add test for rotated SOC constraint and support Gurobi via reformulation

docs/source/constraint.md

@@ -115,7 +115,7 @@ vars = [model.add_variable() for i in range(N)]
 con = model.add_second_order_cone_constraint(vars)
 ```
 
-Some optimizers (COPT, Mosek) support another form of second-order cone constraint called as rotated second-order cone constraint, which is defined as:
+There is another form of second-order cone constraint called as rotated second-order cone constraint, which is defined as:
 
 $$
 variables=(t_{1},t_{2},x) \in \mathbb{R}^{N} : 2 t_1 t_2 \ge \lVert x \rVert_2^2

src/pyoptinterface/_src/constraint_bridge.py

@@ -1,12 +1,17 @@
 from .core_ext import ScalarQuadraticFunction, ConstraintSense
 
 
-def bridge_soc_quadratic_constraint(model, cone_variables, name=""):
+def bridge_soc_quadratic_constraint(model, cone_variables, name="", rotated=False):
     """
     Convert a second order cone constraint to a quadratic constraint.
     x[0] >= sqrt(x[1]^2 + ... + x[n]^2)
     to
     x[0]^2 - x[1]^2 - ... - x[n]^2 >= 0
+    or
+    convert a rotated second order cone constraint to a quadratic constraint.
+    2 * x[0] * x[1] >= x[2]^2 + ... + x[n]^2
+    to
+    2 * x[0] * x[1] - x[2]^2 - ... - x[n]^2 >= 0
     """
     N = len(cone_variables)
     if N < 2:
@@ -17,13 +22,28 @@ def bridge_soc_quadratic_constraint(model, cone_variables, name=""):
     expr = ScalarQuadraticFunction()
     expr.reserve_quadratic(N)
 
-    x0 = cone_variables[0]
-    expr.add_quadratic_term(x0, x0, 1.0)
+    if not rotated:
+        x0 = cone_variables[0]
+        expr.add_quadratic_term(x0, x0, 1.0)
+
+        for i in range(1, N):
+            xi = cone_variables[i]
+            expr.add_quadratic_term(xi, xi, -1.0)
 
-    for i in range(1, N):
-        xi = cone_variables[i]
-        expr.add_quadratic_term(xi, xi, -1.0)
+        con = model.add_quadratic_constraint(
+            expr, ConstraintSense.GreaterEqual, 0.0, name
+        )
+    else:
+        x0 = cone_variables[0]
+        x1 = cone_variables[1]
+        expr.add_quadratic_term(x0, x1, 2.0)
 
-    con = model.add_quadratic_constraint(expr, ConstraintSense.GreaterEqual, 0.0, name)
+        for i in range(2, N):
+            xi = cone_variables[i]
+            expr.add_quadratic_term(xi, xi, -1.0)
+
+        con = model.add_quadratic_constraint(
+            expr, ConstraintSense.GreaterEqual, 0.0, name
+        )
 
     return con

tests/test_soc.py

@@ -45,3 +45,43 @@ def test_soc(model_interface):
     assert z_val == approx(4.0)
     obj_val = model.get_value(obj)
     assert obj_val == approx(9.5)
+
+
+def test_rotated_soc(model_interface):
+    model = model_interface
+
+    if not hasattr(model, "add_second_order_cone_constraint"):
+        pytest.skip("Model interface does not support second order cone constraints")
+
+    x = model.add_variable(lb=0.0, name="x")
+    y = model.add_variable(lb=0.0, name="y")
+    z = model.add_variable(lb=4.0, ub=4.0, name="z")
+
+    obj = x + 2 * y
+    model.set_objective(obj, poi.ObjectiveSense.Minimize)
+
+    con1 = model.add_second_order_cone_constraint([x, y, z], rotated=True)
+    model.optimize()
+    status = model.get_model_attribute(poi.ModelAttribute.TerminationStatus)
+    assert status == poi.TerminationStatusCode.OPTIMAL
+    x_val = model.get_value(x)
+    y_val = model.get_value(y)
+    assert x_val == approx(4.0, rel=1e-3)
+    assert y_val == approx(2.0, rel=1e-3)
+    obj_val = model.get_value(obj)
+    assert obj_val == approx(8.0, rel=1e-3)
+
+    model.delete_constraint(con1)
+    xx = model.add_variable(lb=0.0, name="xx")
+    model.add_linear_constraint(xx - 4 * x, poi.ConstraintSense.Equal, 0.0)
+
+    model.add_second_order_cone_constraint([xx, y, z], rotated=True)
+    model.optimize()
+    status = model.get_model_attribute(poi.ModelAttribute.TerminationStatus)
+    assert status == poi.TerminationStatusCode.OPTIMAL
+    x_val = model.get_value(x)
+    y_val = model.get_value(y)
+    assert x_val == approx(2.0, rel=1e-3)
+    assert y_val == approx(1.0, rel=1e-3)
+    obj_val = model.get_value(obj)
+    assert obj_val == approx(4.0, rel=1e-3)