source: stable/0.4/Couenne/src/heuristics/BonNlpHeuristic.cpp @ 905

Last change on this file since 905 was 905, checked in by pbelotti, 8 years ago

merged change 904 from trunk

File size: 13.1 KB
Line 
1/* $Id: BonNlpHeuristic.cpp 905 2012-09-04 04:13:26Z pbelotti $ */
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
25using namespace Ipopt;
26using namespace Couenne;
27
28NlpSolveHeuristic::NlpSolveHeuristic():
29  CbcHeuristic(),
30  nlp_(NULL),
31  hasCloned_(false),
32  maxNlpInf_(maxNlpInf_0),
33  numberSolvePerLevel_(-1),
34  couenne_(NULL){
35  setHeuristicName("NlpSolveHeuristic");
36}
37 
38NlpSolveHeuristic::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 
47NlpSolveHeuristic::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 
57CbcHeuristic * 
58NlpSolveHeuristic::clone() const{
59  return new NlpSolveHeuristic(*this);
60}
61 
62NlpSolveHeuristic &
63NlpSolveHeuristic::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 
83NlpSolveHeuristic::~NlpSolveHeuristic(){
84  if(hasCloned_)
85    delete nlp_;
86  nlp_ = NULL;
87}
88 
89void
90NlpSolveHeuristic::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 
100void
101NlpSolveHeuristic::setCouenneProblem(CouenneProblem * couenne)
102{couenne_ = couenne;}
103
104
105int
106NlpSolveHeuristic::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_) * (depth - numberSolvePerLevel_))))
155      too_deep = true;
156  }
157
158  if (too_deep)
159    throw maxTime;
160
161  double *lower = new double [couenne_ -> nVars ()];
162  double *upper = new double [couenne_ -> nVars ()];
163
164  CoinFillN (lower, couenne_ -> nVars (), -COUENNE_INFINITY);
165  CoinFillN (upper, couenne_ -> nVars (),  COUENNE_INFINITY);
166
167  CoinCopyN (solver->getColLower(), nlp_ -> getNumCols (), lower);
168  CoinCopyN (solver->getColUpper(), nlp_ -> getNumCols (), upper);
169
170  /*printf ("-- int candidate, before: ");
171    for (int i=0; i<couenne_ -> nOrig (); i++)
172    printf ("[%g %g] ", lower [i], upper [i]);
173    printf ("\n");*/
174
175  const double * solution = solver->getColSolution();
176  OsiBranchingInformation info (solver, true);
177  const int & numberObjects = model_->numberObjects();
178  OsiObject ** objects = model_->objects();
179  double maxInfeasibility = 0;
180
181  bool haveRoundedIntVars = false;
182
183  for (int i = 0 ; i < numberObjects ; i++) {
184
185    CouenneObject * couObj = dynamic_cast <CouenneObject *> (objects [i]);
186
187    if (couObj) {
188      if (too_deep) { // only test infeasibility if BB level is high
189        int dummy;
190        double infeas;
191        maxInfeasibility = CoinMax ( maxInfeasibility, infeas = couObj->infeasibility(&info, dummy));
192
193        if (maxInfeasibility > maxNlpInf_){
194          delete [] lower;
195          delete [] upper;
196          throw noSolution;
197        }
198      }
199    } else {
200
201      OsiSimpleInteger * intObj = dynamic_cast<OsiSimpleInteger *>(objects[i]);
202
203      if (intObj) {
204        const int & i = intObj -> columnNumber ();
205        // Round the variable in the solver
206        double value = solution [i];
207        if (value < lower[i])
208          value = lower[i];
209        else if (value > upper[i])
210          value = upper[i];
211
212        double rounded = floor (value + 0.5);
213
214        if (fabs (value - rounded) > COUENNE_EPS) {
215          haveRoundedIntVars = true;
216          //value = rounded;
217        }
218
219        // fix bounds anyway, if a better candidate is not found
220        // below at least we have an integer point
221        //lower[i] = upper[i] = value;
222      }
223      else{
224
225        // Probably a SOS object -- do not stop here
226        //throw CoinError("Bonmin::NlpSolveHeuristic","solution",
227        //"Unknown object.");
228      }
229    }
230  }
231
232  // if here, it means the infeasibility is not too high. Generate a
233  // better integer point as there are rounded integer variables
234
235  bool skipOnInfeasibility = false;
236
237  double *Y = new double [couenne_ -> nVars ()];
238  CoinFillN (Y, couenne_ -> nVars (), 0.);
239  CoinCopyN (solution, nlp_ -> getNumCols (), Y);
240
241  /*printf ("-- int candidate, upon call: ");
242    for (int i=0; i<couenne_ -> nOrig (); i++)
243    if (couenne_ -> Var (i) -> isInteger ())
244    printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
245    else printf ("%g ", Y [i]);
246    printf ("\n");*/
247
248  if (haveRoundedIntVars) // create "good" integer candidate for Ipopt
249    skipOnInfeasibility = (couenne_ -> getIntegerCandidate (solution, Y, lower, upper) < 0);
250
251  /*printf ("-- int candidate, after: ");
252    for (int i=0; i<couenne_ -> nOrig (); i++)
253    if (couenne_ -> Var (i) -> isInteger ())
254    printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
255    else printf ("%g ", Y [i]);
256    printf ("\n");*/
257
258  bool foundSolution = false;
259
260  if (haveRoundedIntVars && skipOnInfeasibility) 
261    // no integer initial point could be found, make up some random rounding
262
263    for (int i = couenne_ -> nOrigVars (); i--;) 
264
265      if (couenne_ -> Var (i) -> isDefinedInteger ())
266        lower [i] = upper [i] = Y [i] = 
267          (CoinDrand48 () < 0.5) ? 
268          floor (Y [i] + COUENNE_EPS) : 
269          ceil  (Y [i] - COUENNE_EPS);
270
271      else if (lower [i] > upper [i]) { 
272
273        // sanity check (should avoid problems in ex1263 with
274        // couenne.opt.obbt)
275
276        double swap = lower [i];
277        lower [i] = upper [i];
278        upper [i] = swap;
279      }
280
281  {
282    //  printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
283
284    /*printf ("int candidate: ");
285      for (int i=0; i<couenne_ -> nOrig (); i++)
286      if (couenne_ -> Var (i) -> isInteger ())
287      printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
288      else printf ("%g ", Y [i]);
289      printf ("\n");*/
290
291    // Now set column bounds and solve the NLP with starting point
292    double * saveColLower = CoinCopyOfArray (nlp_ -> getColLower (), nlp_ -> getNumCols ());
293    double * saveColUpper = CoinCopyOfArray (nlp_ -> getColUpper (), nlp_ -> getNumCols ());
294
295    for (int i = nlp_ -> getNumCols (); i--;) {
296
297      if (lower [i] > upper [i]) {
298        double swap = lower [i];
299        lower [i] = upper [i];
300        upper [i] = swap;
301      }
302
303      if      (Y [i] < lower [i]) Y [i] = lower [i];
304      else if (Y [i] > upper [i]) Y [i] = upper [i];
305    }
306
307    nlp_ -> setColLower    (lower);
308    nlp_ -> setColUpper    (upper);
309    nlp_ -> setColSolution (Y);
310
311    // apply NLP solver /////////////////////////////////
312    try {
313      nlp_ -> options () -> SetNumericValue ("max_cpu_time", CoinMax (0., couenne_ -> getMaxCpuTime () - CoinCpuTime ()));
314      nlp_ -> initialSolve ();
315    }
316    catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
317
318    double obj = (nlp_ -> isProvenOptimal()) ? nlp_ -> getObjValue (): COIN_DBL_MAX;
319
320    if (nlp_ -> isProvenOptimal () &&
321        couenne_ -> checkNLP (nlp_ -> getColSolution (), obj, true) && // true for recomputing obj
322        (obj < couenne_ -> getCutOff ())) {
323
324      // store solution in Aux info
325
326      const int nVars = solver->getNumCols();
327      double* tmpSolution = new double [nVars];
328      CoinCopyN (nlp_ -> getColSolution(), nlp_ -> getNumCols(), tmpSolution);
329
330      //Get correct values for all auxiliary variables
331      CouenneInterface * couenne = dynamic_cast <CouenneInterface *> (nlp_);
332
333      if (couenne)
334        couenne_ -> getAuxs (tmpSolution);
335
336#ifdef FM_CHECKNLP2
337      if(!couenne_->checkNLP2(tmpSolution, 
338                              0, false, // do not care about obj
339                              true, // stopAtFirstViol
340                              false, // checkAll
341                              couenne_->getFeasTol())) {
342#ifdef FM_USE_REL_VIOL_CONS
343        printf("NlpSolveHeuristic::solution(): ### ERROR: checkNLP(): returns true,  checkNLP2() returns false\n");
344        exit(1);
345#endif
346      }
347      obj = couenne_->getRecordBestSol()->getModSolVal(); 
348      couenne_->getRecordBestSol()->update();
349#else
350      couenne_->getRecordBestSol()->update(tmpSolution, nVars,
351                                           obj, couenne_->getFeasTol());
352#endif
353
354      if (babInfo){
355        babInfo->setNlpSolution (tmpSolution, nVars, obj);
356        babInfo->setHasNlpSolution (true);
357      }
358
359      if (obj < objectiveValue) { // found better solution?
360
361        const CouNumber
362          *lb = solver -> getColLower (),
363          *ub = solver -> getColUpper ();
364
365        // check bounds once more after getAux. This avoids false
366        // asserts in CbcModel.cpp:8305 on integerTolerance violated
367        for (int i=0; i < nVars; i++, lb++, ub++) {
368
369          CouNumber &t = tmpSolution [i];
370          if      (t < *lb) t = *lb;
371          else if (t > *ub) t = *ub;
372        }
373         
374        //printf ("new cutoff %g from BonNlpHeuristic\n", obj);
375        couenne_ -> setCutOff (obj);
376        foundSolution = true;
377        objectiveValue = obj;
378        CoinCopyN (tmpSolution, nVars, newSolution);
379      }
380      delete [] tmpSolution;
381    }
382
383    nlp_->setColLower (saveColLower);
384    nlp_->setColUpper (saveColUpper);
385
386    delete [] saveColLower;
387    delete [] saveColUpper;
388  }
389
390  delete [] Y;
391
392  delete [] lower;
393  delete [] upper;
394
395  if (depth <= 0) {
396
397    if (foundSolution) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
398    else               couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n");
399  }
400
401  return foundSolution;
402
403  }
404  catch (int &e) {
405
406    return 0;
407
408    // no solution available? Use the one from the global cutoff
409
410    if ((couenne_ -> getCutOff () < objectiveValue) &&
411        couenne_ -> getCutOffSol ()) {
412
413      objectiveValue = couenne_ -> getCutOff    ();
414      CoinCopyN       (couenne_ -> getCutOffSol (), couenne_ -> nVars (), newSolution);
415
416      if (depth <= 0)
417        couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
418
419      return 1;
420
421    } else {
422
423      if (depth <= 0 && e==noSolution)
424        couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n", objectiveValue);
425
426      return 0;
427    }
428  }
429}
430
431
432/// initialize options
433void NlpSolveHeuristic::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
434
435  roptions -> AddStringOption2
436    ("local_optimization_heuristic",
437     "Search for local solutions of MINLPs",
438     "yes",
439     "no","",
440     "yes","",
441     "If enabled, a heuristic based on Ipopt is used to find feasible solutions for the problem. "
442     "It is highly recommended that this option is left enabled, as it would be difficult to find feasible solutions otherwise.");
443
444  roptions -> AddLowerBoundedIntegerOption
445    ("log_num_local_optimization_per_level",
446     "Specify the logarithm of the number of local optimizations to perform" 
447     " on average for each level of given depth of the tree.",
448     -1,
449     2, "Solve as many nlp's at the nodes for each level of the tree. "
450     "Nodes are randomly selected. If for a "
451     "given level there are less nodes than this number nlp are solved for every nodes. "
452     "For example if parameter is 8, nlp's are solved for all node until level 8, " 
453     "then for half the node at level 9, 1/4 at level 10.... "
454     "Value -1 specify to perform at all nodes.");
455}
Note: See TracBrowser for help on using the repository browser.