• Gurobi Staff

Hi Carlos,

To model the absolute values, I have tried using gurobipy's abs_() function as well as manually modelling it using auxiliary binary variables, which did not help. (By the way, is abs_() in general "smarter" than just doing the canonical auxiliary variables to model absolute values?)

You did it the correct way. Using gurobipy's abs_() is usually smarter in the sense that Gurobi will internally compute (hopefully) tight big-M values for the binary formulation. If however, you can provide best possible big-M values for the binary formulation of the ABS function, then you can just write the formulation yourself.

Tried to avoid absolute values altogether by modeling a quadratic optimization function where each term is (2 * z_i - 1) * x_i, with z_i binary and x_i in [0, 1]. Also did not seem to help.

This is almost always a very bad idea when x is continuous, because it introduces more room for numerical inaccuracies without gaining any relaxation tightness.

Used the automated parameter tuning tool. Specifically, did it for the k = 2 instances then assumed the same parameters would help the larger instances (is this usually true?). They seem to have some effect but not large enough to solve larger instances.

Usually, parameters found by tuning are also useful for larger instances as long as the model structure stays the same, so what you see is correct behavior.

I had a look at the model and the solutions found for k=1,2, and 3. It looks like you fix one variable to a given value and then in the optimal solution all other x attain the same value in absolute terms, i.e., $$|x_i| = |x_j|$$ except for the one you fixed. Is this some property that may be true?

Altogether, I think that you should try to deduce some more properties of the model, which would help the solver, by making the formulation more precise. Our Tech-Talk on Weak vs Strong MIP formulations might be helpful here.

Best regards,
Jaromił

Hi Jaromił,

Thanks a lot for the help and the tech-talk recommendation! In case you or anyone else can help me further, I investigated the model a bit further and for $k \geq 2$ something similar to your suggestion seems to be true:

1. The value of $x_1$ can be predetermined.
2. It seems to be true that $\lvert x_i \rvert = \lvert x_j \rvert$ for $i = 2, \ldots, 4^k$

The straightforward way of adding this to the model is to have some constraints model.addConstrs((aux[i] == aux[i + 1] for i in range(1, ncoeffs - 1))) (see function mwe_equals below).

• Is there a smarter way to do this?

There is also another approach: Since the coefficients are all equal, I can try to guess the optimal value of the function and then run a binary feasibility program where I only try to find an attributed of $\pm 1$ signs that satisfies all the constraints (an attempt is in function mwe_signs).

• Is this a better idea than before, or are feasibility problems not recommended?

For both cases, I can solve $k = 2, 3$ much faster, but $k=4$ still seems difficult (as in 15min were not enough to get a solution so far).

from itertools import repeatfrom functools import reduceimport numpy as npimport gurobipy as gpfrom gurobipy import GRB# k > 1def mwe_equals(k):    ncoeffs = 4 ** k    first = 1 / ncoeffs    A = reduce(np.kron, repeat(A0, k))    B = reduce(np.kron, repeat(B0, k))    model = gp.Model()    x = model.addMVar((ncoeffs,), lb=-1, ub=1, vtype=GRB.CONTINUOUS, name="x")    model.addConstr(x[0] == first)    y = A @ x    model.addConstrs((y[i] >= 0 for i in range(ncoeffs)))    model.addConstr(y.sum() == 1)    tmp = A @ B @ x    model.addConstrs((tmp[i] >= 0 for i in range(ncoeffs)))    aux = model.addMVar(ncoeffs, vtype=GRB.CONTINUOUS, name="aux")    model.addConstrs((aux[i] == gp.abs_(x[i]) for i in range(ncoeffs)))    model.addConstrs((aux[i] == aux[i + 1] for i in range(1, ncoeffs - 1)))    model.setObjective(aux.sum(), GRB.MAXIMIZE)    model.optimize()    return model# k integer > 1# k = 2 -> try coeff = 0.0208333333# k = 3 -> try coeff = 0.003125# k = 4 -> I believe it should be feasible for coeff = 0.0005361519def mwe_signs(k, coeff):    ncoeffs = 4 ** k    first = 1 / ncoeffs    A = reduce(np.kron, repeat(A0, k))    B = reduce(np.kron, repeat(B0, k))    model = gp.Model()    signs = model.addMVar((ncoeffs - 1,), vtype=GRB.BINARY, name="signs")    x = [first] + [(2 * s - 1) * coeff for s in signs]    y = A @ x    model.addConstrs((y[i] >= 0 for i in range(ncoeffs)))    model.addConstr(y.sum() == 1)    tmp = A @ B @ x    model.addConstrs((tmp[i] >= 0 for i in range(ncoeffs)))    model.setObjective(first + (ncoeffs - 1) * coeff, GRB.MAXIMIZE)    model.optimize()    return model

Thanks in advance for any help!