00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "CouenneCutGenerator.hpp"
00013
00014 #include "BonCouenneInterface.hpp"
00015 #include "CouenneObject.hpp"
00016 #include "CouenneProblem.hpp"
00017 #include "CbcCutGenerator.hpp"
00018
00019 #include "BonAuxInfos.hpp"
00020 #include "CoinHelperFunctions.hpp"
00021 #include "BonOsiTMINLPInterface.hpp"
00022 #include "BonNlpHeuristic.hpp"
00023 #include "CouenneRecordBestSol.hpp"
00024
00025 using namespace Ipopt;
00026 using namespace Couenne;
00027
00028 NlpSolveHeuristic::NlpSolveHeuristic():
00029 CbcHeuristic(),
00030 nlp_(NULL),
00031 hasCloned_(false),
00032 maxNlpInf_(maxNlpInf_0),
00033 numberSolvePerLevel_(-1),
00034 couenne_(NULL){
00035 setHeuristicName("NlpSolveHeuristic");
00036 }
00037
00038 NlpSolveHeuristic::NlpSolveHeuristic(CbcModel & model, Bonmin::OsiTMINLPInterface &nlp, bool cloneNlp, CouenneProblem * couenne):
00039 CbcHeuristic(model), nlp_(&nlp), hasCloned_(cloneNlp),maxNlpInf_(maxNlpInf_0),
00040 numberSolvePerLevel_(-1),
00041 couenne_(couenne){
00042 setHeuristicName("NlpSolveHeuristic");
00043 if(cloneNlp)
00044 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (nlp.clone());
00045 }
00046
00047 NlpSolveHeuristic::NlpSolveHeuristic(const NlpSolveHeuristic & other):
00048 CbcHeuristic(other), nlp_(other.nlp_),
00049 hasCloned_(other.hasCloned_),
00050 maxNlpInf_(other.maxNlpInf_),
00051 numberSolvePerLevel_(other.numberSolvePerLevel_),
00052 couenne_(other.couenne_){
00053 if(hasCloned_ && nlp_ != NULL)
00054 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (other.nlp_->clone());
00055 }
00056
00057 CbcHeuristic *
00058 NlpSolveHeuristic::clone() const{
00059 return new NlpSolveHeuristic(*this);
00060 }
00061
00062 NlpSolveHeuristic &
00063 NlpSolveHeuristic::operator=(const NlpSolveHeuristic & rhs){
00064 if(this != &rhs){
00065 CbcHeuristic::operator=(rhs);
00066 if(hasCloned_ && nlp_)
00067 delete nlp_;
00068
00069 hasCloned_ = rhs.hasCloned_;
00070 if(nlp_ != NULL){
00071 if(hasCloned_)
00072 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (rhs.nlp_->clone());
00073 else
00074 nlp_ = rhs.nlp_;
00075 }
00076 }
00077 maxNlpInf_ = rhs.maxNlpInf_;
00078 numberSolvePerLevel_ = rhs.numberSolvePerLevel_;
00079 couenne_ = rhs.couenne_;
00080 return *this;
00081 }
00082
00083 NlpSolveHeuristic::~NlpSolveHeuristic(){
00084 if(hasCloned_)
00085 delete nlp_;
00086 nlp_ = NULL;
00087 }
00088
00089 void
00090 NlpSolveHeuristic::setNlp (Bonmin::OsiTMINLPInterface &nlp, bool cloneNlp){
00091 if(hasCloned_ && nlp_ != NULL)
00092 delete nlp_;
00093 hasCloned_ = cloneNlp;
00094 if(cloneNlp)
00095 nlp_ = dynamic_cast <Bonmin::OsiTMINLPInterface *> (nlp.clone());
00096 else
00097 nlp_ = &nlp;
00098 }
00099
00100 void
00101 NlpSolveHeuristic::setCouenneProblem(CouenneProblem * couenne)
00102 {couenne_ = couenne;}
00103
00104
00105 int
00106 NlpSolveHeuristic::solution (double & objectiveValue, double * newSolution) {
00107
00108 int noSolution = 1, maxTime = 2;
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120 const int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
00121
00122 if (depth <= 0)
00123 couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "NLP Heuristic: "); fflush (stdout);
00124
00125 try {
00126
00127 if (CoinCpuTime () > couenne_ -> getMaxCpuTime ())
00128 throw maxTime;
00129
00130 OsiSolverInterface * solver = model_ -> solver();
00131
00132 OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
00133 Bonmin::BabInfo * babInfo = dynamic_cast <Bonmin::BabInfo *> (auxInfo);
00134
00135 if(babInfo){
00136 babInfo->setHasNlpSolution(false);
00137 if(babInfo->infeasibleNode()){
00138 throw noSolution;
00139 }
00140 }
00141
00142
00143
00144 bool too_deep = false;
00145
00146
00147 if (numberSolvePerLevel_ > -1) {
00148
00149 if (numberSolvePerLevel_ == 0)
00150 throw maxTime;
00151
00152
00153 if (CoinDrand48 () > 1. / CoinMax
00154 (1., (double) ((depth - numberSolvePerLevel_) * (depth - numberSolvePerLevel_))))
00155 too_deep = true;
00156 }
00157
00158 if (too_deep)
00159 throw maxTime;
00160
00161 double *lower = new double [couenne_ -> nVars ()];
00162 double *upper = new double [couenne_ -> nVars ()];
00163
00164 CoinFillN (lower, couenne_ -> nVars (), -COUENNE_INFINITY);
00165 CoinFillN (upper, couenne_ -> nVars (), COUENNE_INFINITY);
00166
00167 CoinCopyN (solver->getColLower(), nlp_ -> getNumCols (), lower);
00168 CoinCopyN (solver->getColUpper(), nlp_ -> getNumCols (), upper);
00169
00170
00171
00172
00173
00174
00175 const double * solution = solver->getColSolution();
00176 OsiBranchingInformation info (solver, true);
00177 const int & numberObjects = model_->numberObjects();
00178 OsiObject ** objects = model_->objects();
00179 double maxInfeasibility = 0;
00180
00181 bool haveRoundedIntVars = false;
00182
00183 for (int i = 0 ; i < numberObjects ; i++) {
00184
00185 CouenneObject * couObj = dynamic_cast <CouenneObject *> (objects [i]);
00186
00187 if (couObj) {
00188 if (too_deep) {
00189 int dummy;
00190 double infeas;
00191 maxInfeasibility = CoinMax ( maxInfeasibility, infeas = couObj->infeasibility(&info, dummy));
00192
00193 if (maxInfeasibility > maxNlpInf_){
00194 delete [] lower;
00195 delete [] upper;
00196 throw noSolution;
00197 }
00198 }
00199 } else {
00200
00201 OsiSimpleInteger * intObj = dynamic_cast<OsiSimpleInteger *>(objects[i]);
00202
00203 if (intObj) {
00204 const int & i = intObj -> columnNumber ();
00205
00206 double value = solution [i];
00207 if (value < lower[i])
00208 value = lower[i];
00209 else if (value > upper[i])
00210 value = upper[i];
00211
00212 double rounded = floor (value + 0.5);
00213
00214 if (fabs (value - rounded) > COUENNE_EPS) {
00215 haveRoundedIntVars = true;
00216
00217 }
00218
00219
00220
00221
00222 }
00223 else{
00224
00225
00226
00227
00228 }
00229 }
00230 }
00231
00232
00233
00234
00235 bool skipOnInfeasibility = false;
00236
00237 double *Y = new double [couenne_ -> nVars ()];
00238 CoinFillN (Y, couenne_ -> nVars (), 0.);
00239 CoinCopyN (solution, nlp_ -> getNumCols (), Y);
00240
00241
00242
00243
00244
00245
00246
00247
00248 if (haveRoundedIntVars)
00249 skipOnInfeasibility = (couenne_ -> getIntegerCandidate (solution, Y, lower, upper) < 0);
00250
00251
00252
00253
00254
00255
00256
00257
00258 bool foundSolution = false;
00259
00260 if (haveRoundedIntVars && skipOnInfeasibility)
00261
00262
00263 for (int i = couenne_ -> nOrigVars (); i--;)
00264
00265 if (couenne_ -> Var (i) -> isDefinedInteger ())
00266 lower [i] = upper [i] = Y [i] =
00267 (CoinDrand48 () < 0.5) ?
00268 floor (Y [i] + COUENNE_EPS) :
00269 ceil (Y [i] - COUENNE_EPS);
00270
00271 else if (lower [i] > upper [i]) {
00272
00273
00274
00275
00276 double swap = lower [i];
00277 lower [i] = upper [i];
00278 upper [i] = swap;
00279 }
00280
00281 {
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 double * saveColLower = CoinCopyOfArray (nlp_ -> getColLower (), nlp_ -> getNumCols ());
00293 double * saveColUpper = CoinCopyOfArray (nlp_ -> getColUpper (), nlp_ -> getNumCols ());
00294
00295 for (int i = nlp_ -> getNumCols (); i--;) {
00296
00297 if (lower [i] > upper [i]) {
00298 double swap = lower [i];
00299 lower [i] = upper [i];
00300 upper [i] = swap;
00301 }
00302
00303 if (Y [i] < lower [i]) Y [i] = lower [i];
00304 else if (Y [i] > upper [i]) Y [i] = upper [i];
00305 }
00306
00307 nlp_ -> setColLower (lower);
00308 nlp_ -> setColUpper (upper);
00309 nlp_ -> setColSolution (Y);
00310
00311
00312 try {
00313 nlp_ -> options () -> SetNumericValue ("max_cpu_time", CoinMax (0., couenne_ -> getMaxCpuTime () - CoinCpuTime ()));
00314 nlp_ -> initialSolve ();
00315 }
00316 catch (Bonmin::TNLPSolver::UnsolvedError *E) {}
00317
00318 double obj = (nlp_ -> isProvenOptimal()) ? nlp_ -> getObjValue (): COIN_DBL_MAX;
00319
00320 if (nlp_ -> isProvenOptimal () &&
00321 couenne_ -> checkNLP (nlp_ -> getColSolution (), obj, true) &&
00322 (obj < couenne_ -> getCutOff ())) {
00323
00324
00325
00326 const int nVars = solver->getNumCols();
00327 double* tmpSolution = new double [nVars];
00328 CoinCopyN (nlp_ -> getColSolution(), nlp_ -> getNumCols(), tmpSolution);
00329
00330
00331 CouenneInterface * couenne = dynamic_cast <CouenneInterface *> (nlp_);
00332
00333 if (couenne)
00334 couenne_ -> getAuxs (tmpSolution);
00335
00336 #ifdef FM_CHECKNLP2
00337 if(!couenne_->checkNLP2(tmpSolution,
00338 0, false,
00339 true,
00340 false,
00341 couenne_->getFeasTol())) {
00342 #ifdef FM_USE_REL_VIOL_CONS
00343 printf("NlpSolveHeuristic::solution(): ### ERROR: checkNLP(): returns true, checkNLP2() returns false\n");
00344 exit(1);
00345 #endif
00346 }
00347 obj = couenne_->getRecordBestSol()->getModSolVal();
00348 couenne_->getRecordBestSol()->update();
00349 #else
00350 couenne_->getRecordBestSol()->update(tmpSolution, nVars,
00351 obj, couenne_->getFeasTol());
00352 #endif
00353
00354 if (babInfo){
00355 babInfo->setNlpSolution (tmpSolution, nVars, obj);
00356 babInfo->setHasNlpSolution (true);
00357 }
00358
00359 if (obj < objectiveValue) {
00360
00361 const CouNumber
00362 *lb = solver -> getColLower (),
00363 *ub = solver -> getColUpper ();
00364
00365
00366
00367 for (int i=0; i < nVars; i++, lb++, ub++) {
00368
00369 CouNumber &t = tmpSolution [i];
00370 if (t < *lb) t = *lb;
00371 else if (t > *ub) t = *ub;
00372 }
00373
00374
00375 couenne_ -> setCutOff (obj);
00376 foundSolution = true;
00377 objectiveValue = obj;
00378 CoinCopyN (tmpSolution, nVars, newSolution);
00379 }
00380 delete [] tmpSolution;
00381 }
00382
00383 nlp_->setColLower (saveColLower);
00384 nlp_->setColUpper (saveColUpper);
00385
00386 delete [] saveColLower;
00387 delete [] saveColUpper;
00388 }
00389
00390 delete [] Y;
00391
00392 delete [] lower;
00393 delete [] upper;
00394
00395 if (depth <= 0) {
00396
00397 if (foundSolution) couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
00398 else couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n");
00399 }
00400
00401 return foundSolution;
00402
00403 }
00404 catch (int &e) {
00405
00406
00407
00408 if ((couenne_ -> getCutOff () < objectiveValue) &&
00409 couenne_ -> getCutOffSol ()) {
00410
00411 objectiveValue = couenne_ -> getCutOff ();
00412 CoinCopyN (couenne_ -> getCutOffSol (), couenne_ -> nVars (), newSolution);
00413
00414 if (depth <= 0)
00415 couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "solution found, obj. %g\n", objectiveValue);
00416
00417 return 1;
00418
00419 } else {
00420
00421 if (depth <= 0 && e==noSolution)
00422 couenne_ -> Jnlst () -> Printf (J_ERROR, J_COUENNE, "no solution.\n", objectiveValue);
00423
00424 return 0;
00425 }
00426 }
00427 }
00428
00429
00431 void NlpSolveHeuristic::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
00432
00433 roptions -> AddStringOption2
00434 ("local_optimization_heuristic",
00435 "Search for local solutions of MINLPs",
00436 "yes",
00437 "no","",
00438 "yes","",
00439 "If enabled, a heuristic based on Ipopt is used to find feasible solutions for the problem. "
00440 "It is highly recommended that this option is left enabled, as it would be difficult to find feasible solutions otherwise.");
00441
00442 roptions -> AddLowerBoundedIntegerOption
00443 ("log_num_local_optimization_per_level",
00444 "Specify the logarithm of the number of local optimizations to perform"
00445 " on average for each level of given depth of the tree.",
00446 -1,
00447 2, "Solve as many nlp's at the nodes for each level of the tree. "
00448 "Nodes are randomly selected. If for a "
00449 "given level there are less nodes than this number nlp are solved for every nodes. "
00450 "For example if parameter is 8, nlp's are solved for all node until level 8, "
00451 "then for half the node at level 9, 1/4 at level 10.... "
00452 "Value -1 specify to perform at all nodes.");
00453 }