Bug in barrier solver with Gurobi 6.5.2 for QP or QCP model when using piecewise linear objective functions via PWLObj

[Viktor Drnic]

There appears to be an issue with the barrier solver when using the PWLObj functionality.

A pair of reproductions of the issue is attached.

Model A:

Model A0 and model A1 logically describe the exact same QP model, the only difference is that A1 uses PWLObj, while A0 does not.  However, Gurobi returns -15.28 for the objective function for A0, and -12.54 for A1 (See the attached files model-a[0|1].[lp|sol])

Manually computing the value of the objective function for the solutions provided show that -15.28 is indeed the correct value for A0, while -16.63 is the true value of the objective function at the solution provided by Gurobi for model A1, not -12.54 as claimed.

To see that it is the combination of the barrier solver and the PWLObj which is causing the error, note that:

• If the quadratic term from the objective function is removed ([ - 20 xp_1 ^2 - 40 xp_1 * xn_1 - 20 xn_1 ^2  - 2 dp_1 ^2 ] / 2) from both models, then the LP solver is used instead of the barrier solver, and both modified models give the same answer (-15.28) (See the attached files model-a[0|1]-noq.[lp|sol]).  Note that this is the same value of the objective function above, as it turns out that the quadratic term is zero at the solution, although the values of the variables at the solution are different.
  
• Alternatively, if the PWLObj terms are removed instead of the quadratic terms, the barrier solver gives the same answer for both modified models (-13.49) (See the attached files model-a[0|1]-nopwl.[lp|sol])

Model B:

Model B0 and B1 logically describe the same QCP model, again the only difference is that B1 uses PWLObj, while B0 does not.  However, Gurobi returns 0.134 for the objective function for B0, and 0.321 for B1.

Manually computing the value of the objective function for the solutions provided show that 0.134 is indeed the correct value for B0, while -1.092 is the true value of the objective function at the solution for model B1, not 0.321 as claimed.

In fact, the sum of all terms in the objective function except for the piecewise linear terms at the solution for B1 is 0.321, which suggests that the piecewise linear terms are somehow being left out of the optimization.  This is despite the fact that the Gurobi console output states "Model has 58 piecewise-linear objective terms" when loading the model.

If the quadratic constraint is removed from both models, then the LP solver is used and both modified models return 0.134 for the final value of the objective function, indicating again that the issue lies with the barrier solver.

If the piecewise linear terms are removed from both models, then both modified models return 0.321 for the final value of the objective function - which matches the value returned for unmodified B1 model, again indicating that the piecewise linear terms are not being taken into account in the optimization by the barrier algorithm for the unmodified model B1.

[Sonja Mars]

Hi,

Thanks a lot for reporting this. It helped us to identify an issue in Gurobi's presolving routines. Our developers are currently working on a fix. We will get back to you soon.

Thanks and best regards,
Sonja

[Sonja Mars]

Hi,

Our developers have fixed two issues in Gurobi's presolving routines. This fix will be part of our next release.

However, please note that your models with and without PWL parts are equivalent. For example in model-a0.lp you use y_2 to express two different piecewise linear terms.


Thanks again for reporting this,
Sonja

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "Gurobi Optimization" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to gurobi+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Viktor Drnic]

Thanks Sonja.  Out of interest - is there any performance benefit to be obtained using the PWLObj functionality vs manually constructing an equivalent set of constraints (e.g., as is done in model-a0.lp)?  That is, does the underlying implementation use any special optimizations for PWLObj terms that would not otherwise occur?

[Viktor Drnic]

There appears to be another issue with the status code returned by Gurobi in some cases.

Running the attached model with PreSolve = 1 returns a status code of 2 (optimal), while the console output states "Warning: max constraint violation (5.0000e-02) exceeds tolerance  (model may be infeasible or unbounded - try turning presolve off)".  In such a situation, I would expect that Gurobi would return a status code of 4 (INF_OR_UNBD), not 2.  Setting PreSolve=0 and re-optimising gives a very different result without the warning about constraint violations.

Full log of interactive session below:

Python 2.7.8 (default, Jun 30 2014, 16:08:48) [MSC v.1500 64 bit (AMD64)] on win32

Type "help", "copyright", "credits" or "license" for more information.

Gurobi Interactive Shell (win64), Version 6.5.2

Copyright (c) 2016, Gurobi Optimization, Inc.

Type "help()" for help

gurobi> m=read('model.lp')

gurobi> m.setParam('PreSolve', 1)

Changed value of parameter Presolve to 1

   Prev: -1  Min: -1  Max: 2  Default: -1

gurobi> m.optimize()

Optimize a model with 30 rows, 25 columns and 106 nonzeros

Model has 36 quadratic objective terms

Coefficient statistics:

  Matrix range    [3e-02, 5e+01]

  Objective range [2e-01, 6e+01]

  Bounds range    [1e-01, 1e+01]

  RHS range       [6e-07, 9e+00]

Found heuristic solution: objective -1e+10

Presolve removed 11 rows and 7 columns

Presolve time: 0.00s

Presolved: 19 rows, 18 columns, 66 nonzeros

Presolved model has 36 quadratic objective terms

Variable types: 15 continuous, 3 integer (2 binary)

Root relaxation: objective 3.446219e+00, 14 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       3.4462190    3.44622  0.00%     -    0s

Explored 0 nodes (14 simplex iterations) in 0.08 seconds

Thread count was 24 (of 24 available processors)

Optimal solution found (tolerance 1.00e-04)

Warning: max constraint violation (5.0000e-02) exceeds tolerance

         (model may be infeasible or unbounded - try turning presolve off)

Best objective 3.446219034861e+00, best bound 3.446219034861e+00, gap 0.0%

gurobi> m.Status

2

gurobi> m.setParam('PreSolve', 0)

Changed value of parameter Presolve to 0

   Prev: 1  Min: -1  Max: 2  Default: -1

gurobi> m.reset()

gurobi> m.optimize()

Optimize a model with 30 rows, 25 columns and 106 nonzeros

Model has 36 quadratic objective terms

Coefficient statistics:

  Matrix range    [3e-02, 5e+01]

  Objective range [2e-01, 6e+01]

  Bounds range    [1e-01, 1e+01]

  RHS range       [6e-07, 9e+00]

Found heuristic solution: objective -1e+10

Variable types: 21 continuous, 4 integer (1 binary)

Root relaxation: objective 2.579814e+00, 44 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    2.57981    0    1 -1.000e+10    2.57981   100%     -    0s

H    0     0                       2.5462190    2.57981  1.32%     -    0s

Cutting planes:

  MIR: 1

Explored 0 nodes (44 simplex iterations) in 0.08 seconds

Thread count was 24 (of 24 available processors)

Optimal solution found (tolerance 1.00e-04)

Best objective 2.546219034859e+00, best bound 2.546228467961e+00, gap 0.0004%

gurobi> m.Status

2

[Sonja Mars]

Hi,

> Thanks Sonja. Out of interest - is there any performance benefit to be obtained using the PWLObj functionality vs manually constructing an equivalent set of constraints (e.g., as is done in model-a0.lp)? That is, does the underlying implementation use any special optimizations for PWLObj terms that would not otherwise occur?

Piecewise-linear objectives are handled directly by the simplex algorithm, so if you have a continuous model with a convex, piecewise-linear objective, you should see a significant performance boost from using this feature.

Best regards,
Sonja


----
Dr. Sonja Mars
Gurobi Optimization - Technical Support

[Sonja Mars]

Hi Viktor,

Thank you for reporting this second issue. It helped our developers to identify a bug in our QP Simplex algorithm. This is why you get this solution with the big violation ( "Warning: max constraint violation (5.0000e-02) exceeds tolerance."). Using the solution quality attributes you can easily check for violations and find out if you can trust the solution. You can find these attributes on this page: http://www.gurobi.com/documentation/6.5/refman/attributes.html in the paragraph "Solution quality attributes:".

The fix will be part of our next release.

Thanks again,
Sonja

[Viktor Drnic]

Hi Sonja,

I tried the barrier solver with PWL terms again with Gurobi 7.0.1, however there still appears to be an issue.

The two attached models (model0 and model1) logically describe the exactly same problem, however model0 uses multiple constraints to construct the piece-wise linear objective terms, while model1 uses the PWLObj functionality of Gurobi.

I believe that model0 is giving the correct solution, while model1 is behaving very strangely.  In particular:

• If the quadratic constraint is removed from both models (i.e., the constraint called 'variance'), both models give exactly the same answer via the LP solver, suggesting that the issue is with the barrier solver.
• If you take the values of the variables at the solution Gurobi provides for model1 (see model1.sol) and manually calculate what the value of the objective function should be at that point, it comes out to 21.2.  Meanwhile, Gurobi returns 28.7 for the ObjVal property of the model, which is a very large difference.  Oddly, this is closer to the expected solution given by model0 of 28.3.
• The difference between the two solutions is due to the values of the variables tp_0, tn_0, tp_1, and tn_1 at the solutions.  Due to the nature of the model, at the correct solution only one of tp_0 and tn_0 should be nonzero, and similarly only one of tp_1 and tn_1 should be nonzero - as is the case with the solution to model0.   However, in the solution provided by Gurobi for model1 all these values are significantly different to zero.

Regards,

Viktor

[Sonja Mars]

Hi Victor,

Do you also have log files showing this issue? That would be very helpful.

Thanks,
Sonja

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "Gurobi Optimization" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to gurobi+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <model0.sol><model1.lp><model1.sol><gurobi.env><model0.lp>

[Viktor Drnic]

Sure, see attached

[Sonja Mars]

Hi Victor,

Thanks for sending the log and for reporting this. Our developers are currently investigating this.

Thanks,
Sonja

> <gurobi.log>

[Sonja Mars]

Hi Victor,

Thanks a lot for reporting this. Our developers identified an issue and fixed it for our next release.

As a workaround you should use Method=2 whenever a piecewise linear model has a quadratic objective or quadratic constraints.

Thanks again and best regards,

Sonja


----
Dr. Sonja Mars
Gurobi Optimization - Technical Support

[Viktor Drnic]

Thanks Sonja, indeed the workaround seems to do the trick.

[Viktor Drnic]

Hi Sonja, I think there may be another issue with either PWLObj or some change between Gurobi 6.5.2 and 7.0.1.  If you try and run a bunch of QP models with PWLObj objective functions in multiple threads (each solver is single threaded), then the CPU utilization is very low (a few percent).  The Visual Studio profiler seems to indicate that this is heap contention (RtlpAllocateHeapInternal and RtlFreeHeap in ntdll.dll), although I'm not sure.  Doing the exact same thing with 6.5.2 without PWLObj, you will get near 100% CPU utilization.  I can put together a reproduction of the issue if you like, although it may not be until the new year.  Rgds, V.

[Sonja Mars]

Hi Viktor,

I have never seen such a behavior before. It would be great if you can give us some small example code to reproduce the issue. We are also happy to take a look at it next year, no need to rush.

Thanks,
Sonja