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