Hyperboloid of two sheets
Answered
Dear Gurobi-Team,
if I model the following problem, Gurobi can solve it with GRB.Param.NonConvex=-1, although the Q-constraint-matrix is non-convex (due to the positive and negative eigenvalues of Q).
Since the variable c is restricted to positive values, the constraint is interpreted as a second-order cone constraint as described in the gurobi documentation as well as in the solver log:
https://www.gurobi.com/documentation/10.0/refman/constraints.html
import gurobipy as gp
from gurobipy import GRB
m = gp.Model("qcp")
x = m.addVar(name="x", lb=-1e3, ub=1e3)
y = m.addVar(name="y", lb=-1e3, ub=1e3)
z = m.addVar(name="z", lb=0, ub=1e3)
obj = x**2 + y**2 + (z+10)**2
m.setObjective(obj, GRB.MINIMIZE)
m.addConstr(x**2 + y**2 - z**2 <= -1, "qc0")
m.optimize()
It even works with the variable z restricted to negative values:
z = m.addVar(name="z", lb=-1e3, ub=0)
Strictly speaking, I would not call the geometric form of the exemplary quadratic constraint a second-order cone, but as one sheet of a Hyperboloid of two sheets, because the right side of the quadratic constraint is not zero, but negative.
If I now release the feasible domain of the variable z to positive and negative values, I receive the non-PSD GurobiError:
z = m.addVar(name="z", lb=-1e3, ub=1e3)
GurobiError: Q matrix is not positive semi-definite (PSD). Set NonConvex parameter to 2 to solve model.
That's reasonable and with GRB.Param.NonConvex=2 the problem can be solved:
m.setParam(GRB.Param.NonConvex, 2)
With GRB.Param.NonConvex=2, the problem can even be solved with positive values on the right side of the quadratic constraint. Then the geometric form is a (non-convex) Hyperboloid of one sheet.
m.addConstr(x**2 + y**2 - z**2 <= 1, "qc0")
--------------------------------------------------------------------------
Having all that said, I would like to come to my question:
Do I miss something or is it possible to reformulate the quadratic constraint, given as Hyperboloid of two sheets (negative right side of the constraint), as described in the following to solve it with GRB.Param.NonConvex=-1?
Then all advantages of mixed-integer convex optimization would remain.
Due to the symmetry of the geometric form, the algebraic sign of the values of z does not play a significant role in the constraint. The algebraic sign only gives the information in which of the two mirror-image sheets the solution is. So with an auxiliary variable zpos and an ABS-constraint (the ABS-constraint is simplified during presolve) it works for the given example, and I think the approach should work in general for Hyperboloids of two sheets (but not for Hyperboloids of one sheets).
import gurobipy as gp
from gurobipy import GRB
m = gp.Model("qcp")
x = m.addVar(name="x", lb=-1e3, ub=1e3)
y = m.addVar(name="y", lb=-1e3, ub=1e3)
z = m.addVar(name="z", lb=-1e3, ub=1e3)
zpos = m.addVar(name="zpos", lb=0, ub=1e3)
obj = x**2 + y**2 + (z+10)**2
m.setObjective(obj, GRB.MINIMIZE)
m.addConstr(x**2 + y**2 - zpos**2 <= -1, "qc0")
m.addGenConstrAbs(zpos, z, "absconstr")
m.optimize()
Thank you for the interest in this matter.
Thomas
-
Hi Thomas,
You are correct that you can reformulate the problem this way.
With reference to the 2 sheet hyperboloid, the epigraph of the upper surface is a convex set, as is the hypograph of the lower surface. When z is unbounded the optimization is over the union of 2 disjoint convex sets - which is of course nonconvex. If the domain of z is restricted to be nonnegative (or nonpositive) then the optimization is over a single convex set, which allows the problem to be solved without setting the NonConvex parameter to 2 in these cases.
Since the variable c is restricted to positive values, the constraint is interpreted as a second-order cone constraint as described in the gurobi documentation as well as in the solver log. (1)
I would not call the geometric form of the exemplary quadratic constraint a second-order cone (2)
I agree with (2) but I'm not sure I'm onboard with (1). My understanding is that the presolve routine produces a new model that includes second-order cone constraints (even if you try to turn presolve off) but this is not the same as saying that the quadratic constraint is a second-order cone constraint.
Do I miss something or is it possible to reformulate the quadratic constraint ...
Possibly just a matter of semantics but I think it would be more correct to say the constraint remains unchanged, but the restricted domain of z (assigned to zpos in your code) allows us to optimize over a convex set. The absolute value constraint requires a binary variable to be added to model this disjunction over a part of the objective function. You could also model this disjunction explicitly with indicator constraints:
# Define convex epigraph of upper surface of 2-sheet hyperboloid
m = gp.Model()
x = m.addVar(name="x", lb=-1e3, ub=1e3)
y = m.addVar(name="y", lb=-1e3, ub=1e3)
z = m.addVar(name="z", lb=0, ub=1e3)
m.addConstr(x**2 + y**2 - z**2 <= -1, "qc0")
# model disjunction over part of the objective
z_free = m.addVar()
w = m.addVar(vtype="B")
m.addConstr((w == 1) >> (z_free - z == 0))
m.addConstr((w == 0) >> (z_free + z == 0))
m.setObjective(x**2 + y**2 + (z_free+10)**2, GRB.MINIMIZE)Alternatively we could also model the disjunction over the entire objective function. I.e. we are now allowing two objective functions (that could be completely different in structure), and the optimal solution will use whichever gives the better result.
obj = m.addVar()
obj1 = m.addVar()
obj2 = m.addVar()
m.addConstr(obj1 >= x**2 + y**2 + (z+10)**2)
m.addConstr(obj2 >= x**2 + y**2 + (-z+10)**2)
w = m.addVar(vtype="B")
m.addConstr((w == 1) >> (obj - obj1 == 0))
m.addConstr((w == 0) >> (obj - obj2 == 0))
m.setObjective(obj, GRB.MINIMIZE)You could also extend this to a disjunction over the two convex regions of the hyperboloid. That is, model the upper region with x1, y1, z1, with objective obj1 defined in the previous example, then model the lower region with x2, y2, z2 with obj2. This approach is not as elegant as your compact formulation but it leads into the next point which is that we can model the union of convex regions without appealing to symmetry.
Say we want to solve the following problem over a non-convex region made from the union of two convex regions, defined by a circle and parabola:
\begin{eqnarray}
&& \min_{(x,y) \in S_1 \cup S_2} x^2\\
&& S_1 = \{(x,y): y^2 + x^2 - 6x \leq -8\}\\
&& S_2 = \{(x,y): x \leq -y^2 - 2y - 2)\}
\end{eqnarray}
This can be modeled like so:
m = gp.Model()
# circular region
x1 = m.addVar(name="x1", lb=0, ub=1e3)
y1 = m.addVar(name="y1", lb=-1e3, ub=1e3)
m.addConstr(y1*y1 + x1*x1 - 6*x1 <= -8)
# parabolic region
x2 = m.addVar(name="x2", lb=-1e3, ub=0)
y2 = m.addVar(name="y2", lb=-1e3, ub=1e3)
m.addConstr(x2 <= -y2*y2 - 2*y2 - 2)
# objective functions
obj = m.addVar()
obj1 = m.addVar()
obj2 = m.addVar()
m.addConstr(obj1 >= x1*x1)
m.addConstr(obj2 >= x2*x2)
# disjunction
w = m.addVar(vtype="B")
m.addConstr((w == 1) >> (obj - obj1 == 0))
m.addConstr((w == 0) >> (obj - obj2 == 0))
m.setObjective(obj, GRB.MINIMIZE)I hope this helps clarify different approaches to your problem.
- Riley
0 -
Hi Riley,
thank you very much for your detailed answer! Your explanations opened my eyes to further details and interesting possibilities!
I like precision a lot so I would like to thank you for clarifying that presolve produces a modified model with SOC and not the original constraint is or becomes a SOC or is interpreted as such.
And it is also understood that the constraint (consisting of a hyperboloid of two sheets) remains the same, but only the modeling of it is done in a way that optimization can be done over two separate convex domains.
Since we agree that disjoint convex regions can be modeled using convex sets and binary variables, I would like to suggest the feature that Gurobi identifies at least standard cases like the hyperboloid of two sheets and internally produces a new model during presolve that includes convex sets and binary variables in one of the ways discussed here (or otherwise). I think that extended capability of Gurobi would be beneficial for many users working with quadratic functions. And it doesn't look any harder to me than the things Gurobi already does with SOCs or bilinear transformations.
Best regards
Thomas
0 -
Hi Thomas,
You're welcome. I've applied the feature request tag to your post and I'll pass these thoughts onto the developers for consideration.
- Riley
0
Please sign in to leave a comment.
Comments
3 comments