00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "CouenneAggrProbing.hpp"
00012 #include "CouenneProblemElem.hpp"
00013 #include "CouenneExprVar.hpp"
00014 #include "CouenneExprOpp.hpp"
00015
00016 #include "CouenneBab.hpp"
00017 #include "CouenneCutGenerator.hpp"
00018 #include <string>
00019
00020 #define COUENNE_AGGR_PROBING_FINITE_BOUND 1.0e+10
00021 #define COUENNE_AGGR_PROBING_MIN_INTERVAL 1.0e-2
00022 #define COUENNE_AGGR_PROBING_BND_RELAX COUENNE_EPS
00023
00024 using namespace Couenne;
00025
00026 CouenneAggrProbing::CouenneAggrProbing(CouenneSetup *setup,
00027 const Ipopt::SmartPtr<Ipopt::OptionsList> options)
00028 {
00029 couenne_ = setup;
00030 numCols_ = couenne_->couennePtr()->Problem()->nVars();
00031 maxTime_ = COIN_DBL_MAX;
00032 maxFailedSteps_ = 10;
00033 maxNodes_ = 1000;
00034 initCutoff_ = COUENNE_INFINITY;
00035 restoreCutoff_ = false;
00036
00037 }
00038
00039 CouenneAggrProbing::CouenneAggrProbing(const CouenneAggrProbing &rhs){
00040 couenne_ = new CouenneSetup(*rhs.couenne_);
00041 numCols_ = rhs.numCols_;
00042 maxTime_ = rhs.maxTime_;
00043 maxFailedSteps_ = rhs.maxFailedSteps_;
00044 maxNodes_ = rhs.maxNodes_;
00045 initCutoff_ = rhs.initCutoff_;
00046 restoreCutoff_ = rhs.restoreCutoff_;
00047 }
00048
00049 CouenneAggrProbing::~CouenneAggrProbing(){
00050 }
00051
00052 void CouenneAggrProbing::registerOptions(Ipopt::SmartPtr <Bonmin::RegisteredOptions> roptions) {
00053
00054 }
00055
00056 double CouenneAggrProbing::getMaxTime() const {
00057 return maxTime_;
00058 }
00059
00060 void CouenneAggrProbing::setMaxTime(double value){
00061 maxTime_ = value;
00062 }
00063
00064 int CouenneAggrProbing::getMaxFailedSteps() const {
00065 return maxFailedSteps_;
00066 }
00067
00068 void CouenneAggrProbing::setMaxFailedSteps(int value){
00069 maxFailedSteps_ = value;
00070 }
00071
00072 int CouenneAggrProbing::getMaxNodes() const {
00073 return maxNodes_;
00074 }
00075
00076 void CouenneAggrProbing::setMaxNodes(int value){
00077 maxNodes_ = value;
00078 }
00079
00080 bool CouenneAggrProbing::getRestoreCutoff() const {
00081 return restoreCutoff_;
00082 }
00083
00084 void CouenneAggrProbing::setRestoreCutoff(bool value){
00085 restoreCutoff_ = value;
00086 }
00087
00088 double CouenneAggrProbing::probeVariable(int index, bool probeLower){
00089
00090
00091 OsiSolverInterface* nlp = couenne_->nonlinearSolver();
00092 OsiSolverInterface* lp = couenne_->continuousSolver();
00093 CouenneProblem* problem = couenne_->couennePtr()->Problem();
00094
00095
00096 double initUpper = lp->getColUpper()[index];
00097 double initLower = lp->getColLower()[index];
00098
00099 double* initLowerLp = new double[numCols_];
00100 double* initUpperLp = new double[numCols_];
00101
00102 memcpy(initLowerLp, lp->getColLower(), numCols_*sizeof(double));
00103 memcpy(initUpperLp, lp->getColUpper(), numCols_*sizeof(double));
00104
00105 if (initUpper < initLower + COUENNE_EPS){
00106
00107 return ((probeLower) ? initLower : initUpper);
00108 }
00109
00110
00111 int indobj = problem->Obj(0)->Body()->Index();
00112
00113
00114 double initCutoff = problem->Ub()[indobj];
00115
00116 double* initCutoffSol = NULL;
00117
00118 if (restoreCutoff_ && problem->getCutOff() < COUENNE_INFINITY){
00119 initCutoffSol = new double[numCols_];
00120 memcpy(initCutoffSol, problem->getCutOffSol(), numCols_*sizeof(double));
00121 }
00122
00123
00124 Bonmin::BabSetupBase::NodeComparison initNodeComparison =
00125 couenne_->nodeComparisonMethod();
00126 int initMaxNodes = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxNodes);
00127 double initMaxTime = couenne_->getDoubleParameter(Bonmin::BabSetupBase::MaxTime);
00128 int initMaxSol = couenne_->getIntParameter(Bonmin::BabSetupBase::MaxSolutions);
00129 couenne_->setNodeComparisonMethod(Bonmin::BabSetupBase::bestBound);
00130
00131 couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, maxNodes_);
00132 couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, COIN_INT_MAX);
00133 problem->setCheckAuxBounds(true);
00134
00136 Bonmin::BabSetupBase::HeuristicMethods heuristics = couenne_->heuristics();
00137 couenne_->heuristics().clear();
00138
00139 double currentBound = (probeLower) ? initLower : initUpper;
00140 double startTime = CoinCpuTime();
00141 int failedSteps = 0;
00142 double intervalSize = 0.0;
00143 double tryBound = 0.0;
00144
00145 int iter = 0;
00146
00147 if (probeLower)
00148 std::cout << "Probing lower on var " << index << std::endl;
00149 else
00150 std::cout << "Probing upper on var " << index << std::endl;
00151
00152 if ((fabs(currentBound) > COUENNE_AGGR_PROBING_FINITE_BOUND) &&
00153 ((probeLower && initUpper > -COUENNE_AGGR_PROBING_FINITE_BOUND) ||
00154 (!probeLower && initLower < COUENNE_AGGR_PROBING_FINITE_BOUND))){
00155
00156
00157
00158
00159 if (probeLower){
00160 tryBound = -COUENNE_AGGR_PROBING_FINITE_BOUND;
00161 lp->setColLower(index, currentBound);
00162 problem->Lb()[index] = currentBound;
00163 lp->setColUpper(index, tryBound);
00164 problem->Ub()[index] = tryBound;
00165 if (index < problem->nOrigVars()){
00166 nlp->setColLower(index, currentBound);
00167 nlp->setColUpper(index, tryBound);
00168 }
00169 }
00170 else{
00171 tryBound = COUENNE_AGGR_PROBING_FINITE_BOUND;
00172 lp->setColLower(index, tryBound);
00173 problem->Lb()[index] = tryBound;
00174 lp->setColUpper(index, currentBound);
00175 problem->Ub()[index] = currentBound;
00176 if (index < problem->nOrigVars()){
00177 nlp->setColLower(index, tryBound);
00178 nlp->setColUpper(index, currentBound);
00179 }
00180 }
00181
00183 couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime,
00184 CoinMin(maxTime_-(CoinCpuTime()-startTime),
00185 maxTime_*0.5));
00186
00187 if (restoreCutoff_){
00188 problem->resetCutOff(initCutoff);
00189 problem->Ub()[indobj] = initCutoff;
00190 problem->installCutOff();
00191 }
00192
00193 std::cout << "Iteration " << iter << ", current bound " << currentBound
00194 << ", try bound " << tryBound << std::endl;
00195
00197
00198
00199
00200 CouenneBab bb;
00201
00202 bb(couenne_);
00203 if (bb.model().isProvenInfeasible()){
00205 currentBound = tryBound;
00206 std::cout << "Probing succeeded; brought to finite" << std::endl;
00207 }
00208 else{
00210 std::cout << "Probing failed; still infinity, exit" << std::endl;
00211 }
00212 iter++;
00213 }
00214
00215
00216
00217
00218
00219 intervalSize = 0.1;
00220
00221 if (intervalSize < COUENNE_AGGR_PROBING_MIN_INTERVAL){
00222 intervalSize = COUENNE_AGGR_PROBING_MIN_INTERVAL;
00223 }
00224
00225 while ((fabs(currentBound) <= COUENNE_AGGR_PROBING_FINITE_BOUND) &&
00226 ((CoinCpuTime() - startTime) < maxTime_) &&
00227 (failedSteps < maxFailedSteps_) &&
00228 (intervalSize >= COUENNE_AGGR_PROBING_MIN_INTERVAL) &&
00229 iter < 100){
00230
00232 if (probeLower){
00233 tryBound = currentBound + intervalSize;
00234 if (tryBound > initUpper){
00235
00236
00237 tryBound = initUpper;
00238 }
00239 if (lp->isInteger(index)){
00240 tryBound = floor(tryBound);
00241 }
00242
00243 lp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
00244 problem->Lb()[index] = currentBound - COUENNE_AGGR_PROBING_BND_RELAX;
00245 lp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
00246 problem->Ub()[index] = tryBound + COUENNE_AGGR_PROBING_BND_RELAX;
00247 if (index < problem->nOrigVars()){
00248 nlp->setColLower(index, currentBound - COUENNE_AGGR_PROBING_BND_RELAX);
00249 nlp->setColUpper(index, tryBound + COUENNE_AGGR_PROBING_BND_RELAX);
00250 }
00251 }
00252 else{
00253 tryBound = currentBound - intervalSize;
00254 if (tryBound < initLower){
00255
00256
00257 tryBound = initLower;
00258 }
00259 if (lp->isInteger(index)){
00260 tryBound = ceil(tryBound);
00261 }
00262
00263 lp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
00264 problem->Lb()[index] = tryBound - COUENNE_AGGR_PROBING_BND_RELAX;
00265 lp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
00266 problem->Ub()[index] = currentBound + COUENNE_AGGR_PROBING_BND_RELAX;
00267 if (index < problem->nOrigVars()){
00268 nlp->setColLower(index, tryBound - COUENNE_AGGR_PROBING_BND_RELAX);
00269 nlp->setColUpper(index, currentBound + COUENNE_AGGR_PROBING_BND_RELAX);
00270 }
00271 }
00272
00273 lp->resolve();
00274 problem->domain()->push(numCols_, lp->getColSolution(),
00275 lp->getColLower(), lp->getColUpper());
00276
00278 couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime,
00279 CoinMin(maxTime_-(CoinCpuTime()-startTime),
00280 maxTime_*0.5));
00281
00282 if (restoreCutoff_){
00283 problem->Ub()[indobj] = initCutoff;
00284 problem->resetCutOff(initCutoff);
00285 problem->installCutOff();
00286 }
00287
00288 std::cout << "Iteration " << iter << ", current bound " << currentBound
00289 << ", try bound " << tryBound << std::endl;
00290
00292
00293
00294
00295 CouenneBab bb;
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) ||
00306 (restoreCutoff_ &&
00307 problem->getCutOffSol() &&
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
00347
00348 if ((probeLower && fabs(currentBound-initUpper) < COUENNE_EPS) ||
00349 (!probeLower && fabs(currentBound-initLower) < COUENNE_EPS)){
00350 failedSteps = maxFailedSteps_;
00351 }
00352
00353
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
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
00402
00403
00404
00405 OsiSolverInterface* lp = couenne_->continuousSolver();
00406 CouenneProblem* problem = couenne_->couennePtr()->Problem();
00407
00408
00409 double initUpper = lp->getColUpper()[index];
00410 double initLower = lp->getColLower()[index];
00411
00412 if (initUpper < initLower + COUENNE_EPS){
00413
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
00444 expression* extraVar = NULL;
00445
00446 lp->writeLp("before");
00447
00448 if (probeLower){
00449 std::cout << "Probing LOWER" << std::endl;
00450
00451 newLpObj[index] = 1.0;
00452 lp->setObjective(newLpObj);
00453
00454 lp->writeLp("lower");
00455
00456
00457 problem->Obj(0)->Body(problem->Variables()[index]);
00458
00459
00460
00461
00462 }
00463 else{
00464
00465
00466
00467
00468
00469 int extraCol = numCols_;
00470 lp->setObjective(newLpObj);
00471 lp->addCol(0, NULL, NULL, -initUpper, -initLower, 1.0);
00472
00473
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
00480 extraVar = problem->addVariable(lp->isInteger(index), NULL);
00481
00482
00483 problem->Obj(0)->Body(extraVar);
00484
00485
00486
00487
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
00498
00499
00500 CouenneBab bb;
00501 bb(couenne_);
00502
00503 double bestBound = bb.model().getBestPossibleObjValue();
00504
00505 std::cout << "Obtained bound: " << bb.model().getBestPossibleObjValue() << std::endl;
00506
00507
00509 couenne_->setNodeComparisonMethod (initNodeComparison);
00510 couenne_->setIntParameter(Bonmin::BabSetupBase::MaxNodes, initMaxNodes);
00511 couenne_->setDoubleParameter(Bonmin::BabSetupBase::MaxTime, initMaxTime);
00512 couenne_->setIntParameter(Bonmin::BabSetupBase::MaxSolutions, initMaxSol);
00513 couenne_->heuristics() = heuristics;
00514
00515 if (!probeLower){
00516 int extra = lp->getNumCols()-1;
00517 lp->deleteCols(1, &extra);
00518 extra = lp->getNumRows()-1;
00519 lp->deleteRows(1, &extra);
00520 problem->Variables().pop_back();
00521 delete extraVar;
00522
00523 }
00524
00525 lp->setObjective(initLpObj);
00526 problem->Obj(0)->Body(initProbObj);
00527
00528 delete[] initLpObj;
00529 delete[] newLpObj;
00530
00531 return ((probeLower) ? bestBound : -bestBound);
00532
00533 }