00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BonOACutGenerator2.hpp"
00011 #include "BonminConfig.h"
00012
00013 #include "OsiClpSolverInterface.hpp"
00014
00015 #include "CbcModel.hpp"
00016 #include "BonCbcLpStrategy.hpp"
00017 #ifdef COIN_HAS_CPX
00018 #include "OsiCpxSolverInterface.hpp"
00019 #endif
00020 #include "OsiAuxInfo.hpp"
00021 #include "BonSolverHelp.hpp"
00022
00023 #include <climits>
00024
00025 namespace Bonmin
00026 {
00027 static const char * txt_id = "OA decomposition";
00028
00029
00031 OACutGenerator2::OACutGenerator2(BabSetupBase & b):
00032 OaDecompositionBase(b, true, false)
00033 {
00034 std::string bonmin="bonmin.";
00035 std::string prefix = (b.prefix() == bonmin) ? "" : b.prefix();
00036 prefix += "oa_decomposition.";
00037 subMip_ = new SubMipSolver(b, prefix);
00038 double oaTime;
00039 b.options()->GetNumericValue("time_limit",oaTime,prefix);
00040 parameter().maxLocalSearch_ = INT_MAX;
00041 b.options()->GetIntegerValue("solution_limit", parameter().maxSols_,prefix);
00042 parameter().maxLocalSearchTime_ =
00043 std::min(b.getDoubleParameter(BabSetupBase::MaxTime), oaTime);
00044 if(parameter().maxSols_ > b.getIntParameter(BabSetupBase::MaxSolutions))
00045 parameter().maxSols_ = b.getIntParameter(BabSetupBase::MaxSolutions);
00046 }
00047 OACutGenerator2::~OACutGenerator2()
00048 {
00049 delete subMip_;
00050 }
00051
00053 bool
00054 OACutGenerator2::doLocalSearch(BabInfo * babInfo) const
00055 {
00056 return (nLocalSearch_<parameters_.maxLocalSearch_ &&
00057 numSols_ < parameters_.maxSols_ &&
00058 CoinCpuTime() - timeBegin_ < parameters_.maxLocalSearchTime_);
00059 }
00061 double
00062 OACutGenerator2::performOa(OsiCuts &cs,
00063 solverManip &lpManip,
00064 BabInfo * babInfo,
00065 double & cutoff, const CglTreeInfo & info) const
00066 {
00067
00068 double lastPeriodicLog = CoinCpuTime();
00069
00070
00071 double gap_tol = this->parameter().gap_tol_;
00072
00073 bool isInteger = false;
00074
00075 subMip_->setLpSolver(lpManip.si());
00076 OsiSolverInterface * lp = subMip_->solver();
00077 lp->resolve();
00078
00079 if(IsValid(nlp_->linearizer())){
00080 nlp_->linearizer()->get_refined_oa(cs);
00081 installCuts(*lp, cs, cs.sizeRowCuts());
00082 }
00083 lp->resolve();
00084
00085 OsiBranchingInformation branch_info(lp, false);
00086 bool milpOptimal = 1;
00087
00088
00089 double milpBound = -COIN_DBL_MAX;
00090
00091 bool feasible = 1;
00092
00093 subMip_->solve(cutoff, parameters_.subMilpLogLevel_,
00094 parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
00095
00096 milpBound = std::max(milpBound, subMip_->lowBound());
00097 milpOptimal = subMip_->optimal();
00098
00099 feasible = milpBound < cutoff;
00100
00101 isInteger = (subMip_->getLastSolution() != NULL);
00102 nLocalSearch_++;
00103
00104 if (milpOptimal)
00105 handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip_->nodeCount()<<subMip_->iterationCount()<<CoinMessageEol;
00106 else
00107 {
00108 handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip_->nodeCount()<<subMip_->iterationCount()<<CoinMessageEol;
00109 }
00110
00111 int nodeCount = subMip_->nodeCount();
00112
00113 int numberPasses = 0;
00114
00115 #ifdef OA_DEBUG
00116 bool foundSolution = 0;
00117 #endif
00118 double * nlpSol = NULL;
00119 double ub = cutoff;
00120 double gap = 1;
00121 while (isInteger && feasible) {
00122 numberPasses++;
00123
00124 double time = CoinCpuTime();
00125 if (time - lastPeriodicLog > parameters_.logFrequency_) {
00126 handler_->message(PERIODIC_MSG,messages_)
00127 <<time - timeBegin_<<numberPasses<<cutoff
00128 <<milpBound
00129 <<CoinMessageEol;
00130 lastPeriodicLog = CoinCpuTime();
00131 }
00132
00133
00134
00135 int numberCutsBefore = cs.sizeRowCuts();
00136
00137
00138 const double * colsol =
00139 subMip_->getLastSolution();
00140 branch_info.solution_ = colsol;
00141
00142 fixIntegers(*nlp_,branch_info, parameters_.cbcIntegerTolerance_,objects_, nObjects_);
00143
00144 nlp_->resolve(txt_id);
00145 if (post_nlp_solve(babInfo, cutoff)) {
00146
00147
00148 ub = std::min(nlp_->getObjValue(), ub);
00149 cutoff = ub > 0 ? ub *(1 - parameters_.cbcCutoffIncrement_) : ub*(1 + parameters_.cbcCutoffIncrement_);
00150 assert(cutoff < ub);
00151
00152 lp->setDblParam(OsiDualObjectiveLimit, cutoff);
00153 numSols_++;
00154 }
00155
00156 nlpSol = const_cast<double *>(nlp_->getColSolution());
00157
00158
00159 if(parameter().addOnlyViolated_){
00160 nlp_->getOuterApproximation(cs, nlpSol, 1, colsol,
00161 parameter().global_);
00162 }
00163 else {
00164 nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
00165 parameter().global_);
00166 }
00167
00168 int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
00169 assert(numberCuts);
00170 installCuts(*lp, cs, numberCuts);
00171
00172 lp->resolve();
00173
00174 double objvalue = lp->getObjValue();
00175
00176 feasible = (lp->isProvenOptimal() &&
00177 !lp->isDualObjectiveLimitReached() && (objvalue<cutoff)) ;
00178
00179 bool changed = !feasible;
00180 branch_info.solution_ = lp->getColSolution();
00181 if (!changed)
00182 changed = isDifferentOnIntegers(*nlp_, objects_, nObjects_,
00183 parameters_.cbcIntegerTolerance_,
00184 nlp_->getColSolution(), lp->getColSolution());
00185 if (changed) {
00186
00187 isInteger = integerFeasible(*lp, branch_info, parameters_.cbcIntegerTolerance_,
00188 objects_, nObjects_);
00189 }
00190 else {
00191 isInteger = 0;
00192
00193 milpBound = 1e200;
00194 }
00195 #ifdef OA_DEBUG
00196 printf("Obj value after cuts %g %d rows\n",lp->getObjValue(),
00197 numberCuts) ;
00198 #endif
00199
00200 if (CoinCpuTime() - timeBegin_ > parameters_.maxLocalSearchTime_)
00201 break;
00202
00203 if(ub!=0)
00204 gap = (ub - milpBound)/fabs(ub);
00205 else
00206 gap = -milpBound/(1e-10);
00207 if (gap < gap_tol){
00208 milpBound = 1e50;
00209 feasible = 0;
00210 }
00211
00212 if (feasible &&
00213 nLocalSearch_ < parameters_.maxLocalSearch_ &&
00214 numSols_ < parameters_.maxSols_) {
00215
00216 if(ub!=0)
00217 gap = (ub - milpBound)/fabs(ub);
00218 else
00219 gap = -milpBound/(1e-10);
00220 if (gap < gap_tol){
00221 milpBound = 1e50;
00222 feasible = 0;
00223 }
00224 nLocalSearch_++;
00225
00226 subMip_->solve(cutoff, parameters_.subMilpLogLevel_,
00227 parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
00228
00229 milpBound = std::max(milpBound, subMip_->lowBound());
00230
00231 if (subMip_->optimal())
00232 handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip_->nodeCount()<<subMip_->iterationCount()<<CoinMessageEol;
00233 else
00234 handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip_->nodeCount()<<subMip_->iterationCount()<<CoinMessageEol;
00235
00236 nodeCount += subMip_->nodeCount();
00237
00238 colsol = const_cast<double *> (subMip_->getLastSolution());
00239 isInteger = (colsol != 0);
00240
00241 feasible = (milpBound < cutoff);
00242
00243 if (feasible && isInteger) {
00244 bool changed = false;
00245 changed = isDifferentOnIntegers(*nlp_, objects_, nObjects_,
00246 0.1,
00247 nlp_->getColSolution(), colsol);
00248
00249 if (!changed) {
00250 feasible = 0;
00251 milpBound = 1e50;
00252 }
00253
00254 }
00255 if (subMip_->optimal())
00256 milpOptimal = 1;
00257 else {
00258 milpOptimal = 0;
00259 }
00260 if (milpBound < cutoff)
00261 handler_->message(UPDATE_LB, messages_)<<milpBound<<CoinCpuTime() - timeBegin_<<CoinMessageEol;
00262 else {
00263 milpBound = 1e50;
00264 feasible = 0;
00265 }
00266 }
00267
00268 }
00269
00270 if(milpBound >= cutoff){
00271 handler_->message(OASUCCESS, messages_)<<"OA "<<CoinCpuTime() - timeBegin_
00272 <<ub<<milpBound<<CoinMessageEol;
00273 }
00274 else {
00275 handler_->message(OAABORT, messages_)<<"OA "<<CoinCpuTime() - timeBegin_
00276 <<ub<<milpBound<<CoinMessageEol;
00277 }
00278 handler_->message(OA_STATS, messages_)<<numberPasses<<nodeCount
00279 <<CoinMessageEol;
00280 #ifdef OA_DEBUG
00281 debug_.printEndOfProcedureDebugMessage(cs, foundSolution, cutoff, milpBound, isInteger, feasible, std::cout);
00282 #endif
00283 return milpBound;
00284 }
00285
00287 void
00288 OACutGenerator2::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00289 {
00290 roptions->SetRegisteringCategory("Outer Approximation Decomposition (B-OA)", RegisteredOptions::BonminCategory);
00291 roptions->AddStringOption2("oa_decomposition", "If yes do initial OA decomposition",
00292 "no",
00293 "no","",
00294 "yes","",
00295 "");
00296 roptions->setOptionExtraInfo("oa_decomposition",19);
00297
00298 roptions->SetRegisteringCategory("Output and Loglevel", RegisteredOptions::BonminCategory);
00299 roptions->AddBoundedIntegerOption("oa_log_level",
00300 "specify OA iterations log level.",
00301 0,2,1,
00302 "Set the level of output of OA decomposition solver : "
00303 "0 - none, 1 - normal, 2 - verbose"
00304 );
00305 roptions->setOptionExtraInfo("oa_log_level", 25);
00306
00307 roptions->AddLowerBoundedNumberOption("oa_log_frequency",
00308 "display an update on lower and upper bounds in OA every n seconds",
00309 0.,1.,100.,
00310 "");
00311 roptions->setOptionExtraInfo("oa_log_frequency", 25);
00312 }
00313 }
00314