• Gurobi Staff

Hi Michael,

In general, small changes of the model can have a huge effect on the performance, we observe this behavior regularly.

Do you have an idea whether the optimal value for the second model is nearer to the bound of 161 or nearer to the current incumbent of 2055? If you re-compute the objective value of your first solution based on the cost factors, what is the value you get? The question is whether a good solution to the second model is also a good solution to the first model, or whether they have a completely different structure.

If it is nearer to 161, Gurobi might not be able to find good heuristic solutions with the modified objective coefficients. What happens if you use the solution of your first model as a MIP start to your second model?

Usually, it is preferable (for a good performance) to have similar coefficients in the objective function for all involved variables. Coefficients that differ by several orders of magnitudes might lead to numerical issues. Could it be that due to the cost factors the coefficients differ now significantly?

Best regards,
Mario

Hello Mario,

first of, thank you for your timely response.

Do you have an idea whether the optimal value for the second model is nearer to the bound of 161 or nearer to the current incumbent of 2055? If you re-compute the objective value of your first solution based on the cost factors, what is the value you get?

If I would use appropriate cost factors for the given solution of the first model the objective value would be 4853.252$. The bound of 161 in model 2 does not make sense since the fixed cost of opening at least one Distribution Center, which is required to do the subsequent routing, is already 1.500$.

Usually, it is preferable (for a good performance) to have similar coefficients in the objective function for all involved variables. Coefficients that differ by several orders of magnitudes might lead to numerical issues. Could it be that due to the cost factors the coefficients differ now significantly?

The cost for distance/travel time is small in comparison to other cost factors in model 2, whereas in model 1 specifically the coefficient for distance is by far the largest. For reference this is how the coefficients change from model 1 to model 2 (for the same solution):

• Cij (Distance/Transportation Cost): 64895.4 meters in model 1 to 16.2$in model 2 • Tij (Time based Cost): 5513 seconds in model 1 to 22.05$ in model 2
• Fixed Costs for Depots: 3000$same for both models • Variable Depot Cost: 1680$ same for both models
• Fixed Vehicle Cost: 135$same for both models So yes, you are right the coefficients now differ significiantly. If that is the reason I'm running into the performance issues, how can I tackle that problem? What happens if you use the solution of your first model as a MIP start to your second model? I'm not yet familiar with setting MIP starts but will read up and try to do that. Best regards, Michael Renner • Gurobi Staff Hi Michael, If the dual bound of 161 is weak, you might think of improving your model. The fact that you have to open at least one DC could be directly added as constraint, i.e., sum_i y_i >= 1. Additionally, you might force the y variables up by the x variables, e.g., sum_{j \in J} x_{ijk} <= y_i, for all i \in I, k \in K. The coefficient changes might have lead to a shift of focus from solutions with low routing costs (in first model) to solutions with low depot costs (in second model). So, you need to take care that the depot variables are tightly linked in the model. You might also think of additional constraints to push up the z variables? Best regards, Mario Hi Mario, thank you for your suggestions. As you might have noticed I am just beginning to learn about mixed integer linear programming, so I am grateful for the pointers. The fact that you have to open at least one DC could be directly added as constraint, i.e., sum_i y_i >= 1. Additionally, you might force the y variables up by the x variables, e.g., sum_{j \in J} x_{ijk} <= y_i, for all i \in I, k \in K. I have added both those constraints and while the bound is now reasonable the optimization performance is still not improving. The coefficient changes might have lead to a shift of focus from solutions with low routing costs (in first model) to solutions with low depot costs (in second model) That certainly seems to be the case. I'm still struggling to understand as to why the performance seems to be worse in the routing part of the optimization, since model 2 struggles even when only one possible depot is given (so the focus of the optimization should theoretically be on the routing cost?). Could it be that the solver struggles because the routing coefficients are now much closer together (in absolute terms)? The solver finds a reasonable bound and subsequently just stays within a few 0.01$ of that bound.

Do you have any recommendations on how I could deal with these problems, which seem to be mostly numerical?

So, you need to take care that the depot variables are tightly linked in the model. You might also think of additional constraints to push up the z variables?

I'm not sure if I understand what you mean by tightly linked. How would you go about pushing up the z variables?

I'll post the implementation of my constraints just to rule out any mistakes in this step of the optimization.

for j in J:    m.addConstr(quicksum(x[i,j,k] for i in [*I,*J] for k in K if i!=j) == 1)    for k in K:    m.addConstr(quicksum(set_of_all_customers.loc[j].Demand_C * x[i,j,k] for i in [*I,*J] for j in J if i!=j) <= set_of_all_vehicles.loc[k].capacity_V)    for l in J:    for j in J:        if l!=j:            for k in K:                m.addConstr(u[l,k] - u[j,k] + (len(set_of_all_customers) * x[l,j,k]) <= len(set_of_all_customers) -1)                for i in [*I,*J]:    for k in K:        m.addConstr(quicksum(x[i,j,k] for j in [*I,*J] if i!=j) - quicksum(x[j,i,k] for j in [*I,*J] if i!=j) == 0)        for k in K:    m.addConstr(quicksum(x[i,j,k] for i in I for j in J) <= 1)    for i in I:    m.addConstr(quicksum(z[i,j] * set_of_all_customers.loc[j].Demand_C for j in J) - (set_of_all_DC.loc[i].capacity_DC * y[i]) <= 0)        for i in I:    for j in J:        for k in K:            m.addConstr(quicksum(x[i,u,k] + x[u,j,k] for u in [*I,*J]) - z[i,j] <= 1)for j in J:    m.addConstr(quicksum(x[i,j,k] for i in I for k in K) <= y[i])                  m.addConstr(quicksum(y[i] for i in I) >= 1)

Best regards,

Michael Renner

• Gurobi Staff

Hi Michael,

Could you please post the solver log of the model including the newly added constraints? Often this helps to understand what is going on and which Gurobi parameters might be beneficial.

Sometimes, the performance differences of two similar models are not easy to understand, this is a typical phenomenon for modern MIP solvers. The solver might take a completely different path during the optimization, just because of different objective coefficients. Additionally, VRPs are very hard to solve, as you can probably see in the literature.

Regarding tightly linking the z-variables: In an LP relaxation to your model it could happen that many x_ijk variables have, e.g., value 0.5. Then, in constraints (8) it might be feasible to leave the z-variables at zero and still satisfy all the constraints.
The additional constraints I mentioned for linking the y-variables are not necessary for the integer feasibility of your model, but they can strengthen the bounds obtained from the LP relaxation (by avoiding such cases as above) which in turn helps to close the optimality gap.
For the z-variables it is less obvious how to find good additional constraints, here are some examples:

• You know that each customer has to be served from some depot, so "sum_i z_ij = 1, for all j \in J".
• At least for the customers directly connected to a depot in a route, you can set tight links, i.e., "sum_k x_ijk <= z_ij, for all i \in I, j \in J".

Let's see whether those can help.

Mario

Hello Mario,

I revisited the lp file of my model and noticed I made a mistake when adding one of your constraints:

sum_{j \in J} x_{ijk} <= y_i, for all i \in I, k \in K.

After correctly implementing this one and the rest of your recommended constraints the model now runs faster than ever. In fact, the optimization is so quick I'm not sure whether I can trust the solution. I played around with multiple cost factors to see if there is any weird performance impact, and there was none yet.

Following is an instance with 75 Customers and 3 possible Depots, which normally would have taken me upwards of 16 hours to get within this level of optimality (with the faulty distance/time coefficients)

Gurobi 9.1.2 (win64) logging started Wed Nov 3 11:43:14 2021Changed value of parameter LogFile to gurobi_log50.logPrev: Default: Changed value of parameter MIPGap to 0.1Prev: 0.0001 Min: 0.0 Max: inf Default: 0.0001Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)Thread count: 4 physical cores, 8 logical processors, using up to 8 threadsOptimize a model with 17953 rows, 18705 columns and 228651 nonzerosModel fingerprint: 0x6bff7982Variable types: 225 continuous, 18480 integer (18480 binary)Coefficient statistics:Matrix range [1e+00, 2e+02]Objective range [1e+00, 1e+02]Bounds range [1e+00, 1e+00]RHS range [1e+00, 1e+02]Presolve removed 0 rows and 234 columnsPresolve time: 0.25sPresolved: 17953 rows, 18471 columns, 227301 nonzerosVariable types: 225 continuous, 18246 integer (18246 binary)Deterministic concurrent LP optimizer: primal and dual simplexShowing first log only...Concurrent spin time: 0.00sSolved with dual simplexRoot relaxation: objective 4.926199e+02, 644 iterations, 0.11 secondsNodes | Current Node | Objective Bounds | WorkExpl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time0 0 492.61986 0 153 - 492.61986 - - 0s0 0 500.15822 0 180 - 500.15822 - - 3s0 0 501.55832 0 201 - 501.55832 - - 4s0 0 503.07794 0 214 - 503.07794 - - 4s0 0 503.34102 0 196 - 503.34102 - - 4s0 0 503.34102 0 193 - 503.34102 - - 5s0 0 507.61986 0 310 - 507.61986 - - 6s0 0 507.61986 0 264 - 507.61986 - - 6s0 0 507.61986 0 310 - 507.61986 - - 7s0 0 507.61986 0 309 - 507.61986 - - 7s0 0 507.61986 0 213 - 507.61986 - - 9s0 0 507.61986 0 248 - 507.61986 - - 9s0 0 507.61986 0 253 - 507.61986 - - 9s0 0 507.61986 0 263 - 507.61986 - - 9s0 0 507.61986 0 272 - 507.61986 - - 10s0 0 507.61986 0 276 - 507.61986 - - 10s0 0 507.61986 0 182 - 507.61986 - - 11s0 0 507.61986 0 214 - 507.61986 - - 11sH 0 0 635.4093000 507.61986 20.1% - 11s0 0 546.33107 0 182 635.40930 546.33107 14.0% - 12sH 0 0 634.9966400 546.33107 14.0% - 18sH 0 0 625.4411100 546.33107 12.6% - 18s0 0 546.33107 0 238 625.44111 546.33107 12.6% - 18s0 0 546.33107 0 264 625.44111 546.33107 12.6% - 18s0 0 546.33107 0 264 625.44111 546.33107 12.6% - 19s0 0 546.33239 0 189 625.44111 546.33239 12.6% - 19s0 0 546.33239 0 188 625.44111 546.33239 12.6% - 20s0 2 546.59621 0 188 625.44111 546.59621 12.6% - 22s19 22 546.59621 6 255 625.44111 546.59621 12.6% 212 25sH 30 31 613.7123200 546.59621 10.9% 183 27sH 33 31 611.3607300 546.59621 10.6% 172 27s80 95 546.98063 16 294 611.36073 546.59621 10.6% 182 30s399 421 552.40852 65 234 611.36073 546.59621 10.6% 88.2 35sH 471 467 611.0197900 546.59621 10.5% 80.9 49sH 472 467 610.0201500 546.59621 10.4% 80.9 49sH 473 467 609.7646900 546.59621 10.4% 80.7 49sH 474 467 607.6571200 546.59621 10.0% 80.6 49sH 475 471 605.8061600 546.59621 9.77% 80.4 52sCutting planes:Gomory: 18Clique: 5MIR: 119Zero half: 19RLT: 144Explored 479 nodes (52766 simplex iterations) in 52.91 secondsThread count was 8 (of 8 available processors)Solution count 10: 605.806 607.657 609.765 ... 635.409Optimal solution found (tolerance 1.00e-01)Best objective 6.058061600000e+02, best bound 5.465962100000e+02, gap 9.7737%

I manually inspected the solution for the routing and there were no nonsensical connections as far as I can see. Can you see any red flags within the log or does the model just work a lot better now? I will keep on experimenting with different cost factors but hope that this was the solution to my problem.

Thanks to you Mario, I would not have been able to solve this problem myself.

Best Regards,

Michael Renner

• Gurobi Staff

Hi Michael,

Good to see that it works better now!
I cannot see any suspicious lines in the solver log. You might even set a smaller MIPGap to get better solutions.

I would recommend to write a method that checks your final solution for feasibility, e.g., by building a small graph structure to see whether all routes are connected to some depot, whether the capacity is satisfied for each vehicle, etc. This helps finding modeling mistakes.

Best regards,
Mario

Hi Mario,

I will now run the automatic tuning tool again to see from which parameters i can benefit to find solutions for different model sizes, as well as what a reasonable MIPGap is.

I already have a few different solution printers in place, which let me manually check the feasibility of the provided solution. I will get back to you if something should come up within the following experiments. If not, I sincerely thank you for your time and effort.

Best Regards,

Michael Renner