Originally posted by: EhsanN

Warning: Control callbacks may disable some MIP features.

MIP emphasis: balance optimality and feasibility.

MIP search method: traditional branch-and-cut.

Parallel mode: none, using 1 thread.

Root relaxation solution time = 9.31 sec. (6789.24 ticks)

 

        Nodes                                         Cuts/

   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap         Variable B NodeID Parent  Depth

 

      0     0   208329.2048     2                 165164.8395    45047         

      0     2   208329.2048     2                 208329.4789    45047                                 0             0

      1     3   208567.0599     3                 208456.5453    46615                    Z(19) D      1      0      1

...

     27    29   249364.9619     1                 213933.3465   229849                     Z(6) D     27     25     14

*    28    24      integral     0   252172.4612   213933.3465   230821   15.16%           Z(24) U     28     27     15

     29    25   250764.3937     2   252172.4612   213933.3465   234670   15.16%           Z(24) D     29     27     15      

...            

     90    60   228195.5669     5   252172.4612   217436.3989   584301   13.77%            Z(3) D     90     89      8

     91    61   228322.1725     4   252172.4612   217436.3989   584694   13.77%           Z(14) D     91     90      9

     92    62   228797.2027     4   252172.4612   217436.3989   586584   13.77%           Z(27) D     92     91     10

*    93     0      integral     0   -44112.5927   217436.3989   592281 -592.91%            Z(1) D     93     92     11

 

User cuts applied:  152

 

Root node processing (before b&c):

  Real time             =  147.39 sec. (110701.04 ticks)

Sequential b&c:

  Real time             = 43602.99 sec. (48012402.99 ticks)

                          ------------

Total (root+branch&cut) = 43750.38 sec. (48123104.03 ticks)

I'll assume that the [0, 1] are explicitly stated, not just implied. My first suggestion is to turn on the numerical emphasis parameter and/or collect kappa statistics, to see if you have numerical stability issues.

Paul

Originally posted by: EhsanN

Thanks a lot, Paul. It seems there is a numerical instability issue with that specific instance and turning on numerical emphasis would solve the problem. Although, the surprising thing is that the particular instance facing the error is very similar to more than 15 other instances, which are solved perfectly fine. Their difference is just two parameters and there are instances, where those two parameters are larger two or three times, I guess you'll never know it with MIPs for sure.

PS. I turned on the MIP Kappa parameter and set the TuningDisplay parameter to show full information, but the log file didn't show anything on Kappa statistics. What is the correct way to obtain this information?

Ehsan,

I think numerical instability is very hard to predict in advance. The only generalization I can think of is that people who have 12 or more orders of magnitude difference in their constraint coefficients are skating on very thin ice. Other than that, I guess you just have to wait and see if it crops up. In particular, a modest change in some coefficient can shift that path of branch-and-bound search so that an unstable basis that was lurking but never encountered before pops up (or vice versa).

The kappa statistics are part of the "solution quality" metrics (when collected). In the interactive optimizer, just display solution quality and I think you will get them. In the programming APIs, you have to fetch the solution quality results (I'm pretty sure there is a specific class for them) and then drill down. I don't think tuning display has anything to do with them.

Originally posted by: EhsanN

Dear Paul,

Thanks again. I just found the necessary method in C++, which is getQuality. Additional information is available from here.

I would rerun that instance again and post here additional information if something interesting comes up in the Kappa statistics.

Bests,

Ehsan

Originally posted by: EhsanN

I just got a chance to follow Paul's advice and get the Kappa stats for my model. The reports for the model with both numerical emphasis on and off is included below. As I've no prior experience with Kappa stats, I'm not sure how to interpret them. In particular, values of ExcatKappa and Kappa are printed as -1, which seems like an error.

Kappa Stats (with numerical emphasis off)

Kappa Stats (with numerical emphasis on)

PS. For other instances that shows no strange behavior, KappaStable is 1, while other measures are 0 and KappaMax is significantly smaller (like 100 times).

I'm used to the diagnostics printed by the interactive optimizer, rather than the version returned by API calls, so I'm not entirely sure what "Kappa" and "ExactKappa" are. My best guess is that they are stats for a single basis, and the -1 just means that there is no value for them because you're looking at stats for all bases, not the condition number of a single basis.

The Stable/Suspicious/Unstable/Ill-posed results (in increasing order of cause for alarm) look good -- 99% of bases are stable, none are ill-posed -- and the KappaAttention value (which, per the user manual, estimates the likelihood of there being problems) is quite small. So the takeaway is that whatever is going on in your model is likely not the result of numerical instability (although we cannot rule it out with absolute certainty).

One wrinkle in this is that the kappa stats come from the presolved version of the model. I think in some cases it is possible that the transformation from the original model to the presolved version (or perhaps the reverse transformation) introduces some major rounding error. Also, there have been occasional bugs in the presolve routine that, on certain "edge" cases, produce incorrect results. You might want to try the problem on CPLEX 12.6 and see if it behaves differently.

Originally posted by: EhsanN

Thanks, Paul. I just tried solving the instance with CPLEX 12.6,. The results and Kappa stats have not changed except for the KappaMax measure, which has increased about two time when the numerical emphasis option is turned on.

Kappa Stats (with numerical emphasis on)

Since you are already using callbacks, I suggest you monitor new incumbents in an incumbent callback. When the first nonsense solutions shows up, call getSolutionSource() to find out if it is a heuristic solution or a node solution. (I think this is new in 12.6.) Also call getQuality() (several times) to get SumIntInfeas, SumPrimalInfeas, SumScaledPrimalInfeas and any other measures you think might help determine whether CPLEX actually believes the solution is feasible.

If it's a bug in CPLEX, you'll probably need to pass a copy of the model and code to the IBM folks to let them diagnose it.

Originally posted by: EhsanN

I'm speculating here, but it's possible (if not likely) that the quality measures you get in the callback are based on the presolved problem, whereas what you get from getQuality() after the solver returns a "solution" is based on the original model.

Someone in a different thread ran into some dopeyness that was related to the presolver (and what I suspect is a bug in it). You might try turning off presolving entirely. If that cures the problem, turn it back on and try selectively disabling bits of it (primal reductions, dual reductions, ...) until you find which part of presolving is creating the issue. If turning off the master presolve switch does not fix the problem, then you'll need to look elsewhere for the root cause.

Paul

Originally posted by: EhsanN

I've had disabled CPLEX presolve (IloCplex::PreInd) from the beginning of my research (based on previous experience, presolve has little advantage for my model). In the meanwhile, I have tried many combinations of CPLEX presolve-related parameters (including primal/dual reduction and node presolve), but I had no success. I guess, as you've mentioned, there is a bug somewhere else (either in my code or CPLEX). Since, I've a solution for the faulty instance with numerical emphasis turned on, I'm going to move on from this instance. The sad part is that I've to recheck all other instances to make sure nothing is wrong with them. :-(

Anyway, thank again for time and kind attention. You've been a great help as always.

Originally posted by: EdKlotz

Ehsan, can you upload a SAV file of the model that gives you the negative MIP gap?

Originally posted by: EhsanN

Ed: I'm adding the lazy constraints through a callback and a separation method. So, I don't know them in advance. As far as I know, lazy constraints are not included in exported models.

When I tried to export the model using cplex.exportModel() statement before and after the solve statement, the output of the second one was the same as the first one (I solved them using the interactive optimizer). If you let me know how I could export the model with all the lazy constraints added during the solution process, I will be more than happy to upload the resulting file.

Originally posted by: EdKlotz

Ehsan,

              Sorry, I jumped in at the end of the thread and assumed you just were seeing a negative MIP gap on a regular model.    You are correct that when you use exportModel(), CPLEX does not include the lazy constraints that you add dynamically via callbacks during the course of the optimization.   However, you could accumulate those lazy constraints in a separate IloConstraint or IloRange array within the callback that adds them.   Then, at the end of the optimization,  you could add the lazy constraints to the model and export it.

              It's interesting that the negative MIP gap comes from a solution found via branching rather than heuristics.  Given that turning on numerical emphasis causes the problematic behavior to disappear, the issue could be numeric.   But, it might also have just changed the path CPLEX took and caused the real source of the problem to no longer execute.   So, at a minimum, could you please upload the SAV file of the model at the start of the optimization, or post the results of the 'display problem stats' command after reading it into interactive CPLEX?

              It sounds like you have moved on since turning on numerical emphasis resolved the issue.   If you wanted to investigate this further, I would suggest trying to simplify your program as much as possible in a manner that preserves the unexpected behavior.    For example, you could start by disabling the branch and incumbent callback function and see if the behavior persists.    Also, you could try accumulating the lazy constraints as I described above, then put them in an LP file as lazy constraints at the start of the optimization.

              If you upload a SAV file of either the original model or one with the lazy constraints as well, I will have a look and try to assess the likelihood that  turning on numerical emphasis robustly addressed the issue.

Originally posted by: EhsanN

Regarding running the code with currently found lazy constrains, I've tried the iloCplex.importModel() method with no success. The piece of code I've used is as follows:

//The IloCplex object extract the original model.
cplex = IloCplex(myModel); 
//A new model that will contain the lazy cuts.
IloModel lazyModel(myMainEnv); IloObjective lazyObj; 
IloNumVarArray var(myMainEnv);
IloRangeArray rng(myMainEnv);
cplex.importModel(lazyModel, lazyCutsFileAddress, lazyObj, lazyVar, lazyRng);
cplex.addLazyConstraints(lazyRng);

This code results in an Invalid Cut Exception. I would really appreciate it if you could provide me with some reference or sample code on how to add previously-saved constraints in a SAV file as lazy constraints to my model.

Finally, I've put the SAV files for both the initial and final models on Dropbox (you can access them from here and here). I've also included their problem stats below. I must mention that when I solve the Final Model from the scratch (without adding any other lazy cuts as you asked), it results in a normal termination with a feasible solution.

CPLEX> Problem name: InitialModel.sav

Variables            :  810032  [Nneg: 1,  Fix: 1,  Box: 810000,  Binary: 30]

Objective nonzeros   :  810002

Linear constraints   :   27900  [Less: 27000,  Equal: 900]

  Nonzeros           : 2430000

  RHS nonzeros       :     900

 

Variables            : Min LB: 0.0000000        Max UB: 1.000000       

Objective nonzeros   : Min   : 0.01249680       Max   : 31993.53       

Linear constraints   :

  Nonzeros           : Min   : 1.000000         Max   : 1.000000       

  RHS nonzeros       : Min   : 1.000000         Max   : 1.000000     

CPLEX> Problem name: FinalModel.sav

Variables            :  810032  [Nneg: 1,  Fix: 1,  Box: 810000,  Binary: 30]

Objective nonzeros   :  810002

Linear constraints   :   28052  [Less: 27000,  Greater: 152,  Equal: 900]

  Nonzeros           : 54183591

  RHS nonzeros       :     900

 

Variables            : Min LB: 0.0000000        Max UB: 1.000000       

Objective nonzeros   : Min   : 0.01249680       Max   : 31993.53       

Linear constraints   :

  Nonzeros           : Min   : 0.0005661833     Max   : 19805.08       

  RHS nonzeros       : Min   : 1.000000         Max   : 1.000000

Originally posted by: EdKlotz

Certainly the range of coefficients in the lazy constraints is moderately large, and could potentially cause numerical issues.   However, it is not so bad that we should just assume that.   Still, if you know of a simple way to rescale the coefficients of the lazy constraints so that the range of coefficients is smaller, you might as well do that.

Regarding the Invalid Cut Exception you mention above, I think the problem is that the variables you create for the Lazy Constraints, while they are the same in the SAV, MPS or LP file you imported, are not the same in your C++ program.  So, CPLEX thinks the lazy constraints you are adding involve additional variables that are not part of the model to which you want to add cuts.   Adding cuts with variables that were not previously in the model is indeed an invalid type of cut.  The variables in the lazy constraint model needed to be mapped to the corresponding variables in the original model, and the lazyRng constraints need to reference those original variables.   Rather than trying to do that, I recommend that you create an LP file with the lazy constraints and then proceed.

Also, when I looked in the documentation, it indicated that lazy constraints and user cuts added by your program will appear in any subsequently exported SAV or LP file.

I just downloaded the two zip files and will kick of some runs.   For the one with the lazy constraints, can you give me the index number of name of the first lazy constraint (assuming they all follow the regular constraints).

Originally posted by: EdKlotz

I was able to determine the boundary between actual and lazy constraints and have kicked off some runs.   My first run did not reproduce the negative MIP gap reported above.   The log you attached implies that the problematic run has one thread, traditional branch and cut search, and primal only presolve reductions.    Please let me know if any other non default parameter settings are involved.

Originally posted by: EhsanN

Dear Ed, thanks for your feedback. My answers to your questions are as follows:

1. Regarding rescaling of coefficients, I'll look into it more. However, I think there is little chance to rescale them as the lazy cuts coefficients are results of multiplying two or three other parameters with very different scales. 
2. I'll follow your suggestion about how to read cuts as a lp file and post the results when I get them.
3. I use the exportModel() method to export the model as a SAV file. However, exporting the model before and after the solve statement results in the same exported models. I think that's because I'm using a lazy constraint callback, not addLazyConstraint() method. Can you refer me to the specific part of the documentation on this issue?
4. As I had no idea how to read a SAV file to determine the first lazy cut ID, I converted the SAV file into a LP file. However it was to large to read by a text editor (the LP file is 2.8 GB). I've upload a SAV file containing just the lazy cuts (which is available from here). In addition, the initial model has about 27900 constraints. So if CPLEX keeps the order of constraints, the lazy cuts should be after them. In addition, lazy cuts share a decision variable named "yMinMax", that appears only in them and the objective function.
5. Since, I'm using a lazy constraint callback, I've to use one thread and dynamic search is disabled by default. Other parameter settings are as follows:

cplex.setParam(IloCplex::AdvInd, 0);
cplex.setParam(IloCplex::EpAGap, 0.00001);
cplex.setParam(IloCplex::EpInt, 0.0000001);
cplex.setParam(IloCplex::EpRHS, 0.0000001);
cplex.setParam(IloCplex::PreInd, IloFalse);//The problem happens with IloTrue as well.
cplex.setParam(IloCplex::PreLinear, 0);
cplex.setParam(IloCplex::Reduce, 0);
cplex.setParam(IloCplex::CutsFactor, 1);
cplex.setParam(IloCplex::HeurFreq, 0);
cplex.setParam(IloCplex::RepeatPresolve, 0);

Originally posted by: Laci Ladanyi

> Since, I'm using a lazy constraint callback, I've to use one thread and dynamic search is disabled by default.

While it is true that dynamic search will be disabled if you use a control callback, you are *not* restricted to using one thread. CPLEX will by default use only one thread if a control callback is present, but only beacause CPLEX has no way of checking that the code in the callback is thread-safe or not (after all, if multiple threads are used it is quite possible that sometime more then one thread will execute the callback code simultaneously). If you, the programmer of the control callback code, are sure that your code is thread-safe then just set the number of threads to be whatever you desire (the number of cores in the machine is recommended) and cplex will cheerfully use that many threads.

BTW, if you can enumerate all possible lazy constraints in advance (and it is not some ginormous list) then you may be better off using addLazyConstraint() before optimization is invoked. If you do that then CPLEX will automatically scan that list for violated constraints -- this happens at the same time when the callback is invoked. However, with the list (as opposed to the callback) dynamic search is not disabled and CPLEX automatically runs in parallel. Of course, it is always a tradeoff wich method is better. You may have a very fast algorithm to decide whether there are any violated lazy constraints while CPLEX has to check everything on the list one-by-one -- in this case the callback may be better despite dynamic search being disabled. One has to try... :-).

--Laci

Originally posted by: EhsanN

Dear Laci, thanks for your feedback. For the sake of comparisons, I should use one thread. However, as I was interested in seeing a possible improvement in CPU time, I tried it. For one implementation of lazy constraint callback, I couldn't use it as my separation problem is not thread-safe (I've to access an external data structure). For another implementation, it worked with no error. However, the CPU time increase drastically (from 5 seconds to  about 130 seconds). In particular, the process of adding cuts took a long time.

Regarding possibility of using a pool of lazy cuts and adding cuts via addLazyCosntraint(), I cannot use this approach as the number of candidate lazy constraints is of exponential order (for my smallest instance with just 10 nodes, there are about 2.99 * 10^26 possible cuts). Therefore, I've to separate them using a lazy constraint callback.

BTW, my understanding of thread-safe code is based on the following information from the documentation:

To make sure of that condition, the control callback must use only information queried from CPLEX® itself within the callback as the basis for algorithmic decisions. In other words, no information that accumulated in an external data structure over several invocations of the control callback can be used.

In my callback that doesn't work, I use some global data structures to keep track of found solutions and best one so far. Is this the reason why this callback doesn't work? If yes, how can I make it thread-safe. Is there any other guideline or example on how to make the code thread-safe?

Originally posted by: Laci Ladanyi

Indeed, that many constraints would be a bit difficult to enumerate :-).

One way to make things thread-safe is what is in the documentation :-). Note the phrase "no information that accumulated in an external data structure over several invocations of the control callback can be used". This does *not* mean that you cannot use external info, just that you cannot use info that is accumulated over the course of the solve process. That is, if you set up some data before you invoke cplex to solve the instance, and then you access this data (which stays constant while you are solving the problem!) from the callback then you are still thread safe.

Now, you say you keep track of found solutions. That means you are writing into the global data. If two threads would do this simultaneously then you'd get messed up solutions... What you can do is you can use mutexes to temporarily lock the global data to a thread while that thread writes into the data. (A good tutorial for using mutexes is at http://www.cs.kent.edu/~ruttan/sysprog/lectures/multi-thread/multi-thread.html.) Of course, since there is no guarantee about the scheduling of threads, simply using mutexes will most likely result in non-determinism, i.e., running the code twice with the same input on the same machine will result in your callback in different threads writing into the global data in different order in the two runs. And both runs will be correct :-). Making multithreading deterministic is difficult in certain situations. Still, if you only write into the global data but do not use its content for generating the lazy constraints then at least the CPLEX run will be deterministic.

Finally, about the CPU time increase... I would be very interested to know more. If it's really spinning while adding cuts that is bad. On the other hand, it's possible that you are just very unlucky and by increasing the number of threads the search tree has changed in a way that just increased the size of the tree by a factor of 30. Could you post the logs of the single threaded and the multithreaded runs? Or just how many nodes were evaluated in each?

--Laci

Originally posted by: EhsanN

MIP emphasis: balance optimality and feasibility.

MIP search method: traditional branch-and-cut.

Parallel mode: deterministic, using up to 8 threads.

Root relaxation solution time = 0.06 sec. (7.94 ticks)

 

        Nodes                                         Cuts/

   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap         Variable B NodeID Parent  Depth

 

      0     0    38651.1003     2                  25483.7557      684         

      0     2    38651.1003     2                  38651.1003      684                                 0             0

Elapsed time = 5.09 sec. (93.75 ticks, tree = 0.01 MB, solutions = 0)

*     9     9      integral     0    42165.1593    38880.0334     2372    7.79%            Z(1) U     48     24      3

     10     6        cutoff          42165.1593    38880.0334     2867    7.79%            Z(1) D     41     25      4

     13     5    40567.1511     1    42165.1593    38880.0334     3446    7.79%            Z(6) U     49     16      2

*    15     2      integral     0    40529.4060    38880.0334     3801    4.07%            Z(1) D     35      3      4

 

User cuts applied:  156

 

Root node processing (before b&c):

  Real time             =    5.02 sec. (90.63 ticks)

Parallel b&c, 8 threads:

  Real time             =  547.36 sec. (2062.52 ticks)

  Sync time (average)   =  424.42 sec.

  Wait time (average)   =    0.00 sec.

                          ------------

Total (root+branch&cut) =  552.38 sec. (2153.15 ticks)

In addition, I've included time profiler report for both cases:

One Thread:

Total Time: 26.05

Time Spent outside Callback: 5.02

Time Spent in Lazy Constraint Callback (minus the time spent in separation method): 0.51

Time Spent in the Separation Method: 20.51

Time Spent Solving the Separation Problem: 15.90

Eight Threads:

Total Time: 549.87

Time Spent outside Callback: 7.86

Time Spent in Lazy Constraint Callback (minus the time spent in separation method): 102.85

Time Spent in the Separation Method: 439.16

Time Spent Solving the Separation Problem: 27.38

As one can see, in case of using eight threads, a lot of time is spent in the lazy constraint callback and the separation method. Perhaps, there is some kind of confusion on accessing the IloCplex object for the separation problem. In the code, I defined the separation problem problem IloCplex object and its related objects (i.e., variables, model, and constraints) and then only update its objective function each time the separation method is called (a similar approach to the ilobendersatsp example).

Originally posted by: EdKlotz

Ehsan,

If your models typically take only a handful of nodes, and the root node processing is uninteresting and consumes very little time, I doubt you will see any benefit using 8 threads.   The payoff comes either when there is lots of work at the root node, or the model requires a significant number of nodes to finish.   Based on you logs above as well as the model you sent me, that is not that case.    Also, note that most of the run time involved sychronization of the 8 threads:

What happens if you run with 8 threads but chance your lazy constraint callback so that you just enter the callback function and exit immediately without doing anything?

Regarding exporting SAV files, you are correct that since Lazy Constraints added via callback are considered part of the optimization strategy rather than the model, an exportModel() call will not include them.   But, LP format does support them at the start of the optimization; just add a 'Lazy Constraints' section to the LP file after the regular constraints but before the bounds section.   See

for more info.

Regarding the file being too large for a text editor, I was able to use vi on Unix to insert the lazy constraints.  Give that a try.   If that doesn't work, send me a private e-mail to klotz at us dot ibm dot com, and I can ftp the LP file to you (or have you download it from one of our ftp sites).

                                                                                                                                                                    Ed