• Gurobi Staff

Hi Wenjia,

What version of Gurobi are you using? Was it working previously?

Cheers,
David

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 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 nonzerosModel fingerprint: 0xce8be06fModel has 1925 quadratic constraintsCoefficient 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 columnsContinuous model is non-convex -- solving as a MIP I 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 • Gurobi Staff 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? 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 • Gurobi Staff 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, David 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, Wenjia • Gurobi Staff 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 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 • Gurobi Staff 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

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

• Gurobi Staff

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

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,
Wenjia

• Gurobi Staff

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

Hi Riley,

Got it. Thank you so much! Your guidance is exactly what I need to solve the problem.

Best,

Wenjia