/home/coin/SVN-release/OS-2.1.0/Bonmin/src/Algorithms/OaGenerators/BonFpForMinlp.cpp

Go to the documentation of this file.
00001 // (C) Copyright CNRS 2008
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // P. Bonami, CNRS
00007 //
00008 // Date : 02/13/2009
00009 
00010 #include "BonFpForMinlp.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 = "FP for MINLP";
00029   MinlpFeasPump::MinlpFeasPump(BabSetupBase & b):
00030       OaDecompositionBase(b, true, false)
00031   {
00032     int ivalue;
00033     std::string bonmin="bonmin.";
00034     std::string prefix = (b.prefix() == bonmin) ? "" : b.prefix();
00035     prefix += "pump_for_minlp.";
00036     b.options()->GetEnumValue("milp_solver",ivalue,prefix);
00037     if (ivalue <= 0) {//uses cbc
00038       CbcStrategyDefault strategy;
00039       setStrategy(strategy);
00040     }
00041     else if (ivalue == 1) {
00042       CbcStrategyChooseCuts strategy(b, prefix);
00043       setStrategy(strategy);
00044     }
00045     else if (ivalue == 2) {
00046 #ifdef COIN_HAS_CPX
00047       OsiCpxSolverInterface * cpxSolver = new OsiCpxSolverInterface;
00048       b.nonlinearSolver()->extractLinearRelaxation(*cpxSolver);
00049       assignLpInterface(cpxSolver);
00050 #else
00051       std::cerr << "You have set an option to use CPLEX as the milp\n"
00052       << "subsolver in oa decomposition. However, apparently\n"
00053       << "CPLEX is not configured to be used in bonmin.\n"
00054       << "See the manual for configuring CPLEX\n";
00055       throw -1;
00056 #endif
00057     }
00058 
00059     double oaTime;
00060     b.options()->GetNumericValue("time_limit",oaTime,prefix);
00061     parameter().maxLocalSearch_ = INT_MAX;
00062     b.options()->GetIntegerValue("solution_limit", parameter().maxSols_,prefix);
00063     parameter().maxLocalSearchTime_ =
00064     std::min(b.getDoubleParameter(BabSetupBase::MaxTime), oaTime);
00065     if(parameter().maxSols_ > b.getIntParameter(BabSetupBase::MaxSolutions))
00066       parameter().maxSols_ = b.getIntParameter(BabSetupBase::MaxSolutions);
00067        
00068   }
00069 
00070   MinlpFeasPump::~MinlpFeasPump()
00071   {}
00072 
00074   bool
00075   MinlpFeasPump::doLocalSearch(BabInfo * babInfo) const
00076   {
00077     return (nLocalSearch_<parameters_.maxLocalSearch_ &&
00078         CoinCpuTime() - timeBegin_ < parameters_.maxLocalSearchTime_ &&
00079         numSols_ < parameters_.maxSols_);
00080   }
00082   double
00083   MinlpFeasPump::performOa(OsiCuts &cs,
00084       solverManip &lpManip,
00085       SubMipSolver * &subMip,
00086       BabInfo * babInfo,
00087       double & cutoff,const CglTreeInfo &info) const
00088   {
00089 
00090     //bool interuptOnLimit = false;
00091     //double lastPeriodicLog = CoinCpuTime();
00092 
00093     const int numcols = nlp_->getNumCols();
00094     vector<double> savedColLower(nlp_->getNumCols());
00095     CoinCopyN(nlp_->getColLower(), nlp_->getNumCols(), savedColLower());
00096     vector<double> savedColUpper(nlp_->getNumCols());
00097     CoinCopyN(nlp_->getColUpper(), nlp_->getNumCols(), savedColUpper());
00098 
00099 
00100     OsiSolverInterface * lp = lpManip.si();
00101 
00102     vector<int> indices;
00103     for(int i = 0; i < numcols ; i++) {
00104       lp->setObjCoeff(i,0);
00105       if(!lp->isInteger(i)) {
00106       }
00107       else { indices.push_back(i);}
00108     }
00109 
00110     // If objective is linear need to add to lp constraint for objective
00111     const double * colsol = NULL;
00112     lp->resolve();
00113     OsiBranchingInformation branch_info(lp, false);
00114     branch_info.lower_ = savedColLower();
00115     branch_info.upper_ = savedColUpper();
00116     if(lp->getNumCols() == nlp_->getNumCols())
00117       nlp_->addObjectiveFunction(*lp, nlp_->getColSolution());
00118     lp->setObjCoeff(numcols,0);
00119 
00120     bool milpOptimal = false;
00121     nlp_->resolve(txt_id);
00122     if (subMip)//Perform a local search
00123     {
00124       assert(subMip->solver() == lp);
00125       set_fp_objective(*lp, nlp_->getColSolution());
00126       lp->initialSolve();
00127       lp->setColUpper(numcols, cutoff);
00128       subMip->find_good_sol(DBL_MAX, parameters_.subMilpLogLevel_,
00129       //subMip->optimize(DBL_MAX, parameters_.subMilpLogLevel_,
00130           (parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime()) /* time limit */);
00131 
00132       milpOptimal = subMip -> optimal(); 
00133       colsol = subMip->getLastSolution();
00134       nLocalSearch_++;
00135       if(milpOptimal)
00136         handler_->message(SOLVED_LOCAL_SEARCH, messages_)
00137         <<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
00138       else
00139         handler_->message(LOCAL_SEARCH_ABORT, messages_)
00140         <<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
00141     }
00142     int numberPasses = 0;
00143 
00144 #ifdef OA_DEBUG
00145     bool foundSolution = 0;
00146 #endif
00147     double * nlpSol = NULL;
00148     int major_iteration = 0;
00149     double ub = cutoff;
00150     while (colsol) {
00151       numberPasses++;
00152 
00153       //after a prescribed elapsed time give some branch_information to user
00154       //double time = CoinCpuTime();
00155 
00156 
00157       //setup the nlp
00158       int numberCutsBefore = cs.sizeRowCuts();
00159 
00160       //Fix the variable which have to be fixed, after having saved the bounds
00161       branch_info.solution_ = colsol;
00162 
00163       vector<double> x_bar(indices.size());
00164       for(unsigned int i = 0 ; i < indices.size() ; i++){
00165          assert(fabs(colsol[indices[i]] - floor(colsol[indices[i]] + 0.5)) < 1e-5);
00166          x_bar[i] = colsol[indices[i]];
00167       }
00168 
00169       double dist = nlp_->solveFeasibilityProblem(indices.size(), x_bar(), indices(), 1, 0, 2);
00170 
00171       handler_->message(FP_DISTANCE, messages_) 
00172       <<dist<<CoinMessageEol;
00173 
00174       if(dist < 1e-05){
00175          fixIntegers(*nlp_,branch_info, parameters_.cbcIntegerTolerance_, objects_, nObjects_);
00176 
00177          nlp_->resolve(txt_id);
00178          if(!nlp_->isProvenOptimal()){
00179            relaxIntegers(*nlp_,branch_info, parameters_.cbcIntegerTolerance_, objects_, nObjects_);
00180            nlp_->resolve(txt_id);
00181          }
00182          bool restart = false;
00183          if (post_nlp_solve(babInfo, cutoff)) {
00184            restart = true;
00185            //nlp is solved and feasible
00186            // Update the cutoff
00187            ub = std::min(ub, nlp_->getObjValue());
00188            cutoff = ub * (1 - parameters_.cbcCutoffIncrement_);
00189            
00190            numSols_++;
00191          }
00192          else{
00193            //nlp_->setColLower(savedColLower());
00194            //nlp_->setColUpper(savedColUpper());
00195            //dist = nlp_->solveFeasibilityProblem(indices.size(), x_bar(), indices(), 1, 0, 2);
00196          }
00197          nlpSol = const_cast<double *>(nlp_->getColSolution());
00198          nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
00199                                   parameter().global_);
00200          //if(restart){
00201            nlp_->setColLower(savedColLower());
00202            nlp_->setColUpper(savedColUpper());
00203         if(restart){
00204           major_iteration++;
00205           handler_->message(FP_MAJOR_ITERATION, messages_) 
00206           <<major_iteration<<cutoff<<CoinMessageEol;
00207           nlp_->resolve(txt_id);
00208         }
00209 
00210          //}
00211       }
00212       else {
00213          nlpSol = const_cast<double *>(nlp_->getColSolution());
00214          nlp_->getOuterApproximation(cs, nlpSol, 1, NULL,
00215                                   parameter().global_);
00216       }
00217 
00218 
00219 #if 0
00220       handler_->message(FP_MINOR_ITERATION, messages_) 
00221       <<nLocalSearch_<<cutoff<<CoinMessageEol;
00222 #endif
00223       
00224       int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
00225       assert(numberCuts);
00226       installCuts(*lp, cs, numberCuts);
00227       numberCutsBefore = cs.sizeRowCuts();
00228      
00229       //check time
00230       if (CoinCpuTime() - timeBegin_ > parameters_.maxLocalSearchTime_)
00231         break;
00232       //do we perform a new local search ?
00233       if (nLocalSearch_ < parameters_.maxLocalSearch_ &&
00234           numSols_ < parameters_.maxSols_) {
00235 
00237         if (subMip == NULL) subMip = new SubMipSolver(lp, parameters_.strategy());
00238 
00239         nLocalSearch_++;
00240         set_fp_objective(*lp, nlp_->getColSolution());
00241 
00242         lp->setColUpper(numcols, cutoff); 
00243      
00244         subMip->find_good_sol(DBL_MAX, parameters_.subMilpLogLevel_,
00245         // subMip->optimize(DBL_MAX, parameters_.subMilpLogLevel_,
00246                          parameters_.maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
00247         milpOptimal = subMip -> optimal(); 
00248         colsol = subMip->getLastSolution();
00249       if(milpOptimal)
00250         handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
00251       else
00252         handler_->message(LOCAL_SEARCH_ABORT, messages_)<<subMip->nodeCount()<<subMip->iterationCount()<<CoinMessageEol;
00253       if(colsol)
00254         handler_->message(FP_MILP_VAL, messages_) 
00255         <<colsol[nlp_->getNumCols()]<<CoinMessageEol;
00256          
00257       }
00258       else if (subMip!=NULL) {
00259         delete subMip;
00260         subMip = NULL;
00261         colsol = NULL;
00262         milpOptimal = false;
00263       }
00264     }
00265     if(colsol || ! milpOptimal)
00266       return -DBL_MAX;
00267     else{
00268       handler_->message(OASUCCESS, messages_)<<"FP"<<CoinCpuTime() - timeBegin_ 
00269       <<ub<<CoinMessageEol;
00270       return DBL_MAX;
00271     }
00272   }
00273 
00275   void
00276   MinlpFeasPump::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00277   {
00278     roptions->SetRegisteringCategory("Options for feasibility pump", RegisteredOptions::BonminCategory);
00279 
00280     roptions->AddBoundedIntegerOption("fp_log_level",
00281         "specify FP iterations log level.",
00282         0,2,1,
00283         "Set the level of output of OA decomposition solver : "
00284         "0 - none, 1 - normal, 2 - verbose"
00285                                      );
00286     roptions->setOptionExtraInfo("fp_log_level",3);
00287 
00288     roptions->AddLowerBoundedNumberOption("fp_log_frequency",
00289         "display an update on lower and upper bounds in FP every n seconds",
00290         0.,1.,100.,
00291         "");
00292     roptions->setOptionExtraInfo("fp_log_frequency",3);
00293   }
00294 
00296 void
00297 MinlpFeasPump::set_fp_objective(OsiSolverInterface &si, const double * colsol) const{
00298   if (objects_) {
00299     for (int i = 0 ; i < nObjects_ ; i++) {
00300       OsiObject * obj = objects_[i];
00301       int colnum = obj->columnNumber();
00302       if (colnum >= 0) {//Variable branching object
00303         double round = floor(colsol[colnum] + 0.5);
00304         double coeff = (colsol[colnum] - round ) < 0;
00305         si.setObjCoeff(colnum, 1 - 2 * coeff);
00306       }
00307       else {
00308         throw CoinError("OaDecompositionBase::solverManip",
00309                         "setFpObjective",
00310                         "Can not use FP on problem with SOS constraints");
00311       }
00312     }
00313   }
00314   else {
00315     int numcols = nlp_->getNumCols();
00316     for (int i = 0; i < numcols ; i++) {
00317       if (nlp_->isInteger(i)){
00318          double round = floor(colsol[i] + 0.5);
00319          double coeff = (colsol[i] - round ) < 0;
00320          si.setObjCoeff(i, 1 - 2*coeff);
00321       }
00322     }
00323   }
00324   si.initialSolve();
00325 }
00326 
00327 }/* End namespace Bonmin. */

Generated on Tue Mar 30 03:04:33 2010 by  doxygen 1.4.7