pabloazurduy / qp-mip

a quadratic linealization example, using python-mip library and coinOR solver
0 stars 0 forks source link

Instability to approximation parameters #2

Open MaxMayya opened 3 years ago

MaxMayya commented 3 years ago

Thank you for sharing this work. I tried to use get_quadratic_appx to approximate the following quadratic integer problem:

argmin(j) ||u(theta + 2pi * j) - r||ˆ2

where u (3d unit vector) and theta are angle-axis representation of a rotation and r (3d vector) is a groundtruth rotation (with no angle ambiguity). The idea is to pick the best j so that u|theta is as-close-as possible to r. when I change min_interval, max_interval or num_linspace, results change drastically!

here is a minimal code:

        mipModel = mip.Model()
        j = mipModel.add_var(var_type=mip.INTEGER)
        vars = u * (theta + 2 * np.pi * j) - r
        qaVars = {}
        for v, var in enumerate(vars):
            qaVars[v] = get_quadratic_appx(model=mipModel, var=var, min_interval=3, max_interval=3, num_linspace=50)
        mipModel.objective = mip.minimize(mip.xsum(qaVars[v] for v in range(len(vars))))
        status = mipModel.optimize()

I'm I missing something?

Many thanks

pabloazurduy commented 3 years ago

Hi @MaxMayya was looking at your particular example, a few comments from what I noticed overviewing your code:

  1. min_interval and max_interval are the "linearization bounds" that means that given that we want to linearize x^2 it must contain the domine of x = f(y), in your case x = u(theta + 2pi * j) - r therefore we want to study x \in [min_interval, max_interval] (they shouldn’t be equal, min_interval< max_interval )
  2. given that (from what I'm seeing) your problem haven’t linear constraints or other linear definitions, its way more efficient to use a non linear solver (there are plenty of them available for python). I personally recommend the scipy implementation, it even allows you to constraints some variables.
  3. The implementation of this function relies on a number of discrete intervals (num_linspace) this increments a lot the cost of optimization every time you increment that value, the trade-off is that it would be more imprecise if you use fewer intervals (for example 10)

Some additional notes:

(about 2). ** the reason why we might be interested in modeling a QP in a linear solver is that you might have a big optimization linear-problem, and for some reason, you need to add a few quadratic functions. If that's not your case I would recommend choosing a non-linear solver.

(about 3). there is a more efficient way to linearize a function implemented in CBC, called SOS I' didn’t use when I wrote this code because there was a bug on the SOS-CBC /python-mip integration, which is now solved. That implementation should be working now, so it might be productive to use the SOS linearization

Thanks ! Hope it helps, let me know If you need anything else

MaxMayya commented 3 years ago

Man! Thanks a lot for elaborating. Sure x is not within the approximation min/max limits. Shame on me! Following your recommendation, I looked for other python interfaces and found pyomo which provides interface to multiple optimizers (including NL ones like IPOPT). Seems promising, I'll give it a shot. Thanks again