How to formulate a general second-order cone constraint
AnsweredHi, I am inquiring if I can formulate a general second-order cone program (SOCP) as shown in Wikipedia:
min_x q.T * x
s.t. ||A_i * x + b_i || <= c_i.T * x + d_i
Fx = g
I am especially wondering how to formulate the cone constraint ||A_i * x + b_i || ≤ c_i.T * x + d_i. It seems that there is no API directly supporting it. Any help will be appreciated. Thank you!
For starters, here is my code generating a general random socp:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-import numpy as np
import gurobipy as gp
from gurobipy import GRB
import torchdef generate_random_socp(dim,nIneq,n_As,nEq,seed=1):
'''
Generate random general socp:
min_x q.T * x
s.t. ||A_i * x + b_i || <= c_i.T * x + d_i
Fx = gParameters
----------
dim : Dim of variable
nIneq : #Cone constraints
n_As : List, n_As[i]: # rows of A_i
nEq : #Equality constraints
seed : Random seed. The default is 1.Returns
-------
q : Torch tensor with size (dim).
A : List, A[i]: Torch tensor with size (n_As[i],dim).
b : List, b[i]: Torch tensor with size (n_As[i]).
c : List, c[i]: Torch tensor with size (dim).
d : List, d[i]: Torch tensor with size ().
F : Torch tensor with size (nEq,dim).
g : Torch tensor with size (nEq).'''
torch.manual_seed(seed)
q = torch.randn(dim,dtype=torch.float64)
A = []
b = []
c = []
d = []
x0 = torch.randn(dim,dtype=torch.float64)
for i in range(nIneq):
A_i = torch.randn(n_As[i], dim, dtype=torch.float64)
b_i = torch.randn(n_As[i], dtype=torch.float64)
A.append(A_i)
b.append(b_i)
c_i = torch.randn(dim, dtype=torch.float64)
c.append(c_i)
d_i = torch.norm(A[i] @ x0 + b[i], 2) - c[i].T @ x0
d_i = d_i.detach().clone()
d.append(d_i)
if nEq>0:
F = torch.randn(nEq, dim, dtype=torch.float64)
g = F @ x0
else:
F,g = None,None
return q,A,b,c,d,F,g
dim = 5
nIneq = 2
n_As = [dim] * nIneq
nEq = 2
seed = 0
q,A,b,c,d,F,g = generate_random_socp(dim,nIneq,n_As,nEq,seed=seed)
q = q.detach().numpy()
if F is not None:
F,g = F.detach().numpy(),g.detach().numpy()
A = [ Ai.detach().numpy() for Ai in A]
b = [ bi.detach().numpy() for bi in b]
c = [ ci.detach().numpy() for ci in c]
d = [ di.detach().numpy() for di in d]
-
-
Gurobi Staff
Hi Fengyu,
If you want to model
||A_i * x + b_i || <= c_i.T * x + d_ithen it should be sufficient to introduce an auxiliary non-negative variable v := c_i.T * x + d_i then square both sides. You can optionally introduce a free variable (y) to simplify the left hand side, egm = gp.Model() ... v = m.addVar() y = m.addVar(lb=-float("inf")) m.addConstr(v == c_i@x + d_i) m.addConstr(y == A_i@x + b_i) m.addConstr(y*y <= v*v) # soc- Riley
0 -
-
Hi Riley,
Thanks for your help. I tried this way but got “Model is infeasible or unbounded”:
Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (mac64[arm] - Darwin 22.3.0 22D49) CPU model: Apple M2 Thread count: 8 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 14 rows, 9 columns and 82 nonzeros Model fingerprint: 0xb0a81c5d Model has 2 quadratic constraints Coefficient statistics: Matrix range [4e-03, 2e+00] QMatrix range [1e+00, 1e+00] Objective range [3e-01, 2e+00] Bounds range [0e+00, 0e+00] RHS range [2e-02, 7e+00] Presolve time: 0.00s Barrier solved model in 0 iterations and 0.00 seconds (0.00 work units) Model is infeasible or unboundedHowever, the problem is feasible and can be solved by other models like CVXPY. Below is my code. Could you help me find the issue?
import numpy as np import gurobipy as gp from gurobipy import GRB q = np.array([ 1.54099611, -0.29342891, -2.17878938, 0.56843128, -1.08452234]) A = [np.array([[-0.55882344, 0.13559269, 0.2309147 , -0.61956913, 0.8847803 ], [-1.06620742, 0.15351884, -0.63327461, -0.94543633, -1.21768778], [-1.26546531, -1.11953518, 1.16648572, 0.92623246, 0.79701002], [ 0.21166418, 0.13557114, 1.87277638, -0.30991497, 0.24018023], [-1.34876103, 0.24439909, -0.0137207 , -0.19960823, -1.65642469]]), np.array([[ 0.11494421, 0.14466901, 1.07417172, -1.12879864, -0.09376858], [ 0.99834358, -0.08262167, 0.18253195, -0.29286147, -0.39490341], [ 0.44642257, 0.54011354, 0.28233708, -0.514967 , -1.46687432], [ 1.10965351, 0.27565496, -1.31818712, 0.88871446, 1.1200386 ], [ 0.60258537, -0.92537068, -1.3154607 , 1.96756662, -0.04131868]])] b = [np.array([-0.02403159, 0.64799465, 1.73401236, -0.13356864, -0.83510869]), np.array([-0.56646389, -1.15361677, -2.50230094, 0.87555775, -2.6725649 ])] c = [np.array([ 0.98554328, -0.00522718, -0.09443805, -1.12359376, -0.65126342]), np.array([-0.0313326 , 0.49878616, -0.52332327, -0.2514795 , -1.05553181])] d = [np.array(5.29530156), np.array(7.12259806)] F = np.array([[-0.55926127, -0.1197086 , -0.16354572, -0.25046325, -1.08822452], [ 0.00440725, 0.23278888, -0.82014159, -0.67901932, 1.42508336]]) g = np.array([ 1.21591643, -0.6859779 ]) dim = q.shape[0] nIneq = len(A) m = gp.Model() x = m.addMVar(shape=dim) for i in range(nIneq): A_i,b_i,c_i,d_i = A[i],b[i],c[i],d[i] y_i = m.addVar() z_i = m.addVar(lb=0) m.addConstr(z_i==c_i @ x + d_i) m.addConstr(y_i==A_i @ x + b_i) m.addConstr(y_i*y_i <= z_i*z_i) m.addConstr(F @ x == g) m.setObjective(q @ x, GRB.MINIMIZE) m.optimize()0 -
-
-
Hi Fengyu,
It will be good for you to read this article: How do I resolve the error "Model is infeasible or unbounded"?
Then follow it with How do I determine why my model is infeasible?
When you compute an IIS you will find the infeasibility is unrelated to the cones.
Note that the default lower bound for variables in Gurobi is 0 - this is a common cause of unexpected infeasibility.
- Riley
0 -
-
Hi Riley,
Thanks a lot for your help! Now I have fixed the infeasibility issue. I wonder if it is possible to get the dual variables corresponding to the cone constraint. Now I get an error at
nu = np.array([cone.Pi for cone in cones]) File src/gurobipy/mconstr.pxi:52 in gurobipy.MConstr.__getattr__ File src/gurobipy/mqconstr.pxi:169 in gurobipy.MQConstr.getAttr AttributeError: 'gurobipy.MQConstr' object has no attribute 'Pi'Besides, am I right that additional computation is needed for Gurobi to get the dual variable for QCP ? My updated code is attached. Thanks!
dim = q.shape[0] nIneq = len(A) m = gp.Model() x = m.addMVar(shape=dim,lb=-GRB.INFINITY) m.setObjective(q @ x, GRB.MINIMIZE) # m.Params.DualReductions = 0 m.Params.QCPDual = 1 cones = [] for i in range(nIneq): A_i,b_i,c_i,d_i = A[i],b[i],c[i],d[i] y_i = m.addMVar(shape=A_i.shape[0],lb=-GRB.INFINITY) z_i = m.addVar() m.addConstr(z_i == c_i @ x + d_i) m.addConstr(y_i == A_i @ x + b_i) cone_i = m.addConstr(gp.quicksum(y_i[j] * y_i[j] for j in range(y_i.shape[0])) <= z_i * z_i) cones.append(cone_i) if F is not None: eq_constr = m.addConstr(F @ x == g) m.optimize() x = x.X nu = np.array([cone.Pi for cone in cones]) mu = None if F is not None: mu = np.array([c.Pi for c in eq_constr])0 -
-
-
Hi Fengyu,
Yes, dual values are available for quadratic constraints. You can find “QCPi” and other quadratic constraint attributes here: https://docs.gurobi.com/projects/optimizer/en/current/reference/attributes/constraintquadratic.html#qcpi
- Riley
0 -
-
Hi Riley,
Thanks! I ran
nu = np.array([cone.QCPi for cone in cones]), then I got the dual variable for the quadratic constraints. But it is very different from the dual solutions corresponding to the cone conditions from other solvers, e.g. Cvxpy, which exactly formulates the problem as a SOCP. Is the reason Gurobi has to solve SOCP in a QCP format?
And I notice that I got this output when setting QCPDual=1,Solving KKT system to obtain QCP duals...Does it mean Gurobi did additional computation to retrieve the dual solutions?
Fengyu
1 -
-
-
Hi Fengyu,
Is the reason Gurobi has to solve SOCP in a QCP format?
I don't think so. The logs indicate the constraints are recognized as SOCP so then our SOCP solver is used. What is different about the dual values? The same primal solution could be dual degenerate giving rise to different dual values, or there could be multiple optimal primal solutions.
Does it mean Gurobi did additional computation to retrieve the dual solutions?
Correct. From QCPDual docs: “Computing them [QC dual values] can add significant time to the optimization, so you should only set this parameter to 1 if you need them.”
- Riley
0 -
Please sign in to leave a comment.
Comments
7 comments