1 /* $Id$ */
2 // (C) Copyright International Business Machines Corporation 2007
3 // All Rights Reserved.
4 // This code is published under the Eclipse Public License (EPL).
5 //
6 // Authors :
7 // Pierre Bonami, International Business Machines Corporation
8 // Pietro Belotti, Lehigh University
9 //
10 // Date : 04/09/2007
11
12 #include "CouenneCutGenerator.hpp"
13
14 #include "BonCouenneInterface.hpp"
15 #include "CouenneObject.hpp"
16 #include "CouenneProblem.hpp"
17 #include "CbcCutGenerator.hpp"
18 //#include "CbcBranchActual.hpp"
19 #include "BonAuxInfos.hpp"
20 #include "CoinHelperFunctions.hpp"
21 #include "BonOsiTMINLPInterface.hpp"
22 #include "BonNlpHeuristic.hpp"
23 #include "CouenneRecordBestSol.hpp"
24
25 using namespace Ipopt;
26 using namespace Couenne;
27
NlpSolveHeuristic()28 NlpSolveHeuristic::NlpSolveHeuristic():
29 CbcHeuristic(),
30 nlp_(NULL),
31 hasCloned_(false),
32 maxNlpInf_(maxNlpInf_0),
33 numberSolvePerLevel_(-1),
34 couenne_(NULL){
35 setHeuristicName("NlpSolveHeuristic");
36 }
37
NlpSolveHeuristic(CbcModel & model,Bonmin::OsiTMINLPInterface & nlp,bool cloneNlp,CouenneProblem * couenne)38 NlpSolveHeuristic::NlpSolveHeuristic(CbcModel & model, Bonmin::OsiTMINLPInterface &nlp, bool cloneNlp, CouenneProblem * couenne):
39 CbcHeuristic(model), nlp_(&nlp), hasCloned_(cloneNlp),maxNlpInf_(maxNlpInf_0),
40 numberSolvePerLevel_(-1),
41 couenne_(couenne){
42 setHeuristicName("NlpSolveHeuristic");
43 if(cloneNlp)
44 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (nlp.clone());
45 }
46
NlpSolveHeuristic(const NlpSolveHeuristic & other)47 NlpSolveHeuristic::NlpSolveHeuristic(const NlpSolveHeuristic & other):
48 CbcHeuristic(other), nlp_(other.nlp_),
49 hasCloned_(other.hasCloned_),
50 maxNlpInf_(other.maxNlpInf_),
51 numberSolvePerLevel_(other.numberSolvePerLevel_),
52 couenne_(other.couenne_){
53 if(hasCloned_ && nlp_ != NULL)
54 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (other.nlp_->clone());
55 }
56
57 CbcHeuristic *
clone() const58 NlpSolveHeuristic::clone() const{
59 return new NlpSolveHeuristic(*this);
60 }
61
62 NlpSolveHeuristic &
operator =(const NlpSolveHeuristic & rhs)63 NlpSolveHeuristic::operator=(const NlpSolveHeuristic & rhs){
64 if(this != &rhs){
65 CbcHeuristic::operator=(rhs);
66 if(hasCloned_ && nlp_)
67 delete nlp_;
68
69 hasCloned_ = rhs.hasCloned_;
70 if(nlp_ != NULL){
71 if(hasCloned_)
72 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (rhs.nlp_->clone());
73 else
74 nlp_ = rhs.nlp_;
75 }
76 }
77 maxNlpInf_ = rhs.maxNlpInf_;
78 numberSolvePerLevel_ = rhs.numberSolvePerLevel_;
79 couenne_ = rhs.couenne_;
80 return *this;
81 }
82
~NlpSolveHeuristic()83 NlpSolveHeuristic::~NlpSolveHeuristic(){
84 if(hasCloned_)
85 delete nlp_;
86 nlp_ = NULL;
87 }
88
89 void
setNlp(Bonmin::OsiTMINLPInterface & nlp,bool cloneNlp)90 NlpSolveHeuristic::setNlp (Bonmin::OsiTMINLPInterface &nlp, bool cloneNlp){
91 if(hasCloned_ && nlp_ != NULL)
92 delete nlp_;
93 hasCloned_ = cloneNlp;
94 if(cloneNlp)
95 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (nlp.clone());
96 else
97 nlp_ = &nlp;
98 }
99
100 void
setCouenneProblem(CouenneProblem * couenne)101 NlpSolveHeuristic::setCouenneProblem(CouenneProblem * couenne)
102 {couenne_ = couenne;}
103
104
105 int
solution(double & objectiveValue,double * newSolution)106 NlpSolveHeuristic::solution (double & objectiveValue, double * newSolution) {
107
108 int noSolution = 1, maxTime = 2;
109
110 // do heuristic the usual way, but if for any reason (time is up, no
111 // better solution found) there is no improvement, get the best
112 // solution from the GlobalCutOff object in the pointer to the
113 // CouenneProblem and return it instead.
114 //
115 // Although this should be handled by Cbc, very often this doesn't
116 // happen.
117
118 // int nodeDepth = -1;
119
120 const int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
121
122 if (depth <= 0)
123 couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "NLP Heuristic: "); fflush (stdout);
124
125 try {
126
127 if (CoinCpuTime () > couenne_ -> getMaxCpuTime ())
128 throw maxTime;
129
130 OsiSolverInterface * solver = model_ -> solver();
131
132 OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
133 Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (auxInfo);
134
135 if(babInfo){
136 babInfo->setHasNlpSolution(false);
137 if(babInfo->infeasibleNode()){
138 throw noSolution;
139 }
140 }
141
142 // if too deep in the BB tree, only run NLP heuristic if
143 // feasibility is low
144 bool too_deep = false;
145
146 // check depth
147 if (numberSolvePerLevel_ > -1) {
148
149 if (numberSolvePerLevel_ == 0)
150 throw maxTime;
151
152 //if (CoinDrand48 () > pow (2., numberSolvePerLevel_ - depth))
153 if (CoinDrand48 () > 1. / CoinMax
154 (1., (double) ((depth - numberSolvePerLevel_) *
155 (depth - numberSolvePerLevel_))))
156 too_deep = true;
157 }
158
159 if (too_deep)
160 throw maxTime;
161
162 double *lower = new double [couenne_ -> nVars ()];
163 double *upper = new double [couenne_ -> nVars ()];
164
165 CoinFillN (lower, couenne_ -> nVars (), -COUENNE_INFINITY);
166 CoinFillN (upper, couenne_ -> nVars (), COUENNE_INFINITY);
167
168 CoinCopyN (solver->getColLower(), nlp_ -> getNumCols (), lower);
169 CoinCopyN (solver->getColUpper(), nlp_ -> getNumCols (), upper);
170
171 /*printf ("-- int candidate, before: ");
172 for (int i=0; i<couenne_ -> nOrig (); i++)
173 printf ("[%g %g] ", lower [i], upper [i]);
174 printf ("\n");*/
175
176 const double * solution = solver->getColSolution();
177 OsiBranchingInformation info (solver, true);
178 const int & numberObjects = model_->numberObjects();
179 OsiObject ** objects = model_->objects();
180 double maxInfeasibility = 0;
181
182 bool haveRoundedIntVars = false;
183
184 for (int i = 0 ; i < numberObjects ; i++) {
185
186 CouenneObject * couObj = dynamic_cast <CouenneObject *> (objects [i]);
187
188 if (couObj) {
189 if (too_deep) { // only test infeasibility if BB level is high
190 int dummy;
191 double infeas;
192 maxInfeasibility = CoinMax ( maxInfeasibility, infeas = couObj->infeasibility(&info, dummy));
193
194 if (maxInfeasibility > maxNlpInf_){
195 delete [] lower;
196 delete [] upper;
197 throw noSolution;
198 }
199 }
200 } else {
201
202 OsiSimpleInteger * intObj = dynamic_cast<OsiSimpleInteger *>(objects[i]);
203
204 if (intObj) {
205 const int & i = intObj -> columnNumber ();
206 // Round the variable in the solver
207 double value = solution [i];
208 if (value < lower[i])
209 value = lower[i];
210 else if (value > upper[i])
211 value = upper[i];
212
213 double rounded = floor (value + 0.5);
214
215 if (fabs (value - rounded) > COUENNE_EPS) {
216 haveRoundedIntVars = true;
217 //value = rounded;
218 }
219
220 // fix bounds anyway, if a better candidate is not found
221 // below at least we have an integer point
222 //lower[i] = upper[i] = value;
223 }
224 else{
225
226 // Probably a SOS object -- do not stop here
227 //throw CoinError("Bonmin::NlpSolveHeuristic","solution",
228 //"Unknown object.");
229 }
230 }
231 }
232
233 // if here, it means the infeasibility is not too high. Generate a
234 // better integer point as there are rounded integer variables
235
236 bool skipOnInfeasibility = false;
237
238 double *Y = new double [couenne_ -> nVars ()];
239 CoinFillN (Y, couenne_ -> nVars (), 0.);
240 CoinCopyN (solution, nlp_ -> getNumCols (), Y);
241
242 /*printf ("-- int candidate, upon call: ");
243 for (int i=0; i<couenne_ -> nOrig (); i++)
244 if (couenne_ -> Var (i) -> isInteger ())
245 printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
246 else printf ("%g ", Y [i]);
247 printf ("\n");*/
248
249 if (haveRoundedIntVars) // create "good" integer candidate for Ipopt
250 skipOnInfeasibility = (couenne_ -> getIntegerCandidate (solution, Y, lower, upper) < 0);
251
252 /*printf ("-- int candidate, after: ");
253 for (int i=0; i<couenne_ -> nOrig (); i++)
254 if (couenne_ -> Var (i) -> isInteger ())
255 printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
256 else printf ("%g ", Y [i]);
257 printf ("\n");*/
258
259 bool foundSolution = false;
260
261 if (haveRoundedIntVars && skipOnInfeasibility)
262 // no integer initial point could be found, make up some random rounding
263
264 for (int i = couenne_ -> nOrigVars (); i--;)
265
266 if (couenne_ -> Var (i) -> isDefinedInteger ())
267 lower [i] = upper [i] = Y [i] =
268 (CoinDrand48 () < 0.5) ?
269 floor (Y [i] + COUENNE_EPS) :
270 ceil (Y [i] - COUENNE_EPS);
271
272 else if (lower [i] > upper [i]) {
273
274 // sanity check (should avoid problems in ex1263 with
275 // couenne.opt.obbt)
276
277 double swap = lower [i];
278 lower [i] = upper [i];
279 upper [i] = swap;
280 }
281
282 {
283 // printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
284
285 /*printf ("int candidate: ");
286 for (int i=0; i<couenne_ -> nOrig (); i++)
287 if (couenne_ -> Var (i) -> isInteger ())
288 printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
289 else printf ("%g ", Y [i]);
290 printf ("\n");*/
291
292 // Now set column bounds and solve the NLP with starting point
293 double * saveColLower = CoinCopyOfArray (nlp_ -> getColLower (), nlp_ -> getNumCols ());
294 double * saveColUpper = CoinCopyOfArray (nlp_ -> getColUpper (), nlp_ -> getNumCols ());
295
296 for (int i = nlp_ -> getNumCols (); i--;) {
297
298 if (lower [i] > upper [i]) {
299 double swap = lower [i];
300 lower [i] = upper [i];
301 upper [i] = swap;
302 }
303
304 if (Y [i] < lower [i]) Y [i] = lower [i];
305 else if (Y [i] > upper [i]) Y [i] = upper [i];
306 }
307
308 nlp_ -> setColLower (lower);
309 nlp_ -> setColUpper (upper);
310 nlp_ -> setColSolution (Y);
311
312 // apply NLP solver /////////////////////////////////
313 try {
314 nlp_ -> options () -> SetNumericValue ("max_cpu_time", CoinMax (0.1, couenne_ -> getMaxCpuTime () - CoinCpuTime ()));
315 nlp_ -> initialSolve ();
316 }
317 catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
318
319 double obj = (nlp_ -> isProvenOptimal()) ? nlp_ -> getObjValue (): COIN_DBL_MAX;
320
321 if (nlp_ -> isProvenOptimal () &&
322 couenne_ -> checkNLP (nlp_ -> getColSolution (), obj, true) && // true for recomputing obj
323 (obj < couenne_ -> getCutOff ())) {
324
325 // store solution in Aux info
326
327 const int nVars = solver->getNumCols();
328 double* tmpSolution = new double [nVars];
329 CoinCopyN (nlp_ -> getColSolution(), nlp_ -> getNumCols(), tmpSolution);
330
331 //Get correct values for all auxiliary variables
332 CouenneInterface * couenne = dynamic_cast <CouenneInterface *> (nlp_);
333
334 if (couenne)
335 couenne_ -> getAuxs (tmpSolution);
336
337 #ifdef FM_CHECKNLP2
338 if(!couenne_->checkNLP2(tmpSolution,
339 0, false, // do not care about obj
340 true, // stopAtFirstViol
341 false, // checkAll
342 couenne_->getFeasTol())) {
343 #ifdef FM_USE_REL_VIOL_CONS
344 printf("NlpSolveHeuristic::solution(): ### ERROR: checkNLP(): returns true, checkNLP2() returns false\n");
345 exit(1);
346 #endif
347 }
348 obj = couenne_->getRecordBestSol()->getModSolVal();
349 couenne_->getRecordBestSol()->update();
350 #else
351 couenne_->getRecordBestSol()->update(tmpSolution, nVars,
352 obj, couenne_->getFeasTol());
353 #endif
354
355 if (babInfo){
356 babInfo->setNlpSolution (tmpSolution, nVars, obj);
357 babInfo->setHasNlpSolution (true);
358 }
359
360 if (obj < objectiveValue) { // found better solution?
361
362 const CouNumber
363 *lb = solver -> getColLower (),
364 *ub = solver -> getColUpper ();
365
366 // check bounds once more after getAux. This avoids false
367 // asserts in CbcModel.cpp:8305 on integerTolerance violated
368 for (int i=0; i < nVars; i++, lb++, ub++) {
369
370 CouNumber &t = tmpSolution [i];
371 if (t < *lb) t = *lb;
372 else if (t > *ub) t = *ub;
373 }
374
375 //printf ("new cutoff %g from BonNlpHeuristic\n", obj);
376 couenne_ -> setCutOff (obj);
377 foundSolution = true;
378 objectiveValue = obj;
379 CoinCopyN (tmpSolution, nVars, newSolution);
380 }
381 delete [] tmpSolution;
382 }
383
384 nlp_ -> setColLower (saveColLower);
385 nlp_ -> setColUpper (saveColUpper);
386
387 delete [] saveColLower;
388 delete [] saveColUpper;
389 }
390
391 delete [] Y;
392
393 delete [] lower;
394 delete [] upper;
395
396 if (depth <= 0) {
397
398 if (foundSolution) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
399 else couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n");
400 }
401
402 return foundSolution;
403
404 }
405 catch (int &e) {
406
407 // forget about using the global cutoff. That has to trickle up to
408 // Cbc some other way
409
410 if (e==noSolution) {if (depth <= 0) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n"); return 0;}
411 else if (e==maxTime) {if (depth <= 0) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "time limit reached.\n"); return 0;}
412 else {if (depth <= 0) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue); return 1;}
413
414 // // no solution available? Use the one from the global cutoff
415 // if ((couenne_ -> getCutOff () < objectiveValue) &&
416 // couenne_ -> getCutOffSol ()) {
417 // objectiveValue = couenne_ -> getCutOff ();
418 // CoinCopyN (couenne_ -> getCutOffSol (), couenne_ -> nVars (), newSolution);
419 // if (depth <= 0)
420 // couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
421 // return 1;
422 // } else {
423 // if (depth <= 0 && e==noSolution)
424 // couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n", objectiveValue);
425 // return 0;
426 // }
427 }
428 }
429
430
431 /// initialize options
registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)432 void NlpSolveHeuristic::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
433
434 roptions -> AddStringOption2
435 ("local_optimization_heuristic",
436 "Search for local solutions of MINLPs",
437 "yes",
438 "no","",
439 "yes","",
440 "If enabled, a heuristic based on Ipopt is used to find feasible solutions for the problem. "
441 "It is highly recommended that this option is left enabled, as it would be difficult to find feasible solutions otherwise.");
442
443 roptions -> AddLowerBoundedIntegerOption
444 ("log_num_local_optimization_per_level",
445 "Specify the logarithm of the number of local optimizations to perform"
446 " on average for each level of given depth of the tree.",
447 -1,
448 2, "Solve as many nlp's at the nodes for each level of the tree. "
449 "Nodes are randomly selected. If for a "
450 "given level there are less nodes than this number nlp are solved for every nodes. "
451 "For example if parameter is 8, nlp's are solved for all node until level 8, "
452 "then for half the node at level 9, 1/4 at level 10.... "
453 "Value -1 specify to perform at all nodes.");
454 }
455