00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouenneIterativeRounding.hpp"
00012 #include "BonTMINLP2Quad.hpp"
00013 #include "BonTMINLPLinObj.hpp"
00014 #ifdef COIN_HAS_CPX
00015 #include "OsiCpxSolverInterface.hpp"
00016 #include "cplex.h"
00017 #endif
00018
00019 #include "CouenneRecordBestSol.hpp"
00020
00021 #define MILPTIME 5
00022 #define CBCMILPTIME 20
00023
00024 namespace Couenne{
00025
00026 CouenneIterativeRounding::CouenneIterativeRounding():
00027 CbcHeuristic(),
00028 nlp_(NULL), cinlp_(NULL), milp_(NULL),
00029 maxRoundingIter_(10),
00030 maxFirPoints_(5), maxTime_(60), maxTimeFirstCall_(60), numInitialRows_(0),
00031 numSol_(-1), colLower_(NULL), colUpper_(NULL),
00032 colLowerNlp_(NULL), colUpperNlp_(NULL),
00033 omega_(0.2), baseLbRhs_(15),
00034 couenne_(NULL)
00035 {
00036 setHeuristicName("CouenneIterativeRounding");
00037 }
00038
00039 CouenneIterativeRounding::CouenneIterativeRounding(Bonmin::OsiTMINLPInterface* nlp,
00040 OsiSolverInterface* cinlp,
00041 CouenneProblem* couenne,
00042 Ipopt::SmartPtr<Ipopt::OptionsList> options):
00043 CbcHeuristic(), nlp_(NULL),
00044 cinlp_(NULL), milp_(NULL),
00045 numSol_(-1), colLower_(NULL), colUpper_(NULL),
00046 colLowerNlp_(NULL), colUpperNlp_(NULL),
00047 couenne_(couenne)
00048 {
00049
00050 setNlp(nlp, cinlp);
00051
00052
00053 int irAggressiveness;
00054 options->GetIntegerValue("iterative_rounding_aggressiveness",
00055 irAggressiveness, "couenne.");
00056 setAggressiveness(irAggressiveness);
00057 double maxTime, maxTimeInit;
00058 options->GetNumericValue("iterative_rounding_time", maxTime,
00059 "couenne.");
00060 options->GetNumericValue("iterative_rounding_time_firstcall",
00061 maxTimeInit, "couenne.");
00062 if (maxTime >= 0){
00063 setMaxTime(maxTime);
00064 }
00065 if (maxTimeInit >= 0){
00066 setMaxTimeFirstCall(maxTimeInit);
00067 }
00068 int irLbrhs;
00069 options->GetIntegerValue("iterative_rounding_base_lbrhs",
00070 irLbrhs, "couenne.");
00071 setBaseLbRhs(irLbrhs);
00072 int numFirPoints;
00073 options->GetIntegerValue("iterative_rounding_num_fir_points",
00074 numFirPoints, "couenne.");
00075 setMaxFirPoints(numFirPoints);
00076 double omega;
00077 options->GetNumericValue("iterative_rounding_omega",
00078 omega, "couenne.");
00079 setOmega(omega);
00080 }
00081
00082 CouenneIterativeRounding::CouenneIterativeRounding(const CouenneIterativeRounding & other):
00083 CbcHeuristic(other), nlp_(other.nlp_),
00084 cinlp_(other.cinlp_), milp_(other.milp_),
00085 maxRoundingIter_(other.maxRoundingIter_),
00086 maxFirPoints_(other.maxFirPoints_),
00087 maxTime_(other.maxTime_),
00088 maxTimeFirstCall_(other.maxTimeFirstCall_),
00089 numInitialRows_(other.numInitialRows_),
00090 numSol_(other.numSol_),
00091 omega_(other.omega_), baseLbRhs_(other.baseLbRhs_),
00092 couenne_(other.couenne_)
00093 {
00094 if(nlp_ != NULL){
00095 nlp_ = dynamic_cast<Bonmin::OsiTMINLPInterface*>(other.nlp_->clone());
00096 }
00097 if(milp_ != NULL)
00098 #ifdef COIN_HAS_CPX
00099 milp_ = dynamic_cast<OsiCpxSolverInterface*>(other.milp_->clone());
00100 #else
00101 milp_ = dynamic_cast<OsiClpSolverInterface*>(other.milp_->clone());
00102 #endif
00103 if (other.colLower_ != NULL){
00104 if (colLower_ != NULL)
00105 delete colLower_;
00106 colLower_ = new double[milp_->getNumCols()];
00107 CoinCopyN (other.colLower_, milp_->getNumCols(), colLower_);
00108 }
00109 if (other.colUpper_ != NULL){
00110 if (colUpper_ != NULL)
00111 delete colUpper_;
00112 colUpper_ = new double[milp_->getNumCols()];
00113 CoinCopyN (other.colUpper_, milp_->getNumCols(), colUpper_);
00114 }
00115 if (other.colLowerNlp_ != NULL){
00116 if (colLowerNlp_ != NULL)
00117 delete colLowerNlp_;
00118 colLowerNlp_ = new double[nlp_->getNumCols()];
00119 CoinCopyN (other.colLowerNlp_, nlp_->getNumCols(), colLowerNlp_);
00120 }
00121 if (other.colUpperNlp_ != NULL){
00122 if (colUpperNlp_ != NULL)
00123 delete colUpperNlp_;
00124 colUpperNlp_ = new double[nlp_->getNumCols()];
00125 CoinCopyN (other.colUpperNlp_, nlp_->getNumCols(), colLowerNlp_);
00126 }
00127 }
00128
00129 CbcHeuristic *
00130 CouenneIterativeRounding::clone() const{
00131 return new CouenneIterativeRounding(*this);
00132 }
00133
00134 CouenneIterativeRounding &
00135 CouenneIterativeRounding::operator=(const CouenneIterativeRounding & rhs){
00136 if(this != &rhs){
00137 CbcHeuristic::operator=(rhs);
00138 if(nlp_)
00139 delete nlp_;
00140
00141 if(rhs.nlp_ != NULL){
00142 nlp_ = dynamic_cast<Bonmin::OsiTMINLPInterface*>(rhs.nlp_->clone());
00143 }
00144 cinlp_ = rhs.cinlp_;
00145 maxRoundingIter_ = rhs.maxRoundingIter_;
00146 maxFirPoints_ = rhs.maxFirPoints_;
00147 maxTime_ = rhs.maxTime_;
00148 maxTimeFirstCall_ = rhs.maxTimeFirstCall_;
00149 numSol_ = rhs.numSol_;
00150 numInitialRows_ = rhs.numInitialRows_;
00151 omega_ = rhs.omega_;
00152 baseLbRhs_ = rhs.baseLbRhs_;
00153 couenne_ = rhs.couenne_;
00154 if (rhs.colLower_ != NULL){
00155 if (colLower_ != NULL)
00156 delete colLower_;
00157 colLower_ = new double[milp_->getNumCols()];
00158 CoinCopyN (rhs.colLower_, milp_->getNumCols(), colLower_);
00159 }
00160 if (rhs.colUpper_ != NULL){
00161 if (colUpper_ != NULL)
00162 delete colUpper_;
00163 colUpper_ = new double[milp_->getNumCols()];
00164 CoinCopyN (rhs.colUpper_, milp_->getNumCols(), colLower_);
00165 }
00166 if (rhs.colLowerNlp_ != NULL){
00167 if (colLowerNlp_ != NULL)
00168 delete colLowerNlp_;
00169 colLowerNlp_ = new double[nlp_->getNumCols()];
00170 CoinCopyN (rhs.colLowerNlp_, nlp_->getNumCols(), colLowerNlp_);
00171 }
00172 if (rhs.colUpperNlp_ != NULL){
00173 if (colUpperNlp_ != NULL)
00174 delete colUpperNlp_;
00175 colUpperNlp_ = new double[nlp_->getNumCols()];
00176 CoinCopyN (rhs.colUpperNlp_, nlp_->getNumCols(), colLowerNlp_);
00177 }
00178 }
00179 return *this;
00180 }
00181
00182 CouenneIterativeRounding::~CouenneIterativeRounding(){
00183 delete nlp_;
00184 nlp_ = NULL;
00185 if (colLower_)
00186 delete[] colLower_;
00187 if (colUpper_)
00188 delete[] colUpper_;
00189 if (colLowerNlp_)
00190 delete[] colLowerNlp_;
00191 if (colUpperNlp_)
00192 delete[] colUpperNlp_;
00193 if (milp_)
00194 delete milp_;
00195 milp_ = NULL;
00196 }
00197
00198 void
00199 CouenneIterativeRounding::setNlp(Bonmin::OsiTMINLPInterface* nlp,
00200 OsiSolverInterface * cinlp){
00201
00202
00203 if(nlp_ != NULL)
00204 delete nlp_;
00205 nlp_ = new Bonmin::OsiTMINLPInterface;
00206 Ipopt::SmartPtr<Bonmin::TMINLP> tminlp = nlp->model();
00207 if (tminlp->hasLinearObjective()){
00208 Ipopt::SmartPtr<Bonmin::TMINLPLinObj> linObj =
00209 new Bonmin::TMINLPLinObj;
00210 linObj->setTminlp(GetRawPtr(tminlp));
00211 tminlp = GetRawPtr(linObj);
00212 }
00213 nlp_->initialize(nlp->regOptions(), nlp->options(), nlp->solver()->journalist(), "bonmin.", tminlp);
00214 nlp_->use(new Bonmin::TMINLP2TNLPQuadCuts(tminlp));
00215 cinlp_ = cinlp;
00216 }
00217
00218 void
00219 CouenneIterativeRounding::setMilp(){
00220 if(milp_ != NULL)
00221 delete milp_;
00222
00223
00224
00225 OsiSolverInterface * milp = model_->solver();
00226 int n = milp->getNumCols();
00227
00228 #ifdef COIN_HAS_CPX
00229 milp_ = new OsiCpxSolverInterface();
00230 milp_->loadProblem(*(milp->getMatrixByRow()), milp->getColLower(),
00231 milp->getColUpper(), milp->getObjCoefficients(),
00232 milp->getRowLower(), milp->getRowUpper());
00233 for (int i = 0; i < n; ++i){
00234 if (milp->isInteger(i))
00235 milp_->setInteger(i);
00236 }
00237 #else
00238 milp_ = dynamic_cast<OsiClpSolverInterface*>(milp->clone());
00239 #endif
00240
00241 colLower_ = new double[n];
00242 colUpper_ = new double[n];
00243 memcpy(colLower_, milp->getColLower(), n*sizeof(double));
00244 memcpy(colUpper_, milp->getColUpper(), n*sizeof(double));
00245
00246 int nNlp = cinlp_->getNumCols();
00247
00248 colLowerNlp_ = new double[nNlp];
00249 colUpperNlp_ = new double[nNlp];
00250 memcpy(colLowerNlp_, cinlp_->getColLower(), nNlp*sizeof(double));
00251 memcpy(colUpperNlp_, cinlp_->getColUpper(), nNlp*sizeof(double));
00252
00253 numIntegers_ = 0;
00254 for (int i = 0; i < nNlp; ++i){
00255 if (cinlp_->isInteger(i)){
00256 numIntegers_++;
00257 }
00258 }
00259
00260
00261
00262
00263 double swap;
00264 for (int i = 0; i < n; ++i){
00265 if (colUpper_[i] < colLower_[i]){
00266 swap = colUpper_[i];
00267 colUpper_[i] = colLower_[i];
00268 colLower_[i] = swap;
00269 }
00270 }
00271
00272 numInitialRows_ = milp_->getNumRows();
00273
00274
00275 double * tmpArray = new double[n];
00276 CoinFillN(tmpArray, n, 0.0);
00277 milp_->setObjective(tmpArray);
00278 milp_->setObjSense(1);
00279
00280
00281
00282 for (int i = 0; i < n; ++i){
00283 milp_->addCol(0, NULL, NULL, 0.0, COIN_DBL_MAX, 1.0);
00284 }
00285
00286 milp_->setHintParam(OsiDoDualInResolve,true,OsiHintDo);
00287 milp_->setHintParam(OsiDoPresolveInResolve,true,OsiHintDo);
00288 milp_->setHintParam(OsiDoReducePrint,true,OsiHintDo);
00289 milp_->setDblParam(OsiPrimalTolerance, COUENNE_EPS_INT);
00290 milp_->messageHandler()->setLogLevel(0);
00291 milp_->setDblParam(OsiDualTolerance, COUENNE_EPS_INT);
00292
00293 #ifdef COIN_HAS_CPX
00294
00295 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_MIPEMPHASIS, CPX_MIPEMPHASIS_HIDDENFEAS);
00296 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_FPHEUR, 2);
00297 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_NODESEL, 0);
00298 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_LBHEUR, 1);
00299 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_RINSHEUR, 99);
00300 CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_NODELIM, 50);
00301 CPXsetdblparam(milp_->getEnvironmentPtr(), CPX_PARAM_CUTSFACTOR, 1.0);
00302 CPXsetdblparam(milp_->getEnvironmentPtr(), CPX_PARAM_TILIM, MILPTIME);
00303 #else
00304
00305 heuristics_ = new CbcHeuristic*[1];
00306 numHeuristics_ = 1;
00307 CbcHeuristicFPump * feaspump = new CbcHeuristicFPump();
00308 feaspump->setMaximumRetries(2);
00309 feaspump->setMaximumPasses(100);
00310 feaspump->setAccumulate(3);
00311 heuristics_[0] = feaspump;
00312 #endif
00313
00314 delete[] tmpArray;
00315
00316
00317 }
00318
00319 void CouenneIterativeRounding::setAggressiveness(int value){
00320 switch (value){
00321 case 0:
00322 setMaxRoundingIter(5);
00323 setMaxTimeFirstCall(300);
00324 setMaxFirPoints(5);
00325 setMaxTime(60);
00326 break;
00327 case 1:
00328 setMaxRoundingIter(10);
00329 setMaxTimeFirstCall(300);
00330 setMaxFirPoints(5);
00331 setMaxTime(120);
00332 break;
00333 case 2:
00334 setMaxRoundingIter(20);
00335 setMaxTimeFirstCall(1000);
00336 setMaxFirPoints(5);
00337 setMaxTime(300);
00338 break;
00339 default:
00340 std::cerr << "CouenneIterativeRounding::setAggressiveness() : unknown value!\n" << std::endl;
00341 }
00342 }
00343
00344 int
00345 CouenneIterativeRounding::solution(double & objectiveValue, double* newSolution){
00346 if (milp_ == NULL){
00347
00348 setMilp();
00349 return 0;
00350 }
00351
00352 if ((model_->numberIntegers() == 0) ||
00353 (numSol_ == model_->getSolutionCount())){
00354
00355
00356 return 0;
00357 }
00358
00359 numSol_ = model_->getSolutionCount();
00360
00361 std::cout << "Launching IterativeRounding with parameters:" << std::endl;
00362 std::cout << "Max rounding iter: " << maxRoundingIter_ << std::endl;
00363 std::cout << "Max feas point: " << maxFirPoints_ << std::endl;
00364 std::cout << "Base lbrhs: " << baseLbRhs_ << std::endl;
00365 std::cout << "Omega: " << omega_ << std::endl;
00366 std::cout << "Max time firstcall: " << maxTimeFirstCall_ << std::endl;
00367
00368
00369 startTime_ = CoinCpuTime();
00370 endTime_ = ((numSol_ == 0) ? maxTimeFirstCall_ : maxTime_);
00371
00372 const double* bestKnownSolution = model_->bestSolution();
00373 bool found = false;
00374 bool hasSolution = true;
00375 bool improved = true;
00376 if (numSol_ == 0){
00377
00378 hasSolution = feasibilityIR(objectiveValue, newSolution);
00379 if (hasSolution){
00380 bestKnownSolution = newSolution;
00381 found = hasSolution;
00382 }
00383 }
00384 if (!hasSolution){
00385
00386 return found;
00387 }
00388 while (improved && (CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)){
00389
00390 improved = false;
00391 improved = improvementIR(objectiveValue, newSolution, bestKnownSolution);
00392 if (improved){
00393 bestKnownSolution = newSolution;
00394 }
00395 found = (found || improved);
00396 }
00397 if (found){
00398
00399 numSol_++;
00400 }
00401
00402 return found;
00403 }
00404
00405 int
00406 CouenneIterativeRounding::feasibilityIR(double& objectiveValue,
00407 double* newSolution){
00408
00409 std::cout << "starting feasibility IR" << std::endl;
00410
00411 OsiSolverInterface * solver = model_->solver();
00412
00413 OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
00414 Bonmin::BabInfo * babInfo = dynamic_cast<Bonmin::BabInfo *> (auxInfo);
00415
00416 if(babInfo){
00417 babInfo->setHasNlpSolution(false);
00418 if(babInfo->infeasibleNode()){
00419 return 0;
00420 }
00421 }
00422
00423 int n = couenne_->nVars();
00424 int nNlp = cinlp_->getNumCols();
00425
00426
00427 int numIntAtBound = 0;
00428
00429 OsiRowCut cut;
00430 OsiRowCut lbcut1;
00431 OsiRowCut lbcut2;
00432 OsiCuts lbcuts;
00433 double obj;
00434
00435 bool boundsChanged = false;
00436 std::vector<int> previousBranches;
00437
00438
00439
00440
00441 nlp_->resolve();
00442 if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() ||
00443 nlp_->isAbandoned()){
00444 nlp_->resolveForRobustness(3);
00445 }
00446
00447 if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() ||
00448 nlp_->isAbandoned()){
00449 obj = COIN_DBL_MAX;
00450 }
00451 else{
00452 obj = nlp_->getObjValue();
00453 }
00454
00455 if (obj == COIN_DBL_MAX){std::cout << " IR: no feasible solution found " << std::endl;
00456 std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl;
00457 return 0;
00458 }
00459
00460 double* xprime = new double [n];
00461 CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime);
00462 couenne_ -> getAuxs (xprime);
00463
00464
00465
00466
00467 int constrInd[2];
00468 double constrElem[2];
00469
00470
00471 bool foundSolution = false;
00472 double* tmpSolution = new double[n];
00473 OsiCuts revlb_all;
00474 OsiSolverInterface* curr_milp = milp_;
00475
00476 int outerLoop = maxFirPoints_;
00477 for (int h = 0; h < outerLoop &&
00478 ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++h){
00479
00480 OsiCuts cs;
00481 cs.insert(lbcut1);
00482 cs.insert(lbcut2);
00483
00484
00485
00486 for (int i = 0; i < n; i++){
00487 constrInd[0] = i;
00488 constrInd[1] = i+n;
00489 constrElem[0] = -1;
00490 constrElem[1] = -1;
00491 cut.mutableRow().setVector(2, constrInd, constrElem);
00492 cut.setLb(-COIN_DBL_MAX);
00493 cut.setUb(-xprime[i]);
00494 cs.insert(cut);
00495 constrElem[0] = +1;
00496 constrElem[1] = -1;
00497 cut.mutableRow().setVector(2, constrInd, constrElem);
00498 cut.setLb(-COIN_DBL_MAX);
00499 cut.setUb(xprime[i]);
00500 cs.insert(cut);
00501 }
00502 curr_milp->applyCuts(cs);
00503
00504 for (int k = 0; k < maxRoundingIter_ &&
00505 ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++k){
00506
00507 curr_milp->restoreBaseModel(numInitialRows_+cs.sizeRowCuts());
00508 curr_milp->applyCuts(revlb_all);
00509 bool solFound = solveMilp(curr_milp,
00510 endTime_-(CoinCpuTime()-startTime_)-2);
00511 if (!solFound && !boundsChanged){
00512 std::cout << " MILP cannot be solved, terminating LB " << std::endl;
00513
00514 curr_milp->restoreBaseModel(numInitialRows_);
00515 delete[] xprime;
00516 delete[] tmpSolution;
00517 if (boundsChanged){
00518 curr_milp->setColLower(colLower_);
00519 curr_milp->setColUpper(colUpper_);
00520 }
00521 std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl;
00522 return foundSolution;
00523 }
00524 else if (!solFound && boundsChanged){
00525
00526
00527
00528 curr_milp->setColLower(colLower_);
00529 curr_milp->setColUpper(colUpper_);
00530 branchToCut(tmpSolution, curr_milp, previousBranches);
00531 continue;
00532 }
00533
00534 const double * xtilde = curr_milp->getColSolution();
00535
00536
00537
00538 CoinCopyN (xtilde, n, tmpSolution);
00539 for (int i = 0; i < nNlp; ++i){
00540 if (model_->isInteger(i)){
00541 tmpSolution[i] = floor(tmpSolution[i]+0.5);
00542 cinlp_->setColLower(i, tmpSolution[i]);
00543 cinlp_->setColUpper(i, tmpSolution[i]);
00544 }
00545 }
00546 cinlp_->setColSolution(tmpSolution);
00547
00548 cinlp_->messageHandler()->setLogLevel(1);
00549 cinlp_->resolve();
00550 obj = ((cinlp_->isProvenOptimal()) ?
00551 cinlp_->getObjValue():COIN_DBL_MAX);
00552 memcpy(tmpSolution, cinlp_->getColSolution(),
00553 nNlp*sizeof(double));
00554
00555
00556
00557 bool isChecked = false;
00558 #ifdef FM_CHECKNLP2
00559 isChecked = couenne_->checkNLP2(tmpSolution, 0, false,
00560 true,
00561 false,
00562 couenne_->getFeasTol());
00563 if(isChecked) {
00564 obj = couenne_->getRecordBestSol()->getModSolVal();
00565 }
00566 #else
00567 isChecked = couenne_->checkNLP(tmpSolution, obj, true);
00568 #endif
00569
00570 if (cinlp_->isProvenOptimal () &&
00571 isChecked &&
00572 (obj < couenne_->getCutOff())) {
00573
00574 #ifdef FM_CHECKNLP2
00575 #ifdef FM_TRACE_OPTSOL
00576 couenne_->getRecordBestSol()->update();
00577 CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution);
00578 obj = couenne_->getRecordBestSol()->getVal();
00579 #else
00580 CoinCopyN (couenne_->getRecordBestSol()->getModSol(n), n, tmpSolution);
00581 #endif
00582 #else
00583
00584
00585 couenne_ -> getAuxs (tmpSolution);
00586
00587 #ifdef FM_TRACE_OPTSOL
00588 couenne_->getRecordBestSol()->update(tmpSolution, n,
00589 obj, couenne_->getFeasTol());
00590 CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution);
00591 obj = couenne_->getRecordBestSol()->getVal();
00592 #endif
00593 #endif
00594
00595 if (babInfo){
00596 babInfo->setNlpSolution (tmpSolution, n, obj);
00597 babInfo->setHasNlpSolution (true);
00598 }
00599
00600 std::cout << "Final Nlp solution with objective " << obj << " :" << std::endl;
00601
00602 if (obj < objectiveValue - COUENNE_EPS) {
00603 std::cout << "New incumbent found" << std::endl;
00604 const CouNumber
00605 *lb = solver -> getColLower (),
00606 *ub = solver -> getColUpper ();
00607
00608
00609
00610 for (int i=0; i < n; ++i, ++lb, ++ub) {
00611 CouNumber &t = tmpSolution [i];
00612 if (t < *lb) t = *lb;
00613 else if (t > *ub) t = *ub;
00614 }
00615
00616 couenne_ -> setCutOff (obj);
00617 foundSolution = true;
00618 objectiveValue = obj;
00619 CoinCopyN (tmpSolution, n, newSolution);
00620 cinlp_->setColLower(colLowerNlp_);
00621 cinlp_->setColUpper(colUpperNlp_);
00622
00623 curr_milp->restoreBaseModel(numInitialRows_);
00624 delete[] xprime;
00625 delete[] tmpSolution;
00626 if (boundsChanged){
00627 curr_milp->setColLower(colLower_);
00628 curr_milp->setColUpper(colUpper_);
00629 }
00630
00631 double elapsed = CoinCpuTime()-startTime_;
00632 std::cout << "IR: Heuristic: " << objectiveValue << std::endl;
00633 std::cout << "IR: Heuristic time: " << elapsed << std::endl;
00634 return foundSolution;
00635 }
00636 }
00637 cinlp_->setColLower(colLowerNlp_);
00638 cinlp_->setColUpper(colUpperNlp_);
00639 double avgBound;
00640 numIntAtBound = computeIntAtBound(xtilde, avgBound);
00641
00642 if (numIntAtBound >= 50 ||
00643 numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){
00644
00645
00646 avgBound = floor(avgBound + 0.5);
00647 std::cout << "Using reverse LB with rhs " << avgBound << std::endl;
00648 writeLB(cut, xtilde, 'G', avgBound);
00649 revlb_all.insert(cut);
00650 }
00651 else{
00652
00653 branchToCut(xtilde, curr_milp, previousBranches);
00654 boundsChanged = true;
00655 }
00656 }
00657 if (h <= outerLoop -2){
00658
00659
00660 Bonmin::OsiTMINLPInterface * nlp = dynamic_cast<Bonmin::OsiTMINLPInterface *> (nlp_);
00661 Ipopt::SmartPtr< Ipopt::OptionsList > opt = nlp->solver()->options();
00662 nlp->messageHandler()->setLogLevel(10);
00663 double mu_target;
00664 double kappa_d;
00665 opt->GetNumericValue("mu_target", mu_target, "ipopt.");
00666 opt->SetNumericValue("mu_target", omega_*(h+1), "ipopt.");
00667 opt->GetNumericValue("kappa_d", kappa_d, "ipopt.");
00668 opt->SetNumericValue("kappa_d", 0.0, "ipopt.");
00669 nlp_->resolve();
00670 if (nlp_->isProvenPrimalInfeasible() ||
00671 nlp_->isProvenDualInfeasible() ||
00672 nlp_->isAbandoned()){
00673 nlp_->resolveForRobustness(3);
00674 }
00675 opt->SetNumericValue("mu_target", mu_target, "ipopt.");
00676 opt->SetNumericValue("kappa_d", kappa_d, "ipopt.");
00677
00678 if (!nlp->isProvenPrimalInfeasible() &&
00679 !nlp->isProvenDualInfeasible() &&
00680 !nlp->isAbandoned()){
00681 CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime);
00682 couenne_ -> getAuxs (xprime);
00683 curr_milp->restoreBaseModel(numInitialRows_);
00684 if (boundsChanged){
00685 curr_milp->setColLower(colLower_);
00686 curr_milp->setColUpper(colUpper_);
00687 }
00688 }
00689 }
00690 }
00691
00692 curr_milp->restoreBaseModel(numInitialRows_);
00693 delete[] xprime;
00694 delete[] tmpSolution;
00695 if (boundsChanged){
00696 curr_milp->setColLower(colLower_);
00697 curr_milp->setColUpper(colUpper_);
00698 }
00699 double elapsed = CoinCpuTime()-startTime_;
00700 std::cout << "IR: Heuristic: " << COUENNE_INFINITY << std::endl;
00701 std::cout << "IR: Heuristic time: " << elapsed << std::endl;
00702 return foundSolution;
00703 }
00704
00705 int
00706 CouenneIterativeRounding::improvementIR(double& objectiveValue,
00707 double* newSolution,
00708 const double* solution){
00709 std::cout << "starting Improvement IR" << std::endl;
00710
00711 OsiSolverInterface * solver = model_->solver();
00712
00713 OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
00714 Bonmin::BabInfo * babInfo = dynamic_cast<Bonmin::BabInfo *> (auxInfo);
00715
00716 if(babInfo){
00717 babInfo->setHasNlpSolution(false);
00718 if(babInfo->infeasibleNode()){
00719 return 0;
00720 }
00721 }
00722
00723
00724 int n = couenne_->nVars();
00725 int nNlp = cinlp_->getNumCols();
00726
00727
00728 int numIntAtBound = 0;
00729
00730 double lbrhs = CoinMin(baseLbRhs_, CoinMax(1,numIntegers_/2) );
00731 OsiRowCut cut;
00732 OsiRowCut lbcut1;
00733 OsiRowCut lbcut2;
00734 OsiCuts lbcuts;
00735 int currentIndex = 0;
00736 double obj;
00737
00738 bool boundsChanged = false;
00739 std::vector<int> previousBranches;
00740
00741 double avgBound = 0.0;
00742 numIntAtBound = computeIntAtBound(solution, avgBound);
00743
00744 if (numIntAtBound >= 50 ||
00745 numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){
00746
00747 writeLB(lbcut1, solution, 'L', lbrhs + floor(avgBound - 0.5));
00748 lbcuts.insert(lbcut1);
00749 writeLB(lbcut2, solution, 'G', 1);
00750 lbcuts.insert(lbcut2);
00751
00752
00753
00754
00755 nlp_->applyCuts(lbcuts);
00756 }
00757 else{
00758
00759
00760
00761 for (int i = 0; i < nlp_->getNumCols(); ++i){
00762 if (nlp_->isInteger(i)){
00763 nlp_->setColLower(i, colLowerNlp_[i]+(solution[i]-colLower_[i])*0.5);
00764 nlp_->setColUpper(i, colUpperNlp_[i]+(solution[i]-colUpper_[i])*0.5);
00765 }
00766 }
00767 }
00768
00769 nlp_->setColSolution(solution);
00770 nlp_->resolve();
00771 if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() ||
00772 nlp_->isAbandoned()){
00773 nlp_->resolveForRobustness(3);
00774 }
00775
00776 if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() ||
00777 nlp_->isAbandoned()){
00778 obj = COIN_DBL_MAX;
00779 }
00780 else{
00781 obj = nlp_->getObjValue();
00782 }
00783
00784
00785 if (numIntAtBound > 0){
00786 currentIndex = nlp_->getNumRows()-1;
00787 nlp_->deleteRows(1, ¤tIndex);
00788 currentIndex = nlp_->getNumRows()-1;
00789 nlp_->deleteRows(1, ¤tIndex);
00790 }
00791 else{
00792 nlp_->setColLower(colLowerNlp_);
00793 nlp_->setColUpper(colUpperNlp_);
00794 }
00795
00796 if (obj == COIN_DBL_MAX || obj >= objectiveValue - COUENNE_EPS){
00797 std::cout << " IR: no improvement possible " << std::endl;
00798 std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl;
00799 return 0;
00800 }
00801
00802 double* xprime = new double [n];
00803 CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime);
00804 couenne_ -> getAuxs (xprime);
00805
00806
00807
00808
00809 int constrInd[2];
00810 double constrElem[2];
00811
00812
00813 bool foundSolution = false;
00814 double* tmpSolution = new double[n];
00815 OsiCuts revlb_all;
00816 OsiSolverInterface* curr_milp = milp_;
00817
00818 OsiCuts cs;
00819 cs.insert(lbcut1);
00820 cs.insert(lbcut2);
00821
00822
00823
00824 for (int i = 0; i < n; i++){
00825 constrInd[0] = i;
00826 constrInd[1] = i+n;
00827 constrElem[0] = -1;
00828 constrElem[1] = -1;
00829 cut.mutableRow().setVector(2, constrInd, constrElem);
00830 cut.setLb(-COIN_DBL_MAX);
00831 cut.setUb(-xprime[i]);
00832 cs.insert(cut);
00833 constrElem[0] = +1;
00834 constrElem[1] = -1;
00835 cut.mutableRow().setVector(2, constrInd, constrElem);
00836 cut.setLb(-COIN_DBL_MAX);
00837 cut.setUb(xprime[i]);
00838 cs.insert(cut);
00839 }
00840 curr_milp->applyCuts(cs);
00841
00842 for (int k = 0; k < maxRoundingIter_ &&
00843 ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++k){
00844
00845 curr_milp->restoreBaseModel(numInitialRows_+cs.sizeRowCuts());
00846 curr_milp->applyCuts(revlb_all);
00847 bool solFound = solveMilp(curr_milp,
00848 endTime_-(CoinCpuTime()-startTime_)-2);
00849 if (!solFound && !boundsChanged){
00850 std::cout << " MILP cannot be solved, terminating LB " << std::endl;
00851
00852 curr_milp->restoreBaseModel(numInitialRows_);
00853 delete[] xprime;
00854 delete[] tmpSolution;
00855 if (boundsChanged){
00856 curr_milp->setColLower(colLower_);
00857 curr_milp->setColUpper(colUpper_);
00858 }
00859 std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl;
00860 return foundSolution;
00861 }
00862 else if (!solFound && boundsChanged){
00863
00864
00865
00866 curr_milp->setColLower(colLower_);
00867 curr_milp->setColUpper(colUpper_);
00868 branchToCut(tmpSolution, curr_milp, previousBranches);
00869 continue;
00870 }
00871 const double * xtilde = curr_milp->getColSolution();
00872
00873
00874
00875 CoinCopyN (xtilde, n, tmpSolution);
00876 for (int i = 0; i < nNlp; ++i){
00877 if (model_->isInteger(i)){
00878 tmpSolution[i] = floor(tmpSolution[i]+0.5);
00879 cinlp_->setColLower(i, tmpSolution[i]);
00880 cinlp_->setColUpper(i, tmpSolution[i]);
00881 }
00882 }
00883 cinlp_->setColSolution(tmpSolution);
00884
00885 cinlp_->messageHandler()->setLogLevel(1);
00886 cinlp_->resolve();
00887 obj = ((cinlp_->isProvenOptimal()) ? cinlp_->getObjValue():COIN_DBL_MAX);
00888 memcpy(tmpSolution, cinlp_->getColSolution(),
00889 nNlp*sizeof(double));
00890
00891
00892
00893 bool isChecked = false;
00894 #ifdef FM_CHECKNLP2
00895 isChecked = couenne_->checkNLP2(tmpSolution, 0, false,
00896 true,
00897 false,
00898 couenne_->getFeasTol());
00899 if(isChecked) {
00900 obj = couenne_->getRecordBestSol()->getModSolVal();
00901 }
00902 #else
00903 isChecked = couenne_->checkNLP(tmpSolution, obj, true);
00904 #endif
00905
00906 if (cinlp_->isProvenOptimal () &&
00907 isChecked &&
00908 (obj < couenne_->getCutOff())) {
00909
00910 #ifdef FM_CHECKNLP2
00911 #ifdef FM_TRACE_OPTSOL
00912 couenne_->getRecordBestSol()->update();
00913 CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution);
00914 obj = couenne_->getRecordBestSol()->getVal();
00915 #else
00916 CoinCopyN (couenne_->getRecordBestSol()->getModSol(n), n, tmpSolution);
00917 #endif
00918 #else
00919
00920
00921 couenne_ -> getAuxs (tmpSolution);
00922
00923 #ifdef FM_TRACE_OPTSOL
00924 couenne_->getRecordBestSol()->update(tmpSolution, n,
00925 obj, couenne_->getFeasTol());
00926 CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution);
00927 obj = couenne_->getRecordBestSol()->getVal();
00928 #endif
00929 #endif
00930
00931 if (babInfo){
00932 babInfo->setNlpSolution (tmpSolution, n, obj);
00933 babInfo->setHasNlpSolution (true);
00934 }
00935
00936 std::cout << "Final Nlp solution with objective " << obj << " :" << std::endl;
00937
00938 if (obj < objectiveValue - COUENNE_EPS) {
00939 std::cout << "New incumbent found" << std::endl;
00940 const CouNumber
00941 *lb = solver -> getColLower (),
00942 *ub = solver -> getColUpper ();
00943
00944
00945
00946 for (int i=0; i < n; ++i, ++lb, ++ub) {
00947 CouNumber &t = tmpSolution [i];
00948 if (t < *lb) t = *lb;
00949 else if (t > *ub) t = *ub;
00950 }
00951
00952 couenne_ -> setCutOff (obj);
00953 foundSolution = true;
00954 objectiveValue = obj;
00955 CoinCopyN (tmpSolution, n, newSolution);
00956 cinlp_->setColLower(colLowerNlp_);
00957 cinlp_->setColUpper(colUpperNlp_);
00958
00959 curr_milp->restoreBaseModel(numInitialRows_);
00960 delete[] xprime;
00961 delete[] tmpSolution;
00962 if (boundsChanged){
00963 curr_milp->setColLower(colLower_);
00964 curr_milp->setColUpper(colUpper_);
00965 }
00966
00967 double elapsed = CoinCpuTime()-startTime_;
00968 std::cout << "IR: Heuristic: " << objectiveValue << std::endl;
00969 std::cout << "IR: Heuristic time: " << elapsed << std::endl;
00970 return foundSolution;
00971 }
00972 }
00973 cinlp_->setColLower(colLowerNlp_);
00974 cinlp_->setColUpper(colUpperNlp_);
00975 numIntAtBound = computeIntAtBound(xtilde, avgBound);
00976
00977 if (numIntAtBound >= 50 ||
00978 numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){
00979
00980
00981 avgBound = floor(avgBound + 0.5);
00982 std::cout << "Using reverse LB with rhs " << avgBound << std::endl;
00983 writeLB(cut, xtilde, 'G', avgBound);
00984 revlb_all.insert(cut);
00985 }
00986 else{
00987
00988 branchToCut(xtilde, curr_milp, previousBranches);
00989 boundsChanged = true;
00990 }
00991 }
00992
00993 curr_milp->restoreBaseModel(numInitialRows_);
00994 delete[] xprime;
00995 delete[] tmpSolution;
00996 if (boundsChanged){
00997 curr_milp->setColLower(colLower_);
00998 curr_milp->setColUpper(colUpper_);
00999 }
01000 double elapsed = CoinCpuTime()-startTime_;
01001 std::cout << "IR: Heuristic: " << COUENNE_INFINITY << std::endl;
01002 std::cout << "IR: Heuristic time: " << elapsed << std::endl;
01003 return foundSolution;
01004 }
01005
01006 int
01007 CouenneIterativeRounding::computeIntAtBound(const double* x){
01008 int numIntAtBound = 0;
01009 for (int i = 0; i < nlp_->getNumCols(); ++i){
01010 if (nlp_->isInteger(i) && (areEqual(x[i], colLower_[i]) ||
01011 areEqual(x[i], colUpper_[i]))){
01012 numIntAtBound++;
01013 }
01014 }
01015 return numIntAtBound;
01016 }
01017
01018 int
01019 CouenneIterativeRounding::computeIntAtBound(const double* x,
01020 double& avgBoundSize){
01021 int numIntAtBound = 0;
01022 avgBoundSize = 0;
01023 for (int i = 0; i < nlp_->getNumCols(); ++i){
01024 if (nlp_->isInteger(i) && (areEqual(x[i], colLower_[i]) ||
01025 areEqual(x[i], colUpper_[i]))){
01026 numIntAtBound++;
01027 avgBoundSize += colUpper_[i] - colLower_[i];
01028 }
01029 }
01030 avgBoundSize /= numIntAtBound;
01031 return numIntAtBound;
01032 }
01033
01034 void
01035 CouenneIterativeRounding::writeLB(OsiRowCut& cut, const double* x,
01036 char sense, double rhs){
01037 cut.mutableRow().clear();
01038 for (int i = 0; i < nlp_->getNumCols(); ++i){
01039 if (nlp_->isInteger(i)){
01040 if (areEqual(x[i], colUpper_[i])){
01041 cut.mutableRow().insert(i, -1);
01042 rhs -= x[i];
01043 }
01044 else if (areEqual(x[i], colLower_[i])){
01045 cut.mutableRow().insert(i, 1);
01046 rhs += x[i];
01047 }
01048 }
01049 }
01050 if (sense == 'L'){
01051 cut.setLb(-COIN_DBL_MAX);
01052 cut.setUb(rhs);
01053 }
01054 else if (sense == 'G'){
01055 cut.setLb(rhs);
01056 cut.setUb(COIN_DBL_MAX);
01057 }
01058 else{
01059 std::cerr << "### ERROR: wrong sense of LB constraint" << std::endl;
01060 exit(1);
01061 }
01062 }
01063
01064 bool
01065 CouenneIterativeRounding::solveMilp(OsiSolverInterface* milp,
01066 double maxTime){
01067 double start = CoinCpuTime();
01068 #ifdef COIN_HAS_CPX
01069 OsiCpxSolverInterface * solver = dynamic_cast<OsiCpxSolverInterface*>(milp);
01070 CPXENVptr env = solver->getEnvironmentPtr();
01071 CPXLPptr cpxlp = solver->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
01072 int status = 0;
01073 bool solFound = false;
01074 bool infeasible = false;
01075 while (!solFound && !infeasible && ((CoinCpuTime() - start) < maxTime)){
01076 solver->branchAndBound();
01077 status = CPXgetstat(env, cpxlp);
01078 solFound = ((status == CPXMIP_NODE_LIM_FEAS) ||
01079 (status == CPXMIP_TIME_LIM_FEAS) ||
01080 (status == CPXMIP_MEM_LIM_FEAS) ||
01081 (status == CPXMIP_ABORT_FEAS) ||
01082 (status == CPXMIP_OPTIMAL) ||
01083 (status == CPXMIP_OPTIMAL_TOL) ||
01084 (status == CPXMIP_ABORT_FEAS) ||
01085 (status == CPXMIP_FAIL_FEAS) ||
01086 (status == CPXMIP_FAIL_FEAS_NO_TREE) ||
01087 (status == CPXMIP_FEASIBLE));
01088 infeasible = ((status == CPXMIP_INForUNBD) ||
01089 (status == CPXMIP_OPTIMAL_INFEAS) ||
01090 (status == CPXMIP_INFEASIBLE));
01091 }
01092 if (solFound){
01093 return true;
01094 }
01095 return false;
01096 #else
01097 CbcModel cbcModel(*milp);
01098 for (int i = 0; i < numHeuristics_; ++i){
01099 cbcModel.addHeuristic(heuristics_[i]);
01100 }
01101 cbcModel.setMaximumSeconds(CBCMILPTIME);
01102
01103 while ((cbcModel.getSolutionCount() == 0) &&
01104 (!cbcModel.isProvenInfeasible()) &&
01105 (!cbcModel.isProvenDualInfeasible()) &&
01106 (!cbcModel.isAbandoned()) &&
01107 ((CoinCpuTime() - start) < maxTime)){
01108 cbcModel.branchAndBound();
01109 }
01110
01111 milp = cbcModel.solver();
01112 if (cbcModel.getSolutionCount() > 0){
01113 return true;
01114 }
01115 return false;
01116 #endif
01117 }
01118
01119 int
01120 CouenneIterativeRounding::branchToCut(const double* x,
01121 OsiSolverInterface* solver,
01122 std::vector<int>& previousBranches){
01123 int branch;
01124 bool found = false;
01125 while (!found){
01126 branch = rand()%numIntegers_;
01127 found = true;
01128 for (unsigned int i = 0; i < previousBranches.size(); ++i){
01129 if (branch == previousBranches[i]){
01130 found = false;
01131 break;
01132 }
01133 }
01134 if (found){
01135 previousBranches.push_back(branch);
01136 }
01137 else{
01138 continue;
01139 }
01140 for (int i = 0; i < nlp_->getNumCols(); ++i){
01141 if (model_->isInteger(i)){
01142 if (branch == 0){
01143 branch = i;
01144 break;
01145 }
01146 else{
01147 branch--;
01148 }
01149 }
01150 }
01151 }
01152 double sample = rand();
01153 if (sample <= ((x[branch]-colLower_[branch])/(colUpper_[branch]-colLower_[branch]))){
01154 solver->setColUpper(branch, x[branch]-1);
01155 }
01156 else{
01157 solver->setColLower(branch, x[branch]+1);
01158 }
01159 return branch;
01160 }
01161
01162 void
01163 CouenneIterativeRounding::registerOptions (Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions){
01164 roptions -> AddStringOption2
01165 ("iterative_rounding_heuristic",
01166 "Do we use the Iterative Rounding heuristic",
01167 "no",
01168 "no","",
01169 "yes","",
01170 "If enabled, a heuristic based on Iterative Rounding is used "
01171 "to find feasible solutions for the problem. "
01172 "The heuristic may take some time, but usually finds good solutions. "
01173 "Recommended if you want good upper bounds and have Cplex. "
01174 "Not recommended if you do not have Cplex");
01175
01176 roptions -> AddNumberOption
01177 ("iterative_rounding_time",
01178 "Specify the maximum time allowed for the Iterative Rounding heuristic",
01179 -1, "Maximum CPU time employed by the Iterative Rounding heuristic; "
01180 "if no solution found in this time, failure is reported. "
01181 "This overrides the CPU time set by Aggressiveness if positive.");
01182
01183 roptions -> AddNumberOption
01184 ("iterative_rounding_time_firstcall",
01185 "Specify the maximum time allowed for the Iterative Rounding heuristic "
01186 "when no feasible solution is known",
01187 -1, "Maximum CPU time employed by the Iterative Rounding heuristic "
01188 "when no solution is known; if no solution found in this time, "
01189 "failure is reported."
01190 "This overrides the CPU time set by Aggressiveness if posive.");
01191
01192 roptions -> AddBoundedIntegerOption
01193 ("iterative_rounding_aggressiveness",
01194 "Aggressiveness of the Iterative Rounding heuristic",
01195 0, 2, 1,
01196 "Set the aggressiveness of the heuristic; i.e., how many iterations "
01197 "should be run, and with which parameters. The maximum time can be "
01198 "overridden by setting the _time and _time_firstcall options. "
01199 "0 = non aggressive, 1 = standard (default), 2 = aggressive.");
01200
01201 roptions -> AddLowerBoundedIntegerOption
01202 ("iterative_rounding_num_fir_points",
01203 "Max number of points rounded at the beginning of Iterative Rounding",
01204 1, 5,
01205 "Number of different points (obtained solving a log-barrier problem) "
01206 "that the heuristic will try to round at most, during its execution "
01207 "at the root node (i.e. the F-IR heuristic). Default 5.");
01208
01209 roptions -> AddBoundedNumberOption
01210 ("iterative_rounding_omega",
01211 "Omega parameter of the Iterative Rounding heuristic",
01212 0, true, 1, true, 0.2,
01213 "Set the omega parameter of the heuristic, which represents a "
01214 "multiplicative factor for the minimum log-barrier parameter "
01215 "of the NLP which is solved to obtain feasible points. This "
01216 "corresponds to $\\omega'$ in the paper. Default 0.2.");
01217
01218 roptions -> AddLowerBoundedIntegerOption
01219 ("iterative_rounding_base_lbrhs",
01220 "Base rhs of the local branching constraint for Iterative Rounding",
01221 0, 15,
01222 "Base rhs for the local branching constraint that defines a "
01223 "neighbourhood of the local incumbent. The base rhs is modified by "
01224 "the algorithm according to variable bounds. This corresponds to "
01225 "k' in the paper. Default 15.");
01226
01227 }
01228
01229 }