Solve MILP using Big M Method (Yalmip + Gurobi)

[Zhe Yu]

Hello Professor Löfberg,

I am solving a MILP using big M method. I have troubles getting the right values for big M to fast reduce the MIP gap. Currently, what I am doing is by trial and error, which is like a black box, very unreliable.

On the other hand, if not using big M method, there will be bilinear terms (binary variables multiplied by continuous variables) in the objective function. I tried to solve that using CPLEX, it is not tractable as well.

Thank you very much if you could offer some help!

Zhe Yu

[Johan Löfberg]

If you don't have any good a-priori knowlege about bounds, or can solve simple optimization problems to estimate them, you will get a crappy model. That's just a fact. However, you have bounds around 1. Have you really proven them to be sufficiently large. Looks suspiciously small.

BTW, you should vectorize your code

[Johan Löfberg]

btw, is this some kind of manual attempt at modelling an absolute value?

[Johan Löfberg]

this does not model psi_w_plus=z_w_plus*pi as your comments in code indicate they should

        constraint = [constraint, psi_w_plus(i,t) >= -BM_pi*z_w_plus(i,t)];

        constraint = [constraint, psi_w_plus(i,t) >= pi(i,t) - BM_pi*(1-z_w_plus(i,t))];

z_w_plus = 0 leads to psi_w_plus>=0, psi_w_plus >= pi

[Zhe Yu]

Thank you for your quick reply.

In my objective function, it is max - w_plus(i,t)*psi_w_plus(i,t), which is equivalent to min w_plus(i,t)*psi_w_plus(i,t).

[Zhe Yu]

That is just my guessing, I am assuming that sufficiently large M is right. 

Budget_w = sum(out_z_w_plus + out_z_w_minus,2) == Gamma_w;

Budget_l = sum(out_z_l_plus + out_z_l_minus,2) == Gamma_l;

When M = 10000, it is fast but Budget_w  and Budget_l are not all true, which I assume they will be true at optimal solutions.

[Zhe Yu]

I think z_w_plus = 0 leads to psi_w_plus>=0, psi_w_plus >= pi - BM_pi

psi_w_plus will be 0 since it is minimizing w_plus(i,t)*psi_w_plus(i,t)

[Johan Löfberg]

'gurobi.TimeLimit','20' should be 'gurobi.TimeLimit',20 ('20' works, but that is not official behavior)

What I am confused about is that you model x*y for binary y with only 1 constraint. That would require 2 constraints in most situations. Perhaps you have a model where 1 of the constraints will be satisfied by optimality, and you've proven that, but I would definitely keep the full model to give the MILP solver as much information as possible (redundant constraints are good for MILP solvers, as they might lead to improved relaxations)

and of course, you can simplify your code tremendously by dumping all the indexing

So, if you really want psi_w_plus to model pi*z_w_plus, the model would normally be

constraint = [constraint, BM_pi*z_w_plus >= psi_w_plus >= -BM_pi*z_w_plus];

constraint = [constraint, BM_pi*(1-z_w_plus) >= psi_w_plus-pi >= - BM_pi*(1-z_w_plus)];

constraint = [constraint, -BM_pi*z_w_minus <= psi_w_minus <= BM_pi*z_w_minus];

constraint = [constraint, -BM_pi*(1-z_w_minus) <= psi_w_minus-pi <= BM_pi*(1-z_w_minus)];

However, I am still wondering what it is you actually try to model, i.e., what kind of piecewise affine/quadratic structure is it you minimize, I'm slightly worried that you actually have a piecewise convex quadratic that you waste binary variables on

[Johan Löfberg]

[Zhe Yu]

Please find attached for part of the objective function. Basically, there are two bilinear terms I want to linearize: w * pi and L * delta (pi >= 0, delta is unrestricted). Parameters w and L are represented by uncertainty sets.

[Johan Löfberg]

Sure, this is what you are implementing, but I was more wondering if these are low-level representation of something which can be expressed more clearly on a high level, such as a picewise affine structure of some type or something, involving max/sign/etc/switches

[Zhe Yu]

Attached is the model I am trying to solve. All the constraints are linear except for the uncertainty sets (binary constraints). The parameters under max are random and represented by uncertainty sets. Since the third stage objective function (last line of (19a)) has no recourse decisions, I can dualize it and then combine with the second line objective function under min. Next, dualize the second min to max and combine it with the first max, where I got those bilinear terms (uncertain parameters multiplied by dual variables).

[Zhe Yu]

Hi Professor Löfberg, how could I vectorize when sdpvar and nsdpvar both in the same constraints? Also, do all the constraints can be vectorized?

[Johan Löfberg]

You don't have any ndsdpvar so I don't see the relevance, but if you had, the fact that you have sdpvar and ndsdpvar does not change anything (besides  possibly making the vectorization trickier as you have to juggle one more index in your head)

All your current constraints are trivially vectorizable by simply skipping the loops and indexation, except for some where you might have to do some repmat, fixing summation direction, or shifting data as you have both t and t-1.

[Zhe Yu]

Thank you very much!