00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <sstream>
00011 #include <climits>
00012
00013 #include <algorithm>
00014 #include "BonOaDecBase.hpp"
00015
00016
00017 #include "BonminConfig.h"
00018
00019 #include "OsiClpSolverInterface.hpp"
00020
00021 #include "CbcModel.hpp"
00022 #include "CbcStrategy.hpp"
00023 #ifdef COIN_HAS_CPX
00024 #include "OsiCpxSolverInterface.hpp"
00025
00026 #define CHECK_CPX_STAT(a,b) if(b) throw CoinError("Error in CPLEX call",__FILE__,a);
00027
00028 #endif
00029 #include "BonCbc.hpp"
00030 #include "BonSolverHelp.hpp"
00031
00032 extern CbcModel * OAModel;
00033
00034 namespace Bonmin {
00035
00036 OaDecompositionBase::OaDecompositionBase
00037 (OsiTMINLPInterface * nlp)
00038 :
00039 CglCutGenerator(),
00040 nlp_(nlp),
00041 s_(NULL),
00042 nSolve_(0),
00043 lp_(NULL),
00044 objects_(NULL),
00045 nObjects_(0),
00046 nLocalSearch_(0),
00047 handler_(NULL),
00048 leaveSiUnchanged_(true),
00049 reassignLpsolver_(false),
00050 timeBegin_(0),
00051 parameters_(),
00052 currentNodeNumber_(-1)
00053 {
00054 handler_ = new CoinMessageHandler();
00055 handler_ -> setLogLevel(2);
00056 messages_ = OaMessages();
00057 timeBegin_ = CoinCpuTime();
00058 }
00059
00060 OaDecompositionBase::OaDecompositionBase(BabSetupBase &b, bool leaveSiUnchanged,
00061 bool reassignLpsolver):
00062 CglCutGenerator(),
00063 nlp_(b.nonlinearSolver()),
00064 lp_(NULL),
00065 objects_(NULL),
00066 nObjects_(0),
00067 nLocalSearch_(0),
00068 handler_(NULL),
00069 leaveSiUnchanged_(leaveSiUnchanged),
00070 reassignLpsolver_(reassignLpsolver),
00071 timeBegin_(0),
00072 numSols_(0),
00073 parameters_(),
00074 currentNodeNumber_(-1)
00075 {
00076 handler_ = new CoinMessageHandler();
00077 int logLevel;
00078 b.options()->GetIntegerValue("oa_log_level",logLevel,b.prefix());
00079 b.options()->GetNumericValue("oa_log_frequency",parameters_.logFrequency_,b.prefix());
00080
00081 handler_ -> setLogLevel(logLevel);
00082 b.options()->GetIntegerValue("solution_limit", parameters_.maxSols_,b.prefix());
00083
00084 messages_ = OaMessages();
00085 timeBegin_ = CoinCpuTime();
00086 b.options()->GetIntegerValue("milp_log_level",parameters_.subMilpLogLevel_,b.prefix());
00087 b.options()->GetNumericValue("cutoff_decr",parameters_.cbcCutoffIncrement_,b.prefix());
00088 b.options()->GetNumericValue("integer_tolerance",parameters_.cbcIntegerTolerance_,b.prefix());
00089 int ivalue;
00090 b.options()->GetEnumValue("add_only_violated_oa", ivalue,b.prefix());
00091 parameters_.addOnlyViolated_ = ivalue;
00092 b.options()->GetEnumValue("oa_cuts_scope", ivalue,b.prefix());
00093 parameters_.global_ = ivalue;
00094 }
00095
00096 OaDecompositionBase::OaDecompositionBase
00097 (const OaDecompositionBase & other)
00098 :
00099 CglCutGenerator(other),
00100 nlp_(other.nlp_),
00101 lp_(other.lp_),
00102 objects_(other.objects_),
00103 nObjects_(other.nObjects_),
00104 nLocalSearch_(0),
00105 messages_(other.messages_),
00106 leaveSiUnchanged_(other.leaveSiUnchanged_),
00107 reassignLpsolver_(other.reassignLpsolver_),
00108 timeBegin_(0),
00109 numSols_(other.numSols_),
00110 parameters_(other.parameters_),
00111 currentNodeNumber_(other.currentNodeNumber_)
00112 {
00113 timeBegin_ = CoinCpuTime();
00114 handler_ = other.handler_->clone();
00115 }
00117 OaDecompositionBase::Parameters::Parameters():
00118 global_(true),
00119 addOnlyViolated_(false),
00120 cbcCutoffIncrement_(1e-06),
00121 cbcIntegerTolerance_(1e-05),
00122 maxLocalSearch_(0),
00123 maxLocalSearchTime_(3600),
00124 subMilpLogLevel_(0),
00125 maxSols_(INT_MAX),
00126 logFrequency_(1000.),
00127 strategy_(NULL)
00128 {}
00129
00131 OaDecompositionBase::~OaDecompositionBase()
00132 {
00133 delete handler_;
00134 }
00135
00136
00138 OaDecompositionBase::Parameters::Parameters(const Parameters & other):
00139 global_(other.global_),
00140 addOnlyViolated_(other.addOnlyViolated_),
00141 cbcCutoffIncrement_(other.cbcCutoffIncrement_),
00142 cbcIntegerTolerance_(other.cbcIntegerTolerance_),
00143 maxLocalSearch_(other.maxLocalSearch_),
00144 maxLocalSearchTime_(other.maxLocalSearchTime_),
00145 subMilpLogLevel_(other.subMilpLogLevel_),
00146 maxSols_(other.maxSols_),
00147 logFrequency_(other.logFrequency_),
00148 strategy_(NULL)
00149 {
00150 if (other.strategy_)
00151 strategy_ = other.strategy_->clone();
00152 }
00153
00154
00155
00156 OaDecompositionBase::solverManip::solverManip
00157 (OsiSolverInterface * si,
00158 bool saveNumRows,
00159 bool saveBasis,
00160 bool saveBounds,
00161 bool saveCutoff,
00162 bool resolve):
00163 si_(si),
00164 initialNumberRows_(-1),
00165 colLower_(NULL),
00166 colUpper_(NULL),
00167 warm_(NULL),
00168 cutoff_(COIN_DBL_MAX),
00169 deleteSolver_(false),
00170 objects_(NULL),
00171 nObjects_(0)
00172 {
00173 getCached();
00174 if (saveNumRows)
00175 initialNumberRows_ = numrows_;
00176 if (saveBasis)
00177 warm_ = si->getWarmStart();
00178 if (saveBounds) {
00179 colLower_ = new double[numcols_];
00180 colUpper_ = new double[numcols_];
00181 CoinCopyN(si->getColLower(), numcols_ , colLower_);
00182 CoinCopyN(si->getColUpper(), numcols_ , colUpper_);
00183 }
00184 if (saveCutoff)
00185 si->getDblParam(OsiDualObjectiveLimit, cutoff_);
00186 si->messageHandler()->setLogLevel(0);
00187 if (resolve) si->resolve();
00188 }
00189
00190
00191 OaDecompositionBase::solverManip::solverManip
00192 (const OsiSolverInterface & si):
00193 si_(NULL),
00194 initialNumberRows_(-1),
00195 colLower_(NULL),
00196 colUpper_(NULL),
00197 warm_(NULL),
00198 cutoff_(COIN_DBL_MAX),
00199 deleteSolver_(true),
00200 objects_(NULL),
00201 nObjects_(0)
00202 {
00203 si_ = si.clone();
00204 getCached();
00205 }
00206
00207 OaDecompositionBase::solverManip::~solverManip()
00208 {
00209 if (warm_) delete warm_;
00210 if (colLower_) delete [] colLower_;
00211 if (colUpper_) delete [] colUpper_;
00212 if (deleteSolver_) delete si_;
00213 }
00214
00215 void
00216 OaDecompositionBase::solverManip::restore()
00217 {
00218 if (initialNumberRows_ >= 0) {
00219 int nRowsToDelete = si_->getNumRows() - initialNumberRows_;
00220 int * rowsToDelete = new int[nRowsToDelete];
00221 for (int i = 0 ; i < nRowsToDelete ; i++) {
00222 rowsToDelete[i] = i + initialNumberRows_;
00223 }
00224 si_->deleteRows(nRowsToDelete, rowsToDelete);
00225 delete [] rowsToDelete;
00226 numrows_ = si_->getNumRows() ;
00227 }
00228
00229 if (colLower_) {
00230 si_->setColLower(colLower_);
00231 }
00232
00233 if (colUpper_) {
00234 si_->setColUpper(colUpper_);
00235 }
00236
00237 if (cutoff_<COIN_DBL_MAX) {
00238 si_->setDblParam(OsiDualObjectiveLimit, cutoff_);
00239 }
00240
00241 if (warm_) {
00242 if (si_->setWarmStart(warm_)==false) {
00243 throw CoinError("Fail restoring the warm start at the end of procedure",
00244 "restore","OaDecompositionBase::SaveSolverState") ;
00245 }
00246 }
00247 getCached();
00248 }
00249
00250 void
00251 OaDecompositionBase::passInMessageHandler(CoinMessageHandler * handler)
00252 {
00253 int logLevel = handler_->logLevel();
00254 delete handler_;
00255 handler_=handler->clone();
00256 handler_->setLogLevel(logLevel);
00257 }
00258
00259
00260
00261 #if 1
00262
00263 void
00264 OaDecompositionBase::solverManip::cloneOther(const OsiSolverInterface &si){
00265
00266 int numberCutsToAdd = si.getNumRows();
00267 numberCutsToAdd -= numrows_;
00268 if (numberCutsToAdd > 0)
00269 {
00270 CoinPackedVector * * addCuts = new CoinPackedVector *[numberCutsToAdd];
00271 for (int i = 0 ; i < numberCutsToAdd ; i++)
00272 {
00273 addCuts[i] = new CoinPackedVector;
00274 }
00275
00276 const CoinPackedMatrix * mat = si.getMatrixByCol();
00277 const CoinBigIndex * start = mat->getVectorStarts();
00278 const int * length = mat->getVectorLengths();
00279 const double * elements = mat->getElements();
00280 const int * indices = mat->getIndices();
00281 for (int i = 0 ; i < numcols_ ; i++)
00282 for (int k = start[i] ; k < start[i] + length[i] ; k++)
00283 {
00284 if (indices[k] >= numrows_) {
00285 addCuts[ indices[k] - numrows_ ]->insert(i, elements[k]);
00286 }
00287 }
00288 si_->addRows(numberCutsToAdd, (const CoinPackedVectorBase * const *) addCuts, si.getRowLower() + numrows_,
00289 si.getRowUpper() + numrows_);
00290 for (int i = 0 ; i < numberCutsToAdd ; i++){
00291 delete addCuts[i];
00292 }
00293 delete [] addCuts;
00294 }
00295 else if (numberCutsToAdd < 0) {
00296 throw CoinError("Internal error number of cuts wrong",
00297 "generateCuts","OACutGenerator2");
00298 }
00299
00300 si_->setColLower(si.getColLower());
00301 si_->setColUpper(si.getColUpper());
00302
00303 CoinWarmStart * warm = si.getWarmStart();
00304 if (si_->setWarmStart(warm)==false) {
00305 delete warm;
00306 throw CoinError("Fail installing the warm start in the subproblem",
00307 "generateCuts","OACutGenerator2") ;
00308 }
00309 delete warm;
00310
00311 double cutoff;
00312 si.getDblParam(OsiDualObjectiveLimit, cutoff);
00313 si_->setDblParam(OsiDualObjectiveLimit, cutoff);
00314 si_->messageHandler()->setLogLevel(0);
00315 si_->resolve();
00316
00317 numrows_ = si.getNumRows();
00318 #ifdef OA_DEBUG
00319
00320 std::cout<<"Resolve with hotstart :"<<std::endl
00321 <<"number of iterations(should be 0) : "<<si_->getIterationCount()<<std::endl
00322 <<"Objective value and diff to original : "<<si_->getObjValue()<<", "
00323 <<fabs(si_->getObjValue() - si.getObjValue())<<std::endl;
00324 for (int i = 0 ; i <= numcols_ ; i++) {
00325 if (fabs(si.getColSolution()[i]-si_->getColSolution()[i])>1e-08) {
00326 std::cout<<"Diff between solution at node and solution with local solver : "<<fabs(si.getColSolution()[i]-si_->getColSolution()[i])<<std::endl;
00327 }
00328 }
00329 #endif
00330
00331 }
00332 #endif
00333
00334
00336 void
00337 OaDecompositionBase::generateCuts(const OsiSolverInterface &si, OsiCuts & cs,
00338 const CglTreeInfo info) const{
00339 if (nlp_ == NULL) {
00340 throw CoinError("Error in cut generator for outer approximation no NLP ipopt assigned", "generateCuts", "OaDecompositionBase");
00341 }
00342
00343
00344 BabInfo * babInfo = dynamic_cast<BabInfo *> (si.getAuxiliaryInfo());
00345 assert(babInfo);
00346 assert(babInfo->babPtr());
00347 numSols_ = babInfo->babPtr()->model().getSolutionCount ();
00348 CglTreeInfo info_copy = info;
00349 const CbcNode * node = babInfo->babPtr()->model().currentNode();
00350 info_copy.level = (node == NULL) ? 0 : babInfo->babPtr()->model().currentNode()->depth();
00351 if(babInfo->hasSolution()) numSols_ ++;
00352 if (babInfo)
00353 if (!babInfo->mipFeasible())
00354 return;
00355
00356
00357 const double *colsol = si.getColSolution();
00358
00359
00360 vector<double> savedColLower(nlp_->getNumCols());
00361 CoinCopyN(nlp_->getColLower(), nlp_->getNumCols(), savedColLower());
00362 vector<double> savedColUpper(nlp_->getNumCols());
00363 CoinCopyN(nlp_->getColUpper(), nlp_->getNumCols(), savedColUpper());
00364
00365
00366 OsiBranchingInformation brInfo(nlp_, false);
00367 brInfo.solution_ = colsol;
00368
00369 bool isInteger = integerFeasible(*nlp_, brInfo, parameters_.cbcIntegerTolerance_,
00370 objects_, nObjects_);
00371
00372 SubMipSolver * subMip = NULL;
00373
00374
00375 int nodeNumber = babInfo->babPtr()->model().getNodeCount();
00376 if(nodeNumber == currentNodeNumber_){
00377 #ifdef OA_DEBUG
00378 printf("OA decomposition recalled from the same node!\n");
00379 #endif
00380 int numCuts = savedCuts_.sizeRowCuts();
00381 for(int i = 0 ; i < numCuts ; i++){
00382
00383 if(savedCuts_.rowCut(i).violated(colsol) > 0.){
00384 #ifdef OA_DEBUG
00385 printf("A violated cut has been found\n");
00386 #endif
00387 savedCuts_.rowCut(i).setEffectiveness(9.99e99);
00388 cs.insert(savedCuts_.rowCut(i));
00389 savedCuts_.eraseRowCut(i);
00390 return;
00391 i--; numCuts--;
00392 }
00393 }
00394 }
00395 else {
00396 currentNodeNumber_ = nodeNumber;
00397 savedCuts_.dumpCuts();
00398 }
00399
00400 if (!isInteger) {
00401 if (!doLocalSearch(babInfo))
00402 return;
00403 }
00404
00405
00406 double cutoff;
00407 si.getDblParam(OsiDualObjectiveLimit, cutoff);
00408
00409
00410
00411 solverManip * lpManip = NULL;
00412 if (lp_ != NULL) {
00413 if (lp_!=&si) {
00414 #if 1
00415 lpManip = new solverManip(lp_, true, false, false, true, true);
00416 lpManip->cloneOther(si);
00417 #endif
00418 }
00419 else {
00420 #if 0
00421 throw CoinError("Not allowed to modify si in a cutGenerator",
00422 "OACutGenerator2","generateCuts");
00423 #else
00424 lpManip = new solverManip(lp_, true, leaveSiUnchanged_, true, true);
00425 #endif
00426 }
00427 }
00428 else {
00429 lpManip = new solverManip(si);
00430 }
00431 lpManip->setObjects(objects_, nObjects_);
00432 if (!isInteger) {
00433 subMip = new SubMipSolver(lpManip->si(), parameters_.strategy());
00434 }
00435
00436 double milpBound = performOa(cs, *lpManip, subMip, babInfo, cutoff, info_copy);
00437
00438 if(babInfo->hasSolution()){
00439 babInfo->babPtr()->model().setSolutionCount (numSols_ - 1);
00440 }
00441
00442
00443 {
00444 if (milpBound>-1e100)
00445 {
00446
00447 if (babInfo)
00448 babInfo->setMipBound(milpBound);
00449 }
00450 }
00451
00452
00453 if (subMip!= NULL) {
00454 delete subMip;
00455 subMip = NULL;
00456 }
00457
00458
00459 if (leaveSiUnchanged_)
00460 lpManip->restore();
00461 delete lpManip;
00462
00463 nlp_->setColLower(savedColLower());
00464 nlp_->setColUpper(savedColUpper());
00465
00466 return;
00467 }
00468
00469 void
00470 OaDecompositionBase::solverManip::getCached(){
00471 numrows_ = si_->getNumRows();
00472 numcols_ = si_->getNumCols();
00473 siColLower_ = si_->getColLower();
00474 siColUpper_ = si_->getColUpper();
00475 }
00476
00477
00479 bool
00480 OaDecompositionBase::post_nlp_solve(BabInfo * babInfo, double cutoff) const{
00481 nSolve_++;
00482 bool return_value = false;
00483 if (nlp_->isProvenOptimal()) {
00484 handler_->message(FEASIBLE_NLP, messages_)
00485 <<nlp_->getIterationCount()
00486 <<nlp_->getObjValue()<<CoinMessageEol;
00487
00488 #ifdef OA_DEBUG
00489 const double * colsol2 = nlp_->getColSolution();
00490 debug_.checkInteger(*nlp_,std::cerr);
00491 #endif
00492
00493 if ((nlp_->getObjValue() < cutoff) ) {
00494 handler_->message(UPDATE_UB, messages_)
00495 <<nlp_->getObjValue()
00496 <<CoinCpuTime()-timeBegin_
00497 <<CoinMessageEol;
00498
00499 return_value = true;
00500
00501 assert(babInfo);
00502 if (babInfo) {
00503 int numcols = nlp_->getNumCols();
00504 double * lpSolution = new double[numcols + 1];
00505 CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
00506 lpSolution[numcols] = nlp_->getObjValue();
00507 babInfo->setSolution(lpSolution,
00508 numcols + 1, lpSolution[numcols]);
00509 delete [] lpSolution;
00510 }
00511 }
00512 }
00513 else if (nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
00514 (*handler_)<<"Unsolved NLP... exit"<<CoinMessageEol;
00515 }
00516 else {
00517 handler_->message(INFEASIBLE_NLP, messages_)
00518 <<nlp_->getIterationCount()
00519 <<CoinMessageEol;
00520 }
00521 return return_value;
00522 }
00523
00524
00525
00526 #ifdef OA_DEBUG
00527 bool
00528 OaDecompositionBase::OaDebug::checkInteger(const OsiSolverInterface &nlp,
00529 std::ostream & os) const {
00530 const double * colsol = nlp.getColSolution();
00531 int numcols = nlp.getNumCols();
00532 for (int i = 0 ; i < numcols ; i++) {
00533 if (nlp.isInteger(i)) {
00534 if (fabs(colsol[i]) - floor(colsol[i] + 0.5) >
00535 1e-07) {
00536 std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
00537 <<" is, "<<fabs(colsol[i] - floor(colsol[i] + 0.5))<<std::endl;
00538 }
00539 }
00540 return true;
00541 }
00542
00543 }
00544
00545 void
00546 OaDecompositionBase::OaDebug::printEndOfProcedureDebugMessage(const OsiCuts &cs,
00547 bool foundSolution,
00548 double solValue,
00549 double milpBound,
00550 bool isInteger,
00551 bool feasible,
00552 std::ostream & os) const{
00553 std::cout<<"------------------------------------------------------------------"
00554 <<std::endl;
00555 std::cout<<"OA procedure finished"<<std::endl;
00556 std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
00557 if (foundSolution)
00558 std::cout <<"Found NLP-integer feasible solution of value : "<<solValue<<std::endl;
00559 std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
00560 std::cout<<"-------------------------------------------------------------------"<<std::endl;
00561 std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
00562
00563 }
00564
00565
00566
00567 #endif
00568 }
00569