Not getting the dual solution of SOCP model
Answered
Hi,
I'm trying to get the dual solutions of some linear constraints in my SOCP model. The model can be solved to optimal successfully by Gurobi. After setting "QCPDual" to 1, I get the following warning:
Q matrix is not positive semi-definite (PSD). Set NonConvex parameter to -1 or 2 to solve model.
However, SOC is inherently convex. When it is expressed in a general quadratic form, its Q matrix must be positive semi-definite, right? If I set "NonConvex" as 2, the model will be treated as MIP and I still cannot get dual solutions. My model is as follows:
I am wondering if my numerical values are rather complex, so I tried to use auxiliary variables to simplify the above two SOC constraints, as follows:
The new model is completely equivalent to the old model. However, the new model is identified as infeasible by Gurobi. How can I successfully obtain the dual variable values I need? Many thanks!
-
Hi Wenjia,
What version of Gurobi are you using? Was it working previously?
Cheers,
David0 -
Hi David,
My Gurobi was working fine before. Its version is:
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (mac64[arm] - Darwin 23.5.0 23F79)
Additionally, I need to correct a typo: in the first line (on the left hand side) of the original second SOC, the expression should include \( \delta_{. The new model with auxiliary variables is equivalent to the old model.
Best regards,
Wenjia
0 -
Hi David,
I apologize for forgetting to mention this: the initial SOCP problem was also considered non-convex by Gurobi, and it automatically solved it as MIP:
Optimize a model with 22 rows, 24121 columns and 51 nonzeros
Model fingerprint: 0xce8be06f
Model has 1925 quadratic constraints
Coefficient statistics:
Matrix range [1e+00, 3e+02]
QMatrix range [1e-09, 4e+04]
QLMatrix range [1e-05, 4e+04]
Objective range [1e+00, 1e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+03]
QRHS range [3e-11, 7e+00]
Warning: Quadratic constraints contain large coefficient range
Consider reformulating model or setting NumericFocus parameter
to avoid numerical issues.
Presolve removed 21 rows and 22194 columns
Continuous model is non-convex -- solving as a MIPI don't understand why an SOCP would be considered non-convex. Is this related to the numerical issues it mentioned? Thank you in advance for your help.
Best regards,
Wenjia
0 -
This may also come down to the implementation in the code. For example, this is convex:
m = gp.Model()
x = m.addVar()
y = m.addVar()
m.setObjective(x**2 + y**2, gp.GRB.MINIMIZE)
m.optimize()We tend to say this is an equivalent model:
m = gp.Model()
x = m.addVar()
y = m.addVar()
w = m.addVar()
z = m.addVar()
m.addConstr(w == x**2)
m.addConstr(z == y**2)
m.setObjective(w + z, gp.GRB.MINIMIZE)
m.optimize()but ultimately the feasible region is not convex in the (x,y,w,z) space - the projection onto (x,y) is convex - but currently Gurobi does not recognize this (this will likely change in a future version).
Sometimes using epigraphs can work around the issue, eg the following is convex:
m = gp.Model()
x = m.addVar()
y = m.addVar()
w = m.addVar()
z = m.addVar()
m.addConstr(w >= x**2)
m.addConstr(z >= y**2)
m.setObjective(w + z, gp.GRB.MINIMIZE)
m.optimize()Are you using a Gurobi API or 3rd party framework?
0 -
Hi Riley,
Thank you very much for your guidance! Using epigraphs is indeed useful. I am using Gurobi through Python with the gurobipy package.
However, I still hope to know why Gurobi considers the model non-convex when I code according to the initial detailed SOC, even after I eliminate the numerical issues?
Best,
Wenjia
0 -
My Gurobi was working fine before.
Do you mean it was working with a previous version of Gurobi? Or with a previous/equivalent formulation of your problem?
Cheers,
David0 -
Hi David,
Sorry I misunderstood your previous question. I meant that I have only used version 11.0, and in other models, the use of auxiliary variables has not caused infeasibility.
For either the previous or the equivalent models I presented this time, Gurobi does not work as I expect because it keeps indicating non-convexity or infeasibility.
Best,
Wenjia0 -
Hi Wenjia,
I don't understand why an SOCP would be considered non-convex
I think we'd need to see the code to be able to comment further. Feel free to drop a link to file in DropBox, Google Drive, github etc.
- Riley
0 -
Hi Riley,
Please refer to this Google Drive link: https://drive.google.com/file/d/17RFhDFcoQp57Gf3S5F1R5EEGH6XEWtFa/view?usp=share_link
It includes the code ("test.py") and data files. Thank you for your help in advance!
Best,
Wenjia
0 -
Hi Wenjia,
From a quick glance it looks like the issue is this:
\[\left\| \begin{bmatrix} a + b \\ c \end{bmatrix} \right\|^2 = (a+b)^2 + c^2\]
\[\left\| \begin{bmatrix} a + b \\ c \end{bmatrix} \right\|^2 \neq a^2+b^2 + c^2\]
In particular I am referring to the sums over H inside your matrices.
Maybe the simplest fix is to define
\[dd_{ia} = \sum_hdd_{iah}\]
and use \(dd_{ia}\) inside the matrix, and do something similiar for \(aa_{ih}\).
Also I notice you set the lower bound of variables to 0 for a couple of variable families. This raises flag for me - the default is already 0 - are you intending the others to have a lower bound of -infinity?
- Riley
0 -
Hi Riley,
Actually, taking the second SOC as an example, the sum_h part is like this:
Therefore, I should calculate a^2 + b^2 (as I've written in the code), instead of (a+b)^2. I apologize for my previous notation not being clear enough. “{...}_{h \in H}” means that there are H such terms.
Thank you very much for mentioning the lower bound issue. I overlooked this point. Only in the equivalent model, the variables bb and ee may be negative. For all the other variables, they should all be >= 0.
Best,
Wenjia
0 -
Hi Wenjia,
{...}_{h \in H}” means that there are H such terms
Ah I see, this makes the difference, I assumed you meant to sum them.
I'll try and make some time to investigate on the weekend.
- Riley
0 -
Hi Riley,
Due to my unfamiliarity with the default settings of the `model.addVar()` function, the equivalent model was previously deemed infeasible. After I adjusted the lower bounds of `bb` and `ee` to -infinity, the equivalent model is recognized as an SOCP and successfully solved. I am very grateful for your help!
Currently, the remaining issue is only related to the original SOC model. Is it because the original SOC form is too complex for Gurobi to recognize?
Thank you again for your understanding and support. Your guidance is invaluable to me.
Best regards,
Wenjia0 -
Hi Wenjia,
I think that the issue with the original implementation is partly related to the fact that the constraints are being squared on both sides. We can still use this form but we need to restrict the domain.
\[\begin{array}{lr}\{x:||Ax+b|| \leq c^Tx + d\} & \quad\text{convex set}\\\{x:(Ax+b)^T(Ax+b) \leq (c^Tx + d)^2\} & \quad\text{not convex set}\\\{x:(Ax+b)^T(Ax+b) \leq (c^Tx + d)^2, c^Tx + d \geq 0\} &\quad\text{convex set}\end{array}\]
A simpler case shows why:
Let's take the first SOC as an example:
for i in range(I):
model.addQConstr(4 * omega[i, t, j] * np.max(r) * w0 +
quicksum(4 * (omega[i, t, j] * np.max(r) - omega[i, t, j] * r[h]) * w[i, h] * (x[h] ** 2) for h in range(H)) +
(w0 + quicksum(w[i, h] * x[h] for h in range(H)) - omega[i, t, j] * np.max(r) + q[i]) ** 2
<= (w0 + quicksum(w[i, h] * x[h] for h in range(H)) + omega[i, t, j] * np.max(r) - q[i]) ** 2, f"SOC1")In theory it should be enough to add:
model.addConstr(w0 + quicksum(w[i, h] * x[h] for h in range(H)) + omega[i, t, j] * np.max(r) - q[i] >= 0)
to get a convex region, however I have found that I actually need to enforce this through a bound on an auxiliary variable, and using that variable in the constraint, i.e.
for i in range(I):
tmp = model.addVar() # default lb=0, continuous
model.addConstr(tmp == w0 + quicksum(w[i, h] * x[h] for h in range(H)) + omega[i, t, j] * np.max(r) - q[i])
model.addQConstr(4 * omega[i, t, j] * np.max(r) * w0 +
quicksum(4 * (omega[i, t, j] * np.max(r) - omega[i, t, j] * r[h]) * w[i, h] * (x[h] ** 2) for h in range(H)) +
(w0 + quicksum(w[i, h] * x[h] for h in range(H)) - omega[i, t, j] * np.max(r) + q[i]) ** 2
<= tmp ** 2, f"SOC1")If I do a similar thing for the second SOC then Gurobi recognizes the problem as convex.
As for why it is not enough to simply add the RHS >= 0 constraint, I suspect it is because convexity is determined from the Q matrix, and this constraint is separate to that, but using the auxiliary variable allows the bound to affect the PSD calculation.
- Riley
0 -
Hi Riley,
Got it. Thank you so much! Your guidance is exactly what I need to solve the problem.
Best,
Wenjia
0
Please sign in to leave a comment.
Comments
15 comments