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

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines (IBM) 2006, 2007
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // P. Bonami, International Business Machines
00007 //
00008 // Date :  12/20/2006
00009 //#define ECP_DEBUG
00010 #include "BonEcpCuts.hpp"
00011 #include "BonSolverHelp.hpp"
00012 #include "BonBabInfos.hpp"
00013 #include "BonCbc.hpp"
00014 namespace Bonmin
00015 {
00016 
00017 
00018   EcpCuts::EcpCuts(BabSetupBase & b):
00019       OaDecompositionBase(b,false, false)
00020   {
00021     assignLpInterface(NULL);
00022     b.options()->GetIntegerValue("ecp_max_rounds", numRounds_,b.prefix());
00023     b.options()->GetNumericValue("ecp_abs_tol", abs_violation_tol_,b.prefix());
00024     b.options()->GetNumericValue("ecp_rel_tol", rel_violation_tol_,b.prefix());
00025     b.options()->GetNumericValue("ecp_propability_factor", beta_,b.prefix());
00026   }
00027 
00028   double
00029   EcpCuts::doEcpRounds(OsiSolverInterface &si,
00030       bool leaveSiUnchanged,
00031       double* violation)
00032   {
00033     OsiSolverInterface * saveLp = lp_;
00034     lp_ = &si;
00035     OsiCuts cs;
00036     bool saveLeaveSi = leaveSiUnchanged_;
00037     leaveSiUnchanged_ = leaveSiUnchanged;
00038     generateCuts(si, cs);
00039     lp_ = saveLp;
00040     leaveSiUnchanged_ = saveLeaveSi;
00041     if (violation) *violation = violation_;
00042     return objValue_;
00043   }
00044 
00045   void
00046   EcpCuts::generateCuts(const OsiSolverInterface &si,
00047       OsiCuts & cs,
00048       const CglTreeInfo info) const
00049   {
00050     if (beta_ >=0) {
00051       BabInfo * babInfo = dynamic_cast<BabInfo *> (si.getAuxiliaryInfo());
00052       assert(babInfo);
00053       assert(babInfo->babPtr());
00054       const CbcNode * node = babInfo->babPtr()->model().currentNode();
00055       int level = (node == NULL) ? 0 : babInfo->babPtr()->model().currentNode()->depth();
00056       double rand = CoinDrand48();
00057       double score = pow(2.,-level)*beta_;
00058       //printf("depth %i, score %g , rand %g\n", level, score, rand);
00059       if (score <= rand)
00060         return;
00061     }
00062     double orig_violation = nlp_->getNonLinearitiesViolation(
00063           si.getColSolution(), si.getObjValue());
00064 //#define ECP_DEBUG
00065 #ifdef ECP_DEBUG
00066     std::cout<<"Initial Constraint violation: "<<orig_violation<<std::endl;
00067     std::cout<<"Initial objectvie value"<<si.getObjValue()<<std::endl;
00068 #endif
00069     if (orig_violation <= abs_violation_tol_)
00070       return;
00071 #ifdef ECP_DEBUG
00072     std::cout<<"Generating ECP cuts"<<std::endl;
00073 #endif
00074     solverManip * lpManip = NULL;
00075     bool infeasible = false;
00076     violation_ = orig_violation;
00077     for (int i = 0 ; i < numRounds_ ; i++) {
00078       if ( violation_ > abs_violation_tol_ &&
00079           violation_ > rel_violation_tol_*orig_violation) {
00080         int numberCuts =  - cs.sizeRowCuts();
00081         const double * toCut = parameter().addOnlyViolated_ ?
00082             si.getColSolution():NULL;
00083         const OsiSolverInterface &localSi = (lpManip == NULL) ?
00084             si : *(lpManip->si());
00085         nlp_->getOuterApproximation(cs, localSi.getColSolution(), 1, toCut, parameter().global_);
00086         numberCuts += cs.sizeRowCuts();
00087         if (numberCuts > 0 && i + 1 < numRounds_ ) {
00088           if (lpManip==NULL) {
00089             assert(lp_ == NULL);
00090             if (lp_ == NULL)
00091               lpManip = new solverManip(si);
00092             else
00093               lpManip = new solverManip(lp_, true,true,
00094                   false,false);
00095           }
00096           installCuts(*lpManip->si(), cs,numberCuts);
00097 #ifdef ECP_DEBUG
00098           std::cerr<<"Installed "<<numberCuts<<"cuts in lp"<<std::endl;
00099 #endif
00100           lpManip->si()->resolve();
00101 #ifdef ECP_DEBUG
00102           std::cerr<<"New objective "<<lpManip->si()->getObjValue()<<std::endl;
00103 #endif
00104           if (lpManip->si()->isProvenPrimalInfeasible()) {
00105             infeasible = true;
00106 #ifdef ECP_DEBUG
00107             std::cout<<"Stopping Ecp generation because problem became infeasible"<<std::endl;
00108 #endif
00109             break;
00110           }
00111           violation_ =  nlp_->getNonLinearitiesViolation(
00112                 lpManip->si()->getColSolution(),
00113                 lpManip->si()->getObjValue());
00114 #ifdef ECP_DEBUG
00115           std::cout<<"Constraint violation: "<<violation_<<std::endl;
00116 #endif
00117         }
00118         else break;
00119       }
00120       else break;
00121     }
00122     if (!infeasible) {
00123       if (lpManip != NULL) {
00124         lpManip->si()->resolve();
00125         if (lpManip->si()->isProvenPrimalInfeasible())
00126           objValue_ = COIN_DBL_MAX;
00127         else
00128           objValue_ = lpManip->si()->getObjValue();
00129       }
00130     }
00131     else objValue_ = COIN_DBL_MAX;
00132     if (lpManip) {
00133       if (lp_ != NULL && lpManip != NULL) {
00134         lpManip->restore();
00135       }
00136 
00137       delete lpManip;
00138     }
00139 #ifdef ECP_DEBUG
00140     std::cout<<"End ecp cut generation"<<std::endl;
00141 #endif
00142     return;
00143   }
00144 
00145   void
00146   EcpCuts::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00147   {
00148     roptions->SetRegisteringCategory("Options for ecp cuts generation", RegisteredOptions::BonminCategory);
00149     roptions->AddLowerBoundedIntegerOption("filmint_ecp_cuts",
00150         "Specify the frequency (in terms of nodes) at which some a la filmint ecp cuts are generated.",
00151         0,0,
00152         "A frequency of 0 amounts to to never solve the NLP relaxation.");
00153     roptions->setOptionExtraInfo("filmint_ecp_cuts",3);
00154     roptions->AddLowerBoundedIntegerOption
00155     ("ecp_max_rounds",
00156      "Set the maximal number of rounds of ECP cuts.",
00157      0,5,
00158      "");
00159     roptions->setOptionExtraInfo("ecp_max_rounds",3);
00160     roptions->AddLowerBoundedNumberOption
00161     ("ecp_abs_tol",
00162      "Set the absolute termination tolerance for ECP rounds.",
00163      0,false,1e-6,
00164      "");
00165     roptions->setOptionExtraInfo("ecp_abs_tol",3);
00166     roptions->AddLowerBoundedNumberOption
00167     ("ecp_rel_tol",
00168      "Set the relative termination tolerance for ECP rounds.",
00169      0,false,0.,
00170      "");
00171     roptions->setOptionExtraInfo("ecp_rel_tol",3);
00172     roptions->AddNumberOption
00173     ("ecp_propability_factor",
00174      "Factor appearing in formula for skipping ECP cuts.",
00175      10.,
00176      "Choosing -1 disables the skipping.");
00177     roptions->setOptionExtraInfo("ecp_propability_factor",3);
00178   }
00179 } // end namespace bonmin.

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