Two years ago, I
posted an example of how to implement
Benders decomposition in
CPLEX using the Java API. At the time, I believe the current version of CPLEX was 12.4; as of this writing, it is 12.6.0.1. Around version 12.5, IBM refactored the Java API for CPLEX and, in the process, made one or more non-backward-compatible changes that broke the sample code. This morning, I posted an updated version.
While I was at it, I decided to create a Git repository for the code on Bitbucket. You can access the repository at https://bitbucket.org/prubin/bendersexample. If you click the "Downloads" link in the navigation menu, you will be led to a page from which you can download the source code in a Zip archive. (There is no executable for this project, just the source.) There is also a link to an issue tracker in case you run into a bug (not that I
ever write code with bugs in it) and wish to report it.
UPDATE: I moved the code to a new repository. The URL is
https://gitlab.msu.edu/orobworld/BendersExample. Click "Files" in the left side navigation menu to get to the source code. There will be a link (upper right) on that page to download the source in a Zip archive.
The code and issue tracker should be accessible to anonymous users (no login required). If not, hopefully someone will let me know. As with the blog, the code is released under a rather permissive Creative Commons license.
Professor, I was checking the aggressive version. In the code where you generate and test the optimality cut, I think the portion of the code that specifies if a violation exists is missing. Or am I missing something?
ReplyDeleteProfessor, you say that the "aggressive" approach adds cuts at all nodes "by rounding fractional solutions". I don't quite see how you do this in the code. You take a solution at a node whether it is fractional or not, store these values as a HashMap, then create the cut using the rhs's of the constraints and the duals associated with them. Then you test this optimality cut to see if there is a violation. To my understanding we round up or down the fractional solution to the nearest integer and see if the cut is violated by the lhs and rhs obtained from these rounded solutions and you do this with
while (it.hasNext()) {
IloNumVar v = it.nextNumVar();
lhs += it.getValue()*msol.get(v);
}
// if a violation occurs, add an optimality cut
if (lhs > rhs + FUZZ) {
System.out.println("!!! Adding user optimality cut: " + r);
add(r);
}
I am also trying to understand what we exactly do inside this while loop. I would really appreciate if you could help me with these issues.
Best regards
Sorry, you did indeed find a bug in the code. I'm not sure how it got there, as I had previously tested the aggressive version. The bug caused the same cut to be generated and added repeatedly (although I think that the CPLEX cut pool manager would detect the duplication and not literally add it over and over). Any time a callback adds a cut, CPLEX immediately calls the callback again, to see if it wants to add more cuts. So the bug put CPLEX in a loop where it would poll the callback, get the same cut, add it (or not), then call the callback again, get the same cut, ...
Delete> // if a violation occurs, add an optimality cut
> if (lhs > rhs + FUZZ) {
> System.out.println("!!! Adding user optimality cut: " + r);
> add(r);
> }
This code looks correct (and is nearly identical to the patch I uploaded a few minutes ago), so I assume you added the if statement. Thanks for catching this!
> I am also trying to understand what we exactly do inside this while loop. I would really appreciate if you could help me with these issues.
Variable rexpr contains the linear expression in the cut we just generated. This cut is violated by the rounded values of the integer variables from the master problem. I want to plug the unrounded values in and see if it is still violated. If yes, we add the cut to the master; if no, we skip it.
There is no CPLEX method (or at least none I know) to evaluate a constraint on a vector of values, so we do this manually. Variable it iterates over the expression in the while loop. The method it.nextNumVar() advances the iterator to the next term in the expression and returns the variable in that term. The call to msol.get(v) retrieves the unrounded value of that variable in the current master solution (which I stored in a map in line 71). The call to it.getValue() returns the coefficient of that term in the linear expression. Inside the while loop, I multiply those and add the result to the running total in lhs. When the while loop terminates (after iterating over all terms in the expression, lhs contains the value of the left side of the new cut, and we compare that to the constant term (rhs). The cut says lhs <= rhs; if this is violated, we add the cut.
Did that clear things up?
Thank you Professor. It is crystal clear.
DeleteRespectfully,
Vedat Bayram
Dear Professor Robin,
DeleteFirst of all, I would appreciate your endeavors. I am actually a C++ user and am blind to Java. However, I have tried to understand the way the example has been coded. Unfortunately, I am still have problem with a small part of it. I do not understand what does the term "rhs" mean in the aggressive version. Would you please explain what does line 97 mean, i.e., (double rhs = r.getUB();)? What is the upper bound of r?
It has been mentioned that "lhs" contains the left side of the generated cut. Therefore, is "lhs" equal to "expr - flowCost"? or just "expr" excluding the constant term corresponding to the multiplication of duals of demand constraints (lambda) by RHS of demands constraints (Those ones which are )
Thank you in advance for your consideration.
r is a feasibility cut. Initially, it takes the form $a'x - b \le 0$, where $x$ is the vector of variables, $a$ is a coefficient vector and $b$ is a constant. CPLEX will automatically massage this into the form $-\infty \le a'x \le b$. The term "rhs" is the right-hand side ($b$) of this constraint, which CPLEX considers to be the upper bound of the range (hence the call to getUB() to recover it). If you asked CPLEX for the lower bound of the range, it would be $-\infty$.
DeleteThe left-hand side (lhs) would be everything other than the constant term.
May thanks for your prompt response. Clear as usual!
DeleteJust for the sake of clarification! In case of optimality cut, i,e, $a'x - b \le flowCost$, can one suppose "lhs" be equal to $a'x - flowCost$? I believe the answer to this question is "yes"; however, I just would like to make sure that my answer is correct.
Yes, that is correct.
DeleteProfessor,
ReplyDeleteI've been following your blog OR in OB World for a while and find it
very helpful. For that reason I would like to thank you for your kind
help.
I coded a two stage stochastic model using Java Cplex. In the first
stage facilities are located (binary) and in the second stage depending
on the scenario realised (four scenarios in my small example), facility
and route assignments are made. I do not need to check for the
feasibility because in the first stage (master problem) I ensure that
the subproblem is feasible by using induced constraints. So I just
create optimality cuts.
As I have read in many papers in the literaure along with your blog,
using lazy constraint callback has many advantages compared to the
classical approach where you solve the master problem to optimality at
every iteration. In the callback approach, which I use in my code, the
master problem is not solved to optimality and once an "integer
incumbent" solution is found it is passed to subproblem. So callback
class is called every time an integer incumbent solution found by the
master problem.
I have two files regarding this question, one is the output on my console when I run
the code and the other is the print out of the cplex node log. I actually sent an email to rubin@msu.edu attaching these two files, but then I thought maybe you are not using this email adress anymore. Anyway, I was checking these two files and contrary to what I stated in the above paragraph, I saw that the callback visits nodes that are not specified as
"integral" (or with a "*") by the node log as well. More interestingly
although those nodes are not marked as integral in node log, the binary
variables are all integral in the console printout. So there are two
possibilities, either callback interestingly visits some nodes although
they are not integral, or cplex node log somehow does not mark them as
integral. There is a third possibility actually that I am missing
something. I would really appreciate your comments.
Thank you for sharing all these ideas and thank you for your help in
advance.
Best regards
I did receive your email with the logs. I'm having a bit of trouble synchronizing the output with the node log, because the output contains a couple of "upper bound" entries early on that do not match any "best integer" values in the log. That aside, I think I can point out two things that may be contributing to your confusion.
DeleteFirst, the node log only prints an asterisk next to an _accepted_ integer-feasible solution. If CPLEX finds an integer-feasible solution at node N but a lazy constraint callback adds a cut that cuts off that solution, CPLEX resolves the node LP. When the cut is an optimality cut, it is possible that the same solution repeats with a more accurate objective value, and (after possibly more optimality cuts) is accepted and shows up in the node log with an asterisk. It is also possible, however, that the optimality cut shows the solution to actually be suboptimal, in which case no asterisk. Another possibility is that the solution is remains integer-feasible, the new objective value would be better than the current incumbent, but worse than some non-integer corner of the node LP. In that case, the new node LP solution will be non-integer, and I don't know whether CPLEX will spot the new incumbent.
Second, an integer-feasible node LP is only one possible trigger for a lazy constraint callback. Another possibility is that the current node's LP solution is non-integer, but CPLEX applies a heuristic that locates an integer-feasible solution. If the lazy constraint callback adds an optimality cut that causes that solution's objective value to become suboptimal, there will again be no asterisk.
I hope that helps.
Professor, I've been working on my code to see what could be the reason why the upper-bound that I compute does not match the best integer values in the log, as you mentioned in the first paragraph of your answer above.
ReplyDeleteI compute the upper-bound taking the expectation of the sub-problem objective values "sum(w, probOfScen[w]*sub[w].getObjectiveValue())". At the optimum solution also, the objective function value of the master problem "master.getObjectiveValue()" does not match the expectation of the objective values of sub-problems (scenarios). The difference is not too big, but not small either, it changes (increases as I noticed) as the number of scenarios increase. To check these results I also ran a separate model with extended two-stage stochastic formulation where I do not use Benders decomposition but extend the constraints and variables as many as the number of scenarios, as its name suggests. The results from this formulation matches with that of the expectation of sub-problems but not with that of the objective value of the master problem. Although I went over the code many times and tried so many different things, I can't find out what the reason for that could be. I would really appreciate your valuable comments.
Best regards,
Vedat Bayram
Is the master problem optimal value better or worse than the subproblem expectation? If the master objective value is better, is it better by more than the constraint tolerance (EpRhs)? If worse, is it worse by more than both the absolute and relative gaps (EpGap, EpAGap)?
DeleteIt is worse by more than absolute and relative gaps (my relative gap is 0 actually). And I think because of this problem, in some instances the optimum solution I get (opened shelter sites and naturally total evacuation time) is different from that of the extended formulation.
DeleteFirst, just to be clear, you are using a lazy constraint callback to add the Benders cuts? (Adding them in a user cut callback would be an error.)
DeleteMy guess is that some of your Benders cuts are incorrect. To test this, export the master problem to a file. Run the Benders code and either print or save to a file all the Benders cuts. Using a text editor, add those cuts to the master problem file, and also modify the lower and upper bounds of the integer variables to match the optimal solution of the extended formulation. Load that problem file into the CPLEX interactive optimizer and solve it. If CPLEX says it is infeasible, run the conflict refiner and see which Benders cuts are listed as contributing to the conflict, then take a close look at those.
This comment has been removed by the author.
ReplyDeleteProfessor, how fast is the aggressive benders compared to your version of benders with lazy constraint callback? And have you ever tried applying usercut callback and lazy constraint callback in the same class and adding the same type of optimality cut whether the incumbent solution is fractional or integer? Thank you.
ReplyDeleteBest regards,
Vedat Bayram
Vedat,
DeleteI have not actually tested the aggressive version on any real problems; I just coded it as a proof of concept. As far as using both callbacks, no, I have not done it on a real problem, but it would make some sense given the way CPLEX functions. The user cut callback will add the cut at nodes that do not have an integer LP solution. I'm not sure if it is called at nodes where the LP solution produces a new incumbent. Since user cuts are only supposed to tighten bounds, not cut off solutions, it would make sense to me to skip the user cut callback when the LP solution is integer-feasible. So the lazy constraint callback adds the cut in that case, and also in cases where the proposed incumbent is found by a heuristic rather than as an LP solution.
Dear Professor Rubin,
ReplyDeleteFirst of all, thanks for this very useful post. It helped me a lot implementing both the classical way for Benders decomposition and now also the new way using lazy constraint callbacks.
I ran into a problem. After the first cut is added using a lazy constraint callback, Cplex finds a new integral solution. However, instead of calling the lazy constraint callback, the node is cut off and Cplex quits with status "Unknown".
For debugging purposes, I printed the masterproblem to an .lp file. After comparing with the masterproblem after adding the first cut in the classical implementation, I conclude that is exactly the same problem. Also, when solving this printed problem using the interactive optimizer, Cplex finds an incumbent solution and has the status "Optimal". Therefore, I do not think I made a mistake in the way the cut is added or the cut itself.
I found this comment that might be related to the issue I am facing:
"A tricky situation occurs if the solution of one of the LP sub-problems in the branch-and-bound tree appears to be integer. In that case CPLEX will fathom the corresponding node (unless you create a child node using a branch callback) while you may want CPLEX to continue with the node."
(https://groups.google.com/forum/#!topic/aimms/7V0yppvL3Ng)
Do you think I ran into this tricky situation and do you happen to know a way of how to deal with this? Your help would be really appreciated.
Best regards,
Manon
Manon,
DeleteI don't think the comment you mention would apply to you, unless you are using a very old version of CPLEX. The comment was posted in 2011. That might be far enough back to be before CPLEX added lazy constraint callbacks. Prior to lazy constraint callbacks, you had to use an incumbent callback (which would let you accept or reject an incumbent but would not add a cut) combined with a cut callback (to add the new Benders cut globally) and a branch callback (to avoid the situation described in the comment, where the incumbent was an integer solution to the node LP, you rejected it, and CPLEX did not know what to do with the node). All that complexity was wiped away with the advent of the lazy constraint callback.
I suggest you look further at your lazy constraint callback. There may be a bug in it. Even if CPLEX incorrectly pruned the node without accepting the new incumbent, it would be odd to end with status "Unknown".
Good luck!
Dear professor Rubin,
ReplyDeleteThank you so much for your help and insight. You are right. I made a mistake in the lazy constraint callback and thanks to your comment I searched in the right direction. Instead of calling the add(IloRange r) function in the LazyConstraintCallback class, I called cplex.addLe(). It runs now and and the result is stunning. The computation time went down from 250 seconds using the classic implementation to only 6.6 seconds using callbacks! Thank you!!
I'm glad you were able to find the bug -- and thanks for posting the times. It's helpful to accumulate empirical evidence that the callback approach beats the classic implementation sometimes (usually?).
Delete