/home/coin/SVN-release/OS-2.4.0/Couenne/src/bound_tightening/CouenneAggrProbing.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneAggrProbing.cpp 530 2011-03-12 16:52:06Z pbelotti $
00002  *
00003  * Name:    CouenneAggrProbing.cpp
00004  * Author:  Giacomo Nannicini
00005  * Purpose: Aggressive probing
00006  *
00007  * (C) Giacomo Nannicini, 2010.
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CouenneAggrProbing.hpp"
00012 #include "CouenneProblemElem.hpp"
00013 #include "CouenneExprVar.hpp"
00014 #include "CouenneExprOpp.hpp"
00015 #include "BonCbc.hpp"
00016 #include "CouenneCutGenerator.hpp"
00017 #include <string>
00018 
00019 #define COUENNE_AGGR_PROBING_FINITE_BOUND 1.0e+10
00020 #define COUENNE_AGGR_PROBING_MIN_INTERVAL 1.0e-2
00021 #define COUENNE_AGGR_PROBING_BND_RELAX COUENNE_EPS
00022 
00023 using namespace Couenne;
00024 
00025 CouenneAggrProbing::CouenneAggrProbing(CouenneSetup *setup,
00026                                        const Ipopt::SmartPtr<Ipopt::OptionsList> options)
00027 {
00028   couenne_ = setup;
00029   numCols_ = couenne_->couennePtr()->Problem()->nVars();
00030   maxTime_ = COIN_DBL_MAX;
00031   maxFailedSteps_ = 10;
00032   maxNodes_ = 1000;
00033   initCutoff_ = COUENNE_INFINITY;
00034   restoreCutoff_ = false;
00035 
00036 }
00037 
00038 CouenneAggrProbing::CouenneAggrProbing(const CouenneAggrProbing &rhs){
00039   couenne_ = new CouenneSetup(*rhs.couenne_);
00040   numCols_ = rhs.numCols_;
00041   maxTime_ = rhs.maxTime_;
00042   maxFailedSteps_ = rhs.maxFailedSteps_;
00043   maxNodes_ = rhs.maxNodes_;
00044   initCutoff_ = rhs.initCutoff_;
00045   restoreCutoff_ = rhs.restoreCutoff_;
00046 }
00047 
00048 CouenneAggrProbing::~CouenneAggrProbing(){
00049 }
00050 
00051 void CouenneAggrProbing::registerOptions(Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
00052   // Nothing for the moment, but will be added later as needed
00053 }
00054 
00055 void CouenneAggrProbing::generateCuts(const OsiSolverInterface& solver,
00056                                       OsiCuts& cuts, 
00057                                       const CglTreeInfo info) const {
00058   // Empty for the moment (cannot be used automatically in Branch-and-Bound)
00059 }
00060 
00061 double CouenneAggrProbing::getMaxTime() const {
00062   return maxTime_;
00063 }
00064 
00065 void CouenneAggrProbing::setMaxTime(double value){
00066   maxTime_ = value;
00067 }
00068 
00069 int CouenneAggrProbing::getMaxFailedSteps() const {
00070   return maxFailedSteps_;
00071 }
00072 
00073 void CouenneAggrProbing::setMaxFailedSteps(int value){
00074   maxFailedSteps_ = value;
00075 }
00076 
00077 int CouenneAggrProbing::getMaxNodes() const {
00078   return maxNodes_;
00079 }
00080 
00081 void CouenneAggrProbing::setMaxNodes(int value){
00082   maxNodes_ = value;
00083 }
00084 
00085 bool CouenneAggrProbing::getRestoreCutoff() const {
00086   return restoreCutoff_;
00087 }
00088 
00089 void CouenneAggrProbing::setRestoreCutoff(bool value){
00090   restoreCutoff_ = value;
00091 }
00092 
00093 double CouenneAggrProbing::probeVariable(int index, bool probeLower){
00094 
00095   // Useful objects for easy access
00096   OsiSolverInterface* nlp = couenne_->nonlinearSolver();
00097   OsiSolverInterface* lp = couenne_->continuousSolver();
00098   CouenneProblem* problem = couenne_->couennePtr()->Problem();
00099 
00100   // Save initial bounds
00101   double initUpper = lp->getColUpper()[index];
00102   double initLower = lp->getColLower()[index];
00103 
00104   double* initLowerLp = new double[numCols_];
00105   double* initUpperLp = new double[numCols_];
00106 
00107   memcpy(initLowerLp, lp->getColLower(), numCols_*sizeof(double));
00108   memcpy(initUpperLp, lp->getColUpper(), numCols_*sizeof(double));
00109 
00110   if (initUpper < initLower + COUENNE_EPS){
00111     // Variable is fixed, so we can't tighten
00112     return ((probeLower) ? initLower : initUpper);
00113   }
00114 
00115   // Index of the aux variable representing the objective function
00116   int indobj = problem->Obj(0)->Body()->Index();
00117 
00118   // Initial cutoff value
00119   double initCutoff = problem->Ub()[indobj];
00120 
00121   double* initCutoffSol = NULL; 
00122 
00123   if (restoreCutoff_ && problem->getCutOff() < COUENNE_INFINITY){
00124     initCutoffSol = new double[numCols_];
00125     memcpy(initCutoffSol, problem->getCutOffSol(), numCols_*sizeof(double));
00126   }
00127 
00128   // Save parameters
00129   Bonmin::BabSetupBase::NodeComparison initNodeComparison = 
00130     couenne_->nodeComparisonMethod();
00131   int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
00132   double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
00133   int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
00134   couenne_->setNodeComparisonMethod(Bonmin::BabSetupBase::bestBound);
00135   //couenne_->nodeComparisonMethod() = Bonmin::BabSetupBase::bestBound;
00136   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
00137   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
00138   problem->setCheckAuxBounds(true);
00139 
00141   Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
00142   couenne_->heuristics().clear();
00143 
00144   double currentBound = (probeLower) ? initLower : initUpper;
00145   double startTime = CoinCpuTime();
00146   int failedSteps = 0;
00147   double intervalSize = 0.0;
00148   double tryBound = 0.0;
00149 
00150   int iter = 0;
00151 
00152   if (probeLower)
00153     std::cout << "Probing lower on var " << index << std::endl;
00154   else
00155     std::cout << "Probing upper on var " << index << std::endl;
00156 
00157   if ((fabs(currentBound) > COUENNE_AGGR_PROBING_FINITE_BOUND) &&
00158       ((probeLower && initUpper > -COUENNE_AGGR_PROBING_FINITE_BOUND) ||
00159        (!probeLower && initLower < COUENNE_AGGR_PROBING_FINITE_BOUND))){
00160     // The bound is too large to apply the standard probing method;
00161     // try to reduce it to a finite value. We only do this if we want
00162     // to probe a variable on an infinite (or close to) bound, and the
00163     // other bound of the variable is sufficiently far away
00164     if (probeLower){
00165       tryBound = -COUENNE_AGGR_PROBING_FINITE_BOUND;
00166       lp->setColLower(index, currentBound);
00167       problem->Lb()[index] = currentBound;
00168       lp->setColUpper(index, tryBound);
00169       problem->Ub()[index] = tryBound;
00170       if (index < problem->nOrigVars()){
00171         nlp->setColLower(index, currentBound);
00172         nlp->setColUpper(index, tryBound);
00173       }
00174     }
00175     else{
00176       tryBound = COUENNE_AGGR_PROBING_FINITE_BOUND;
00177       lp->setColLower(index, tryBound);
00178       problem->Lb()[index] = tryBound;
00179       lp->setColUpper(index, currentBound);
00180       problem->Ub()[index] = currentBound;
00181       if (index < problem->nOrigVars()){
00182         nlp->setColLower(index, tryBound);
00183         nlp->setColUpper(index, currentBound);
00184       }
00185     }
00186 
00188     couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
00189                                  CoinMin(maxTime_-(CoinCpuTime()-startTime),
00190                                          maxTime_*0.5));
00191 
00192     if (restoreCutoff_){
00193       problem->resetCutOff(initCutoff);
00194       problem->Ub()[indobj] = initCutoff;
00195       problem->installCutOff();
00196     }
00197 
00198     std::cout << "Iteration " << iter << ", current bound " << currentBound
00199               << ", try bound " << tryBound << std::endl;
00200 
00202     Bonmin::Bab bb;
00203     bb.setUsingCouenne(true);
00204     bb(couenne_);
00205     if (bb.model().isProvenInfeasible()){
00207       currentBound = tryBound;
00208       std::cout << "Probing succeeded; brought to finite" << std::endl;
00209     }
00210     else{
00212       std::cout << "Probing failed; still infinity, exit" << std::endl;
00213     }    
00214     iter++;
00215   }
00216 
00217   // Now that we have a finite bound, pick size of the probing interval
00218 
00219   // Override (for testing - will be chosen automatically in final
00220   // implementation)
00221   intervalSize = 0.1;
00222 
00223   if (intervalSize < COUENNE_AGGR_PROBING_MIN_INTERVAL){
00224     intervalSize = COUENNE_AGGR_PROBING_MIN_INTERVAL;
00225   }
00226   
00227   while ((fabs(currentBound) <= COUENNE_AGGR_PROBING_FINITE_BOUND) &&
00228          ((CoinCpuTime() - startTime) < maxTime_) &&
00229          (failedSteps < maxFailedSteps_) &&
00230          (intervalSize >= COUENNE_AGGR_PROBING_MIN_INTERVAL) && 
00231          iter < 100){
00232 
00234     if (probeLower){
00235       tryBound = currentBound + intervalSize;
00236       if (tryBound > initUpper){
00237         // It does not make sense to use bounds larger than the initial
00238         // ones
00239         tryBound = initUpper;
00240       }
00241       if (lp->isInteger(index)){
00242         tryBound = floor(tryBound);
00243       }
00244       // Relax bounds a little bit
00245       lp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
00246       problem->Lb()[index] = currentBound - COUENNE_AGGR_PROBING_BND_RELAX;
00247       lp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
00248       problem->Ub()[index] = tryBound + COUENNE_AGGR_PROBING_BND_RELAX;
00249       if (index < problem->nOrigVars()){
00250         nlp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
00251         nlp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
00252       }
00253     }
00254     else{
00255       tryBound = currentBound - intervalSize;
00256       if (tryBound < initLower){
00257         // It does not make sense to use bounds larger than the initial
00258         // ones
00259         tryBound = initLower;
00260       }
00261       if (lp->isInteger(index)){
00262         tryBound = ceil(tryBound);
00263       }
00264       // Relax bounds a little bit
00265       lp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
00266       problem->Lb()[index] = tryBound - COUENNE_AGGR_PROBING_BND_RELAX;
00267       lp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
00268       problem->Ub()[index] = currentBound + COUENNE_AGGR_PROBING_BND_RELAX;
00269       if (index < problem->nOrigVars()){
00270         nlp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
00271         nlp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
00272       }
00273     }
00274 
00275     lp->resolve();
00276     problem->domain()->push(numCols_, lp->getColSolution(),
00277                             lp->getColLower(), lp->getColUpper());
00278 
00280     couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
00281                                  CoinMin(maxTime_-(CoinCpuTime()-startTime),
00282                                          maxTime_*0.5));
00283 
00284     if (restoreCutoff_){
00285       problem->Ub()[indobj] = initCutoff;
00286       problem->resetCutOff(initCutoff);
00287       problem->installCutOff();
00288     }
00289 
00290     std::cout << "Iteration " << iter << ", current bound " << currentBound
00291               << ", try bound " << tryBound << std::endl;
00292 
00294     Bonmin::Bab bb;
00295     bb.setUsingCouenne(true);
00296     bb(couenne_);
00297 
00298     problem->domain()->pop();
00299 
00300     double obj = 0.0;
00302     bool intervalSearched = (bb.model().isProvenOptimal() || 
00303                              bb.model().isProvenInfeasible());
00304 
00305     if ((!intervalSearched) || // If the search is not complete
00306         (restoreCutoff_ && // or we don't want to accept new solutions
00307          problem->getCutOffSol() && // and we have a new feasible solution
00308          problem->checkNLP(problem->getCutOffSol(), obj, true))){
00310       if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
00312         failedSteps = maxFailedSteps_;
00313       }
00314       else{
00315         intervalSize /= 2;
00316       }
00317       failedSteps++;
00318       std::cout << "Probing failed; shrinking interval" << std::endl;
00319     }
00320     else{
00325       if (lp->isInteger(index) && fabs(tryBound-currentBound) < 0.5){
00328         intervalSize = 1.0;
00329       }
00330       else{
00331         intervalSize *= 2;
00332       }
00333       currentBound = tryBound;
00334       if (lp->isInteger(index)){
00335         if (probeLower){
00336           currentBound += 1.0;
00337         }
00338         else {
00339           currentBound -= 1.0;
00340         }
00341       }
00342       failedSteps = 0;
00343       std::cout << "Probing succeeded; enlarging interval" << std::endl;
00344     }
00345 
00346     // Check early termination condition: if we manage to fix the
00347     // variable (unlikely), there is nothing more we can do
00348     if ((probeLower && fabs(currentBound-initUpper) < COUENNE_EPS) ||
00349         (!probeLower && fabs(currentBound-initLower) < COUENNE_EPS)){
00350       failedSteps = maxFailedSteps_;
00351     }
00352 
00353     // Reset cutoff
00354     if (restoreCutoff_){
00355       problem->Ub()[indobj] = initCutoff;
00356       problem->resetCutOff(initCutoff);
00357       problem->installCutOff();
00358     }
00359 
00360     problem->domain()->pop();
00361 
00362     iter++;
00363   }
00364 
00367   lp->setColLower(initLowerLp);
00368   lp->setColUpper(initUpperLp);
00369   nlp->setColLower(initLowerLp);
00370   nlp->setColUpper(initUpperLp);
00371   memcpy(problem->Lb(), initLowerLp, numCols_*sizeof(double));
00372   memcpy(problem->Ub(), initUpperLp, numCols_*sizeof(double));
00373 
00375   problem->setCheckAuxBounds(false);
00376   //couenne_->nodeComparisonMethod() = initNodeComparison;
00377   couenne_->setNodeComparisonMethod(initNodeComparison);
00378   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
00379   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
00380   couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
00381   couenne_->heuristics() = heuristics;
00382 
00384   if (restoreCutoff_){
00385     problem->resetCutOff();
00386     problem->setCutOff(initCutoff, initCutoffSol);
00387     if (initCutoffSol){
00388       delete[] initCutoffSol;
00389     }
00390   }
00391   
00392   delete[] initLowerLp;
00393   delete[] initUpperLp;
00394 
00396   return currentBound;
00397 
00398 }
00399 
00400 double CouenneAggrProbing::probeVariable2(int index, bool probeLower){
00401   // Does not work! It's impossible to get Maximization problems working...
00402   // Adding extra variables doesn't seem to do the trick
00403 
00404   // Useful objects for easy access
00405   OsiSolverInterface* lp = couenne_->continuousSolver();
00406   CouenneProblem* problem = couenne_->couennePtr()->Problem();
00407 
00408   // Save initial bounds
00409   double initUpper = lp->getColUpper()[index];
00410   double initLower = lp->getColLower()[index];
00411 
00412   if (initUpper < initLower + COUENNE_EPS){
00413     // Variable is fixed, so we can't probe anything
00414     return ((probeLower) ? initLower : initUpper);
00415   }
00416 
00420   Bonmin::BabSetupBase::NodeComparison initNodeComparison = 
00421     couenne_->nodeComparisonMethod();
00422   int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
00423   double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
00424   int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
00425   couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
00426   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
00427   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
00428 
00432   Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
00433   couenne_->heuristics().clear();
00434 
00436   double* initLpObj = new double[numCols_];
00437   memcpy(initLpObj, lp->getObjCoefficients(), numCols_*sizeof(double));
00438   expression* initProbObj = problem->Obj(0)->Body();
00439 
00440   double* newLpObj = new double[numCols_];
00441   memset(newLpObj, 0, numCols_*sizeof(double));
00442 
00443   //expression* exprObj = NULL; // PB: unused
00444   expression* extraVar = NULL;
00445 
00446   lp->writeLp("before");
00447 
00448   if (probeLower){
00449     std::cout << "Probing LOWER" << std::endl;
00450     // Set LP objective function
00451     newLpObj[index] = 1.0;
00452     lp->setObjective(newLpObj);
00453 
00454     lp->writeLp("lower");
00455 
00456     // Set CouenneProblem objective
00457     problem->Obj(0)->Body(problem->Variables()[index]);
00458     // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
00459     // problem->setCutOff(initUpper, lp->getColUpper());
00460     // problem->installCutOff();
00461 
00462   }
00463   else{
00464     // We cannot maximize an objective function in Couenne; so, we
00465     // have to introduce an additional variable, equal to the opposite
00466     // of the variable that we want to maximize, and minimize that.
00467 
00468     // Add one column and one row to the LP
00469     int extraCol = numCols_;
00470     lp->setObjective(newLpObj);
00471     lp->addCol(0, NULL, NULL, -initUpper, -initLower, 1.0);
00472 
00473     // Create the row x_{extraCol} = -x_{index}
00474     int rowIndices[2] = {index, extraCol};
00475     double rowElements[2] = {1.0, 1.0};
00476     lp->addRow(2, rowIndices, rowElements, 0.0, 0.0);
00477     lp->resolve();
00478 
00479     // Create the expression x_{extraCol} = -x_{index} in CouenneProblem
00480     extraVar = problem->addVariable(lp->isInteger(index), NULL);
00481     // exprObj = new exprOpp(problem->Variables()[index]->clone());
00482     // problem->addEQConstraint(extraVar, exprObj);
00483     problem->Obj(0)->Body(extraVar);
00484 
00485     // couenne_->setDoubleParameter(Bonmin::BabSetupBase::Cutoff, COIN_DBL_MAX);
00486     // problem->setCutOff(-initLower, lp->getColLower());
00487     // problem->installCutOff();
00488 
00489     lp->writeLp("upper");
00490   }
00491 
00492   couenne_->setNodeComparisonMethod (Bonmin::BabSetupBase::bestBound);
00493   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
00494   couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, 
00495                                maxTime_);
00496 
00497   Bonmin::Bab bb;
00498   bb.setUsingCouenne(true);
00499   bb(couenne_);
00500 
00501   double bestBound = bb.model().getBestPossibleObjValue();
00502 
00503   std::cout << "Obtained bound: " << bb.model().getBestPossibleObjValue() << std::endl;
00504 
00505 
00507   couenne_->setNodeComparisonMethod (initNodeComparison);
00508   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
00509   couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
00510   couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
00511   couenne_->heuristics() = heuristics;
00512 
00513   if (!probeLower){
00514     int extra = lp->getNumCols()-1;
00515     lp->deleteCols(1, &extra);
00516     extra = lp->getNumRows()-1;
00517     lp->deleteRows(1, &extra);
00518     problem->Variables().pop_back();
00519     delete extraVar;
00520     // delete exprObj;
00521   }
00522 
00523   lp->setObjective(initLpObj);
00524   problem->Obj(0)->Body(initProbObj);
00525 
00526   delete[] initLpObj;
00527   delete[] newLpObj;
00528 
00529   return ((probeLower) ? bestBound : -bestBound);
00530 
00531 }

Generated on Thu Sep 22 03:05:55 2011 by  doxygen 1.4.7