// (C) Copyright International Business Machines Corporation 2010 // and others. // All Rights Reserved. // This code is published under the Eclipse Public License (EPL). // // Authors : // Giacomo Nannicini, Tepper School of Business, Carnegie Mellon University // // Date : 07/25/2010 #include "CouenneIterativeRounding.hpp" #include "BonTMINLP2Quad.hpp" #include "BonTMINLPLinObj.hpp" #include "CouenneRecordBestSol.hpp" // if we use OsiCpx, then we need to have cplex.h, too #if defined(COUENNE_HAS_OSICPX) && !defined(COUENNE_HAS_CPLEX) #undef COUENNE_HAS_OSICPX #endif #ifdef COUENNE_HAS_OSICPX #include "OsiCpxSolverInterface.hpp" #include "cplex.h" #endif #define MILPTIME 5 #define CBCMILPTIME 20 namespace Couenne{ CouenneIterativeRounding::CouenneIterativeRounding(): CbcHeuristic(), nlp_(NULL), cinlp_(NULL), milp_(NULL), maxRoundingIter_(10), maxFirPoints_(5), maxTime_(60), maxTimeFirstCall_(60), numInitialRows_(0), numSol_(-1), colLower_(NULL), colUpper_(NULL), colLowerNlp_(NULL), colUpperNlp_(NULL), omega_(0.2), baseLbRhs_(15), couenne_(NULL) { setHeuristicName("CouenneIterativeRounding"); } CouenneIterativeRounding::CouenneIterativeRounding(Bonmin::OsiTMINLPInterface* nlp, OsiSolverInterface* cinlp, CouenneProblem* couenne, Ipopt::SmartPtr options): CbcHeuristic(), nlp_(NULL), cinlp_(NULL), milp_(NULL), numSol_(-1), colLower_(NULL), colUpper_(NULL), colLowerNlp_(NULL), colUpperNlp_(NULL), couenne_(couenne) { // Initialize Dynamic NLP setNlp(nlp, cinlp); // Read options int irAggressiveness; options->GetIntegerValue("iterative_rounding_aggressiveness", irAggressiveness, "couenne."); setAggressiveness(irAggressiveness); double maxTime, maxTimeInit; options->GetNumericValue("iterative_rounding_time", maxTime, "couenne."); options->GetNumericValue("iterative_rounding_time_firstcall", maxTimeInit, "couenne."); if (maxTime >= 0){ setMaxTime(maxTime); } if (maxTimeInit >= 0){ setMaxTimeFirstCall(maxTimeInit); } int irLbrhs; options->GetIntegerValue("iterative_rounding_base_lbrhs", irLbrhs, "couenne."); setBaseLbRhs(irLbrhs); int numFirPoints; options->GetIntegerValue("iterative_rounding_num_fir_points", numFirPoints, "couenne."); setMaxFirPoints(numFirPoints); double omega; options->GetNumericValue("iterative_rounding_omega", omega, "couenne."); setOmega(omega); } CouenneIterativeRounding::CouenneIterativeRounding(const CouenneIterativeRounding & other): CbcHeuristic(other), nlp_(other.nlp_), cinlp_(other.cinlp_), milp_(other.milp_), maxRoundingIter_(other.maxRoundingIter_), maxFirPoints_(other.maxFirPoints_), maxTime_(other.maxTime_), maxTimeFirstCall_(other.maxTimeFirstCall_), numInitialRows_(other.numInitialRows_), numSol_(other.numSol_), omega_(other.omega_), baseLbRhs_(other.baseLbRhs_), couenne_(other.couenne_) { if(nlp_ != NULL){ nlp_ = dynamic_cast(other.nlp_->clone()); } if(milp_ != NULL) #ifdef COUENNE_HAS_OSICPX milp_ = dynamic_cast(other.milp_->clone()); #else milp_ = dynamic_cast(other.milp_->clone()); #endif if (other.colLower_ != NULL){ if (colLower_ != NULL) delete colLower_; colLower_ = new double[milp_->getNumCols()]; CoinCopyN (other.colLower_, milp_->getNumCols(), colLower_); } if (other.colUpper_ != NULL){ if (colUpper_ != NULL) delete colUpper_; colUpper_ = new double[milp_->getNumCols()]; CoinCopyN (other.colUpper_, milp_->getNumCols(), colUpper_); } if (other.colLowerNlp_ != NULL){ if (colLowerNlp_ != NULL) delete colLowerNlp_; colLowerNlp_ = new double[nlp_->getNumCols()]; CoinCopyN (other.colLowerNlp_, nlp_->getNumCols(), colLowerNlp_); } if (other.colUpperNlp_ != NULL){ if (colUpperNlp_ != NULL) delete colUpperNlp_; colUpperNlp_ = new double[nlp_->getNumCols()]; CoinCopyN (other.colUpperNlp_, nlp_->getNumCols(), colLowerNlp_); } } CbcHeuristic * CouenneIterativeRounding::clone() const{ return new CouenneIterativeRounding(*this); } CouenneIterativeRounding & CouenneIterativeRounding::operator=(const CouenneIterativeRounding & rhs){ if(this != &rhs){ CbcHeuristic::operator=(rhs); if(nlp_) delete nlp_; if(rhs.nlp_ != NULL){ nlp_ = dynamic_cast(rhs.nlp_->clone()); } cinlp_ = rhs.cinlp_; maxRoundingIter_ = rhs.maxRoundingIter_; maxFirPoints_ = rhs.maxFirPoints_; maxTime_ = rhs.maxTime_; maxTimeFirstCall_ = rhs.maxTimeFirstCall_; numSol_ = rhs.numSol_; numInitialRows_ = rhs.numInitialRows_; omega_ = rhs.omega_; baseLbRhs_ = rhs.baseLbRhs_; couenne_ = rhs.couenne_; if (rhs.colLower_ != NULL){ if (colLower_ != NULL) delete colLower_; colLower_ = new double[milp_->getNumCols()]; CoinCopyN (rhs.colLower_, milp_->getNumCols(), colLower_); } if (rhs.colUpper_ != NULL){ if (colUpper_ != NULL) delete colUpper_; colUpper_ = new double[milp_->getNumCols()]; CoinCopyN (rhs.colUpper_, milp_->getNumCols(), colLower_); } if (rhs.colLowerNlp_ != NULL){ if (colLowerNlp_ != NULL) delete colLowerNlp_; colLowerNlp_ = new double[nlp_->getNumCols()]; CoinCopyN (rhs.colLowerNlp_, nlp_->getNumCols(), colLowerNlp_); } if (rhs.colUpperNlp_ != NULL){ if (colUpperNlp_ != NULL) delete colUpperNlp_; colUpperNlp_ = new double[nlp_->getNumCols()]; CoinCopyN (rhs.colUpperNlp_, nlp_->getNumCols(), colLowerNlp_); } } return *this; } CouenneIterativeRounding::~CouenneIterativeRounding(){ delete nlp_; nlp_ = NULL; if (colLower_) delete[] colLower_; if (colUpper_) delete[] colUpper_; if (colLowerNlp_) delete[] colLowerNlp_; if (colUpperNlp_) delete[] colUpperNlp_; if (milp_) delete milp_; milp_ = NULL; } void CouenneIterativeRounding::setNlp(Bonmin::OsiTMINLPInterface* nlp, OsiSolverInterface * cinlp){ // Create a dynamic NLP (i.e. one that can be used to add/remove // linear inequalities) from the initial one if(nlp_ != NULL) delete nlp_; nlp_ = new Bonmin::OsiTMINLPInterface; Ipopt::SmartPtr tminlp = nlp->model(); if (tminlp->hasLinearObjective()){ Ipopt::SmartPtr linObj = new Bonmin::TMINLPLinObj; linObj->setTminlp(GetRawPtr(tminlp)); tminlp = GetRawPtr(linObj); } nlp_->initialize(nlp->regOptions(), nlp->options(), nlp->solver()->journalist(), "bonmin.", tminlp); nlp_->use(new Bonmin::TMINLP2TNLPQuadCuts(tminlp)); cinlp_ = cinlp; } void CouenneIterativeRounding::setMilp(){ if(milp_ != NULL) delete milp_; // Save the LP relaxation of the root node OsiSolverInterface * milp = model_->solver(); int n = milp->getNumCols(); #ifdef COUENNE_HAS_OSICPX milp_ = new OsiCpxSolverInterface(); milp_->loadProblem(*(milp->getMatrixByRow()), milp->getColLower(), milp->getColUpper(), milp->getObjCoefficients(), milp->getRowLower(), milp->getRowUpper()); for (int i = 0; i < n; ++i){ if (milp->isInteger(i)) milp_->setInteger(i); } #else milp_ = dynamic_cast(milp->clone()); #endif colLower_ = new double[n]; colUpper_ = new double[n]; memcpy(colLower_, milp->getColLower(), n*sizeof(double)); memcpy(colUpper_, milp->getColUpper(), n*sizeof(double)); int nNlp = cinlp_->getNumCols(); colLowerNlp_ = new double[nNlp]; colUpperNlp_ = new double[nNlp]; memcpy(colLowerNlp_, cinlp_->getColLower(), nNlp*sizeof(double)); memcpy(colUpperNlp_, cinlp_->getColUpper(), nNlp*sizeof(double)); numIntegers_ = 0; for (int i = 0; i < nNlp; ++i){ if (cinlp_->isInteger(i)){ numIntegers_++; } } // workaround for inconsistent bounds: // sometimes using getColLower and getColUpper we get *slightly* // inconsistent bounds, which make everything crash afterwards double swap; for (int i = 0; i < n; ++i){ if (colUpper_[i] < colLower_[i]){ swap = colUpper_[i]; colUpper_[i] = colLower_[i]; colLower_[i] = swap; } } numInitialRows_ = milp_->getNumRows(); // this array contains zeroes; we use it to set empty elements double * tmpArray = new double[n]; CoinFillN(tmpArray, n, 0.0); milp_->setObjective(tmpArray); milp_->setObjSense(1); // now there is no objective function // create new variables with objective coefficient 1 for (int i = 0; i < n; ++i){ milp_->addCol(0, NULL, NULL, 0.0, COIN_DBL_MAX, 1.0); } milp_->setHintParam(OsiDoDualInResolve,true,OsiHintDo); milp_->setHintParam(OsiDoPresolveInResolve,true,OsiHintDo); milp_->setHintParam(OsiDoReducePrint,true,OsiHintDo); milp_->setDblParam(OsiPrimalTolerance, COUENNE_EPS_INT); milp_->messageHandler()->setLogLevel(0); milp_->setDblParam(OsiDualTolerance, COUENNE_EPS_INT); #ifdef COUENNE_HAS_OSICPX // Set options if we have Cplex CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_MIPEMPHASIS, CPX_MIPEMPHASIS_HIDDENFEAS); CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_FPHEUR, 2); CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_NODESEL, 0); CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_LBHEUR, 1); CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_RINSHEUR, 99); CPXsetintparam(milp_->getEnvironmentPtr(), CPX_PARAM_NODELIM, 50); CPXsetdblparam(milp_->getEnvironmentPtr(), CPX_PARAM_CUTSFACTOR, 1.0); CPXsetdblparam(milp_->getEnvironmentPtr(), CPX_PARAM_TILIM, MILPTIME); #else // Set up heuristics if we use Cbc (not recommended) heuristics_ = new CbcHeuristic*[1]; numHeuristics_ = 1; CbcHeuristicFPump * feaspump = new CbcHeuristicFPump(); feaspump->setMaximumRetries(2); feaspump->setMaximumPasses(100); feaspump->setAccumulate(3); heuristics_[0] = feaspump; #endif delete[] tmpArray; } void CouenneIterativeRounding::setAggressiveness(int value){ switch (value){ case 0: setMaxRoundingIter(5); setMaxTimeFirstCall(300); setMaxFirPoints(5); setMaxTime(60); break; case 1: setMaxRoundingIter(10); setMaxTimeFirstCall(300); setMaxFirPoints(5); setMaxTime(120); break; case 2: setMaxRoundingIter(20); setMaxTimeFirstCall(1000); setMaxFirPoints(5); setMaxTime(300); break; default: std::cerr << "CouenneIterativeRounding::setAggressiveness() : unknown value!\n" << std::endl; } } int CouenneIterativeRounding::solution(double & objectiveValue, double* newSolution){ if (milp_ == NULL){ // On the very first call, we set the MILP setMilp(); return 0; } if ((model_->numberIntegers() == 0) || (numSol_ == model_->getSolutionCount())){ // There are no integers, or we already tried to improve current // solution. Exit immediately. return 0; } numSol_ = model_->getSolutionCount(); std::cout << "Launching IterativeRounding with parameters:" << std::endl; std::cout << "Max rounding iter: " << maxRoundingIter_ << std::endl; std::cout << "Max feas point: " << maxFirPoints_ << std::endl; std::cout << "Base lbrhs: " << baseLbRhs_ << std::endl; std::cout << "Omega: " << omega_ << std::endl; std::cout << "Max time firstcall: " << maxTimeFirstCall_ << std::endl; // write down starting time and apply heuristics startTime_ = CoinCpuTime(); endTime_ = ((numSol_ == 0) ? maxTimeFirstCall_ : maxTime_); const double* bestKnownSolution = model_->bestSolution(); bool found = false; bool hasSolution = true; bool improved = true; if (numSol_ == 0){ // No solution available; start with our feasibility heuristic hasSolution = feasibilityIR(objectiveValue, newSolution); if (hasSolution){ bestKnownSolution = newSolution; found = hasSolution; } } if (!hasSolution){ // Still do not have a starting solution, we cannot run improvementIR return found; } while (improved && (CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)){ // Keep improving current solution as long as it works improved = false; improved = improvementIR(objectiveValue, newSolution, bestKnownSolution); if (improved){ bestKnownSolution = newSolution; } found = (found || improved); } if (found){ // make sure we do not run improvementIR again on the same solution numSol_++; } return found; } int CouenneIterativeRounding::feasibilityIR(double& objectiveValue, double* newSolution){ std::cout << "starting feasibility IR" << std::endl; OsiSolverInterface * solver = model_->solver(); OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo(); Bonmin::BabInfo * babInfo = dynamic_cast (auxInfo); if(babInfo){ babInfo->setHasNlpSolution(false); if(babInfo->infeasibleNode()){ return 0; } } int n = couenne_->nVars(); int nNlp = cinlp_->getNumCols(); // get best solution and count the number of binary variables int numIntAtBound = 0; OsiRowCut cut; OsiRowCut lbcut1; OsiRowCut lbcut2; OsiCuts lbcuts; double obj; bool boundsChanged = false; std::vector previousBranches; // Try to find a feasible solution // apply NLP solver to obtain point xprime nlp_->resolve(); if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() || nlp_->isAbandoned()){ nlp_->resolveForRobustness(3); } if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() || nlp_->isAbandoned()){ obj = COIN_DBL_MAX; } else{ obj = nlp_->getObjValue(); } if (obj == COIN_DBL_MAX){std::cout << " IR: no feasible solution found " << std::endl; std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl; return 0; } double* xprime = new double [n]; CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime); couenne_ -> getAuxs (xprime); // now prepare the MILP that we will solve // the objective is : minimize sum_{j = 1}^n w_j int constrInd[2]; double constrElem[2]; // main loop: solve a sequence of MILPs and NLPs bool foundSolution = false; double* tmpSolution = new double[n]; OsiCuts revlb_all; OsiSolverInterface* curr_milp = milp_; int outerLoop = maxFirPoints_; for (int h = 0; h < outerLoop && ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++h){ OsiCuts cs; cs.insert(lbcut1); cs.insert(lbcut2); // write constraints: // xbar[j] - x[j] <= w_j // x[j] - xbar[j] <= w_j for (int i = 0; i < n; i++){ constrInd[0] = i; constrInd[1] = i+n; constrElem[0] = -1; constrElem[1] = -1; cut.mutableRow().setVector(2, constrInd, constrElem); cut.setLb(-COIN_DBL_MAX); cut.setUb(-xprime[i]); cs.insert(cut); constrElem[0] = +1; constrElem[1] = -1; cut.mutableRow().setVector(2, constrInd, constrElem); cut.setLb(-COIN_DBL_MAX); cut.setUb(xprime[i]); cs.insert(cut); } curr_milp->applyCuts(cs); for (int k = 0; k < maxRoundingIter_ && ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++k){ // Solve MILP to obtain an integral feasible point near xprime curr_milp->restoreBaseModel(numInitialRows_+cs.sizeRowCuts()); curr_milp->applyCuts(revlb_all); bool solFound = solveMilp(curr_milp, endTime_-(CoinCpuTime()-startTime_)-2); if (!solFound && !boundsChanged){ std::cout << " MILP cannot be solved, terminating LB " << std::endl; // we cannot solve the MILP and bounds were not changed; exit curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl; return foundSolution; } else if (!solFound && boundsChanged){ // the MILP is infeasible but bounds were changed; // restore original bounds and branch on a random variable // of the solution at the previous iteration (so we do not cycle) curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); branchToCut(tmpSolution, curr_milp, previousBranches); continue; } const double * xtilde = curr_milp->getColSolution(); // now fix integer variables and solve NLP // xtilde has more columns than the NLP, so this should be ok CoinCopyN (xtilde, n, tmpSolution); for (int i = 0; i < nNlp; ++i){ if (model_->isInteger(i)){ tmpSolution[i] = floor(tmpSolution[i]+0.5); cinlp_->setColLower(i, tmpSolution[i]); cinlp_->setColUpper(i, tmpSolution[i]); } } cinlp_->setColSolution(tmpSolution); cinlp_->messageHandler()->setLogLevel(1); cinlp_->resolve(); obj = ((cinlp_->isProvenOptimal()) ? cinlp_->getObjValue():COIN_DBL_MAX); memcpy(tmpSolution, cinlp_->getColSolution(), nNlp*sizeof(double)); // check solution of the NLP; // if we have a new incumbent we are done, otherwise we iterate bool isChecked = false; isChecked = couenne_ -> checkNLP0 (tmpSolution, obj, true, false, true, false); // #ifdef FM_CHECKNLP2 // isChecked = couenne_->checkNLP2(tmpSolution, 0, false, // do not care about obj // true, // stopAtFirstViol // false, // checkALL // couenne_->getFeasTol()); // if(isChecked) { // obj = couenne_->getRecordBestSol()->getModSolVal(); // } // #else /* not FM_CHECKNLP2 */ // isChecked = couenne_->checkNLP(tmpSolution, obj, true); // #endif /* not FM_CHECKNLP2 */ if (cinlp_->isProvenOptimal () && isChecked && (obj < couenne_->getCutOff())) { #ifdef FM_CHECKNLP2 #ifdef FM_TRACE_OPTSOL couenne_->getRecordBestSol()->update(); CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution); obj = couenne_->getRecordBestSol()->getVal(); #else /* not FM_TRACE_OPTSOL */ CoinCopyN (couenne_->getRecordBestSol()->getModSol(n), n, tmpSolution); #endif /* not FM_TRACE_OPTSOL */ #else /* not FM_CHECKNLP2 */ //Get correct values for all auxiliary variables couenne_ -> getAuxs (tmpSolution); #ifdef FM_TRACE_OPTSOL couenne_->getRecordBestSol()->update(tmpSolution, n, obj, couenne_->getFeasTol()); CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution); obj = couenne_->getRecordBestSol()->getVal(); #endif /* FM_TRACE_OPTSOL */ #endif /* not FM_CHECKNLP2 */ if (babInfo){ babInfo->setNlpSolution (tmpSolution, n, obj); babInfo->setHasNlpSolution (true); } std::cout << "Final Nlp solution with objective " << obj << " :" << std::endl; if (obj < objectiveValue - COUENNE_EPS) { // found better solution? std::cout << "New incumbent found" << std::endl; const CouNumber *lb = solver -> getColLower (), *ub = solver -> getColUpper (); // check bounds once more after getAux. This avoids false // asserts in CbcModel.cpp:8305 on integerTolerance violated for (int i=0; i < n; ++i, ++lb, ++ub) { CouNumber &t = tmpSolution [i]; if (t < *lb) t = *lb; else if (t > *ub) t = *ub; } couenne_ -> setCutOff (obj); foundSolution = true; objectiveValue = obj; CoinCopyN (tmpSolution, n, newSolution); cinlp_->setColLower(colLowerNlp_); cinlp_->setColUpper(colUpperNlp_); curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } //delete curr_milp; double elapsed = CoinCpuTime()-startTime_; std::cout << "IR: Heuristic: " << objectiveValue << std::endl; std::cout << "IR: Heuristic time: " << elapsed << std::endl; return foundSolution; } } cinlp_->setColLower(colLowerNlp_); cinlp_->setColUpper(colUpperNlp_); double avgBound; numIntAtBound = computeIntAtBound(xtilde, avgBound); if (numIntAtBound >= 50 || numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){ // write reverse local branching constraint // this avoids finding the same xtilde again avgBound = floor(avgBound + 0.5); std::cout << "Using reverse LB with rhs " << avgBound << std::endl; writeLB(cut, xtilde, 'G', avgBound); revlb_all.insert(cut); } else{ // Branch on a random variable to obtain a different integer point branchToCut(xtilde, curr_milp, previousBranches); boundsChanged = true; } } if (h <= outerLoop -2){ // Compute a NLP-feasible point by solving a log-barrier problem // with a given minimum value for the log-barrier parameter mu Bonmin::OsiTMINLPInterface * nlp = dynamic_cast (nlp_); Ipopt::SmartPtr< Ipopt::OptionsList > opt = nlp->solver()->options(); nlp->messageHandler()->setLogLevel(10); double mu_target; double kappa_d; opt->GetNumericValue("mu_target", mu_target, "ipopt."); opt->SetNumericValue("mu_target", omega_*(h+1), "ipopt."); opt->GetNumericValue("kappa_d", kappa_d, "ipopt."); opt->SetNumericValue("kappa_d", 0.0, "ipopt."); nlp_->resolve(); if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() || nlp_->isAbandoned()){ nlp_->resolveForRobustness(3); } opt->SetNumericValue("mu_target", mu_target, "ipopt."); opt->SetNumericValue("kappa_d", kappa_d, "ipopt."); if (!nlp->isProvenPrimalInfeasible() && !nlp->isProvenDualInfeasible() && !nlp->isAbandoned()){ CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime); couenne_ -> getAuxs (xprime); curr_milp->restoreBaseModel(numInitialRows_); if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } } } } curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } double elapsed = CoinCpuTime()-startTime_; std::cout << "IR: Heuristic: " << COUENNE_INFINITY << std::endl; std::cout << "IR: Heuristic time: " << elapsed << std::endl; return foundSolution; } int CouenneIterativeRounding::improvementIR(double& objectiveValue, double* newSolution, const double* solution){ std::cout << "starting Improvement IR" << std::endl; OsiSolverInterface * solver = model_->solver(); OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo(); Bonmin::BabInfo * babInfo = dynamic_cast (auxInfo); if(babInfo){ babInfo->setHasNlpSolution(false); if(babInfo->infeasibleNode()){ return 0; } } int n = couenne_->nVars(); int nNlp = cinlp_->getNumCols(); // count the number of binary variables int numIntAtBound = 0; double lbrhs = CoinMin(baseLbRhs_, CoinMax(1,numIntegers_/2) ); OsiRowCut cut; OsiRowCut lbcut1; OsiRowCut lbcut2; OsiCuts lbcuts; int currentIndex = 0; double obj; bool boundsChanged = false; std::vector previousBranches; double avgBound = 0.0; numIntAtBound = computeIntAtBound(solution, avgBound); if (numIntAtBound >= 50 || numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){ // write local branching constraint writeLB(lbcut1, solution, 'L', lbrhs + floor(avgBound - 0.5)); lbcuts.insert(lbcut1); writeLB(lbcut2, solution, 'G', 1); lbcuts.insert(lbcut2); // We cannot add a row because OsiTMINLP will not let us do that! // Therefore we add a cut -- does the same thing. // We must delete the cut afterwards! nlp_->applyCuts(lbcuts); } else{ // we cannot use a local branching constraint, so we simply // reduce the size of the box bounds, to have a smaller NLP and // stay somehow "close" to the current point for (int i = 0; i < nlp_->getNumCols(); ++i){ if (nlp_->isInteger(i)){ nlp_->setColLower(i, colLowerNlp_[i]+(solution[i]-colLower_[i])*0.5); nlp_->setColUpper(i, colUpperNlp_[i]+(solution[i]-colUpper_[i])*0.5); } } } // apply NLP solver to obtain point xprime nlp_->setColSolution(solution); nlp_->resolve(); if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() || nlp_->isAbandoned()){ nlp_->resolveForRobustness(3); } if (nlp_->isProvenPrimalInfeasible() || nlp_->isProvenDualInfeasible() || nlp_->isAbandoned()){ obj = COIN_DBL_MAX; } else{ obj = nlp_->getObjValue(); } // Restore NLP if (numIntAtBound > 0){ currentIndex = nlp_->getNumRows()-1; nlp_->deleteRows(1, ¤tIndex); currentIndex = nlp_->getNumRows()-1; nlp_->deleteRows(1, ¤tIndex); } else{ nlp_->setColLower(colLowerNlp_); nlp_->setColUpper(colUpperNlp_); } if (obj == COIN_DBL_MAX || obj >= objectiveValue - COUENNE_EPS){ std::cout << " IR: no improvement possible " << std::endl; std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl; return 0; } double* xprime = new double [n]; CoinCopyN (nlp_->getColSolution(), nlp_->getNumCols(), xprime); couenne_ -> getAuxs (xprime); // now prepare the MILP that we will solve // the objective is : minimize sum_{j = 1}^n w_j int constrInd[2]; double constrElem[2]; // main loop: solve a sequence of MILPs and NLPs bool foundSolution = false; double* tmpSolution = new double[n]; OsiCuts revlb_all; OsiSolverInterface* curr_milp = milp_; OsiCuts cs; cs.insert(lbcut1); cs.insert(lbcut2); // write constraints: // xbar[j] - x[j] <= w_j // x[j] - xbar[j] <= w_j for (int i = 0; i < n; i++){ constrInd[0] = i; constrInd[1] = i+n; constrElem[0] = -1; constrElem[1] = -1; cut.mutableRow().setVector(2, constrInd, constrElem); cut.setLb(-COIN_DBL_MAX); cut.setUb(-xprime[i]); cs.insert(cut); constrElem[0] = +1; constrElem[1] = -1; cut.mutableRow().setVector(2, constrInd, constrElem); cut.setLb(-COIN_DBL_MAX); cut.setUb(xprime[i]); cs.insert(cut); } curr_milp->applyCuts(cs); for (int k = 0; k < maxRoundingIter_ && ((CoinCpuTime()-startTime_) < (endTime_ - MILPTIME)); ++k){ // Solve MILP to obtain an integral feasible point near xprime curr_milp->restoreBaseModel(numInitialRows_+cs.sizeRowCuts()); curr_milp->applyCuts(revlb_all); bool solFound = solveMilp(curr_milp, endTime_-(CoinCpuTime()-startTime_)-2); if (!solFound && !boundsChanged){ std::cout << " MILP cannot be solved, terminating LB " << std::endl; // we cannot solve the MILP and bounds were not changed; exit curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } std::cout << " IR: elapsed time " << CoinCpuTime()-startTime_ << std::endl; return foundSolution; } else if (!solFound && boundsChanged){ // the MILP is infeasible but bounds were changed; // restore original bounds and branch on a random variable // of the solution at the previous iteration (so we do not cycle) curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); branchToCut(tmpSolution, curr_milp, previousBranches); continue; } const double * xtilde = curr_milp->getColSolution(); // now fix integer variables and solve NLP // xtilde has more columns than the NLP, so this should be ok CoinCopyN (xtilde, n, tmpSolution); for (int i = 0; i < nNlp; ++i){ if (model_->isInteger(i)){ tmpSolution[i] = floor(tmpSolution[i]+0.5); cinlp_->setColLower(i, tmpSolution[i]); cinlp_->setColUpper(i, tmpSolution[i]); } } cinlp_->setColSolution(tmpSolution); cinlp_->messageHandler()->setLogLevel(1); cinlp_->resolve(); obj = ((cinlp_->isProvenOptimal()) ? cinlp_->getObjValue():COIN_DBL_MAX); memcpy(tmpSolution, cinlp_->getColSolution(), nNlp*sizeof(double)); // check solution of the NLP; // if we have a new incumbent we are done, otherwise we iterate bool isChecked = false; isChecked = couenne_ -> checkNLP0 (tmpSolution, obj, true, false, true, false); // #ifdef FM_CHECKNLP2 // isChecked = couenne_->checkNLP2(tmpSolution, 0, false, // do not care about obj // true, // stopAtFirstViol // false, // checkALL // couenne_->getFeasTol()); // if(isChecked) { // obj = couenne_->getRecordBestSol()->getModSolVal(); // } // #else /* not FM_CHECKNLP2 */ // isChecked = couenne_->checkNLP(tmpSolution, obj, true); // #endif /* not FM_CHECKNLP2 */ if (cinlp_->isProvenOptimal () && isChecked && (obj < couenne_->getCutOff())) { #ifdef FM_CHECKNLP2 #ifdef FM_TRACE_OPTSOL couenne_->getRecordBestSol()->update(); CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution); obj = couenne_->getRecordBestSol()->getVal(); #else /* not FM_TRACE_OPTSOL */ CoinCopyN (couenne_->getRecordBestSol()->getModSol(n), n, tmpSolution); #endif /* not FM_TRACE_OPTSOL */ #else /* not FM_CHECKNLP2 */ //Get correct values for all auxiliary variables couenne_ -> getAuxs (tmpSolution); #ifdef FM_TRACE_OPTSOL couenne_->getRecordBestSol()->update(tmpSolution, n, obj, couenne_->getFeasTol()); CoinCopyN (couenne_->getRecordBestSol()->getSol(), n, tmpSolution); obj = couenne_->getRecordBestSol()->getVal(); #endif /* FM_TRACE_OPTSOL */ #endif /* not FM_CHECKNLP2 */ if (babInfo){ babInfo->setNlpSolution (tmpSolution, n, obj); babInfo->setHasNlpSolution (true); } std::cout << "Final Nlp solution with objective " << obj << " :" << std::endl; if (obj < objectiveValue - COUENNE_EPS) { // found better solution? std::cout << "New incumbent found" << std::endl; const CouNumber *lb = solver -> getColLower (), *ub = solver -> getColUpper (); // check bounds once more after getAux. This avoids false // asserts in CbcModel.cpp:8305 on integerTolerance violated for (int i=0; i < n; ++i, ++lb, ++ub) { CouNumber &t = tmpSolution [i]; if (t < *lb) t = *lb; else if (t > *ub) t = *ub; } couenne_ -> setCutOff (obj); foundSolution = true; objectiveValue = obj; CoinCopyN (tmpSolution, n, newSolution); cinlp_->setColLower(colLowerNlp_); cinlp_->setColUpper(colUpperNlp_); curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } //delete curr_milp; double elapsed = CoinCpuTime()-startTime_; std::cout << "IR: Heuristic: " << objectiveValue << std::endl; std::cout << "IR: Heuristic time: " << elapsed << std::endl; return foundSolution; } } cinlp_->setColLower(colLowerNlp_); cinlp_->setColUpper(colUpperNlp_); numIntAtBound = computeIntAtBound(xtilde, avgBound); if (numIntAtBound >= 50 || numIntAtBound >= ((numIntegers_*0.1 > 5) ? numIntegers_*0.1 : 5)){ // write reverse local branching constraint // this avoids finding the same xtilde again avgBound = floor(avgBound + 0.5); std::cout << "Using reverse LB with rhs " << avgBound << std::endl; writeLB(cut, xtilde, 'G', avgBound); revlb_all.insert(cut); } else{ // Branch (randomly) so that we can obtain a new integer point branchToCut(xtilde, curr_milp, previousBranches); boundsChanged = true; } } curr_milp->restoreBaseModel(numInitialRows_); delete[] xprime; delete[] tmpSolution; if (boundsChanged){ curr_milp->setColLower(colLower_); curr_milp->setColUpper(colUpper_); } double elapsed = CoinCpuTime()-startTime_; std::cout << "IR: Heuristic: " << COUENNE_INFINITY << std::endl; std::cout << "IR: Heuristic time: " << elapsed << std::endl; return foundSolution; } int CouenneIterativeRounding::computeIntAtBound(const double* x){ int numIntAtBound = 0; for (int i = 0; i < nlp_->getNumCols(); ++i){ if (nlp_->isInteger(i) && (areEqual(x[i], colLower_[i]) || areEqual(x[i], colUpper_[i]))){ numIntAtBound++; } } return numIntAtBound; } int CouenneIterativeRounding::computeIntAtBound(const double* x, double& avgBoundSize){ int numIntAtBound = 0; avgBoundSize = 0; for (int i = 0; i < nlp_->getNumCols(); ++i){ if (nlp_->isInteger(i) && (areEqual(x[i], colLower_[i]) || areEqual(x[i], colUpper_[i]))){ numIntAtBound++; avgBoundSize += colUpper_[i] - colLower_[i]; } } avgBoundSize /= numIntAtBound; return numIntAtBound; } void CouenneIterativeRounding::writeLB(OsiRowCut& cut, const double* x, char sense, double rhs){ cut.mutableRow().clear(); for (int i = 0; i < nlp_->getNumCols(); ++i){ if (nlp_->isInteger(i)){ if (areEqual(x[i], colUpper_[i])){ cut.mutableRow().insert(i, -1); rhs -= x[i]; } else if (areEqual(x[i], colLower_[i])){ cut.mutableRow().insert(i, 1); rhs += x[i]; } } } if (sense == 'L'){ cut.setLb(-COIN_DBL_MAX); cut.setUb(rhs); } else if (sense == 'G'){ cut.setLb(rhs); cut.setUb(COIN_DBL_MAX); } else{ std::cerr << "### ERROR: wrong sense of LB constraint" << std::endl; exit(1); } } bool CouenneIterativeRounding::solveMilp(OsiSolverInterface* milp, double maxTime){ double start = CoinCpuTime(); #ifdef COUENNE_HAS_OSICPX OsiCpxSolverInterface * solver = dynamic_cast(milp); CPXENVptr env = solver->getEnvironmentPtr(); CPXLPptr cpxlp = solver->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL); int status = 0; bool solFound = false; bool infeasible = false; while (!solFound && !infeasible && ((CoinCpuTime() - start) < maxTime)){ solver->branchAndBound(); status = CPXgetstat(env, cpxlp); solFound = ((status == CPXMIP_NODE_LIM_FEAS) || (status == CPXMIP_TIME_LIM_FEAS) || (status == CPXMIP_MEM_LIM_FEAS) || (status == CPXMIP_ABORT_FEAS) || (status == CPXMIP_OPTIMAL) || (status == CPXMIP_OPTIMAL_TOL) || (status == CPXMIP_ABORT_FEAS) || (status == CPXMIP_FAIL_FEAS) || (status == CPXMIP_FAIL_FEAS_NO_TREE) || (status == CPXMIP_FEASIBLE)); infeasible = ((status == CPXMIP_INForUNBD) || (status == CPXMIP_OPTIMAL_INFEAS) || (status == CPXMIP_INFEASIBLE)); } if (solFound){ return true; } return false; #else CbcModel cbcModel(*milp); for (int i = 0; i < numHeuristics_; ++i){ cbcModel.addHeuristic(heuristics_[i]); } cbcModel.setMaximumSeconds(CBCMILPTIME); //CglPreProcess * prep = cbcModel.preProcess(0,2,5); while ((cbcModel.getSolutionCount() == 0) && (!cbcModel.isProvenInfeasible()) && (!cbcModel.isProvenDualInfeasible()) && (!cbcModel.isAbandoned()) && ((CoinCpuTime() - start) < maxTime)){ cbcModel.branchAndBound(); } //cbcModel.postProcess(prep); milp = cbcModel.solver(); if (cbcModel.getSolutionCount() > 0){ return true; } return false; #endif } int CouenneIterativeRounding::branchToCut(const double* x, OsiSolverInterface* solver, std::vector& previousBranches){ int branch; bool found = false; while (!found){ branch = rand()%numIntegers_; found = true; for (unsigned int i = 0; i < previousBranches.size(); ++i){ if (branch == previousBranches[i]){ found = false; break; } } if (found){ previousBranches.push_back(branch); } else{ continue; } for (int i = 0; i < nlp_->getNumCols(); ++i){ if (model_->isInteger(i)){ if (branch == 0){ branch = i; break; } else{ branch--; } } } } double sample = rand(); if (sample <= ((x[branch]-colLower_[branch])/(colUpper_[branch]-colLower_[branch]))){ solver->setColUpper(branch, x[branch]-1); } else{ solver->setColLower(branch, x[branch]+1); } return branch; } void CouenneIterativeRounding::registerOptions (Ipopt::SmartPtr roptions){ roptions -> AddStringOption2 ("iterative_rounding_heuristic", "Do we use the Iterative Rounding heuristic", "no", "no","", "yes","", "If enabled, a heuristic based on Iterative Rounding is used " "to find feasible solutions for the problem. " "The heuristic may take some time, but usually finds good solutions. " "Recommended if you want good upper bounds and have Cplex. " "Not recommended if you do not have Cplex"); roptions -> AddNumberOption ("iterative_rounding_time", "Specify the maximum time allowed for the Iterative Rounding heuristic", -1, "Maximum CPU time employed by the Iterative Rounding heuristic; " "if no solution found in this time, failure is reported. " "This overrides the CPU time set by Aggressiveness if positive."); roptions -> AddNumberOption ("iterative_rounding_time_firstcall", "Specify the maximum time allowed for the Iterative Rounding heuristic " "when no feasible solution is known", -1, "Maximum CPU time employed by the Iterative Rounding heuristic " "when no solution is known; if no solution found in this time, " "failure is reported." "This overrides the CPU time set by Aggressiveness if posive."); roptions -> AddBoundedIntegerOption ("iterative_rounding_aggressiveness", "Aggressiveness of the Iterative Rounding heuristic", 0, 2, 1, "Set the aggressiveness of the heuristic; i.e., how many iterations " "should be run, and with which parameters. The maximum time can be " "overridden by setting the _time and _time_firstcall options. " "0 = non aggressive, 1 = standard (default), 2 = aggressive."); roptions -> AddLowerBoundedIntegerOption ("iterative_rounding_num_fir_points", "Max number of points rounded at the beginning of Iterative Rounding", 1, 5, "Number of different points (obtained solving a log-barrier problem) " "that the heuristic will try to round at most, during its execution " "at the root node (i.e. the F-IR heuristic). Default 5."); roptions -> AddBoundedNumberOption ("iterative_rounding_omega", "Omega parameter of the Iterative Rounding heuristic", 0, true, 1, true, 0.2, "Set the omega parameter of the heuristic, which represents a " "multiplicative factor for the minimum log-barrier parameter " "of the NLP which is solved to obtain feasible points. This " "corresponds to $\\omega'$ in the paper. Default 0.2."); roptions -> AddLowerBoundedIntegerOption ("iterative_rounding_base_lbrhs", "Base rhs of the local branching constraint for Iterative Rounding", 0, 15, "Base rhs for the local branching constraint that defines a " "neighbourhood of the local incumbent. The base rhs is modified by " "the algorithm according to variable bounds. This corresponds to " "k' in the paper. Default 15."); } }/** Ends namespace Couenne.*/