00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include "BonminConfig.h"
00018
00019 #if defined(_MSC_VER)
00020
00021 # pragma warning(disable:4786)
00022 #endif
00023
00024 #include <cassert>
00025 #include <iomanip>
00026 #include <sstream>
00027
00028
00029
00030
00031 #include "OsiSolverInterface.hpp"
00032 #include "OsiClpSolverInterface.hpp"
00033 #include "CbcModel.hpp"
00034 #include "CbcBranchUser.hpp"
00035 #include "CbcCompareUser.hpp"
00036 #include "CbcCompareActual.hpp"
00037 #include "CbcCutGenerator.hpp"
00038
00039 #include "BonAmplInterface.hpp"
00040 #include "BonDummyHeuristic.hpp"
00041 #include "BonOACutGenerator2.hpp"
00042
00043
00044 #include "CbcHeuristicFPump.hpp"
00045 #include "CbcHeuristicGreedy.hpp"
00046
00047 #include "BonAmplTMINLP.hpp"
00048
00049 #include "CglGomory.hpp"
00050
00051
00052
00053
00054
00055
00056
00057
00058 #ifdef COIN_HAS_CPX
00059 #include "OsiCpxSolverInterface.hpp"
00060 #else
00061 #include "OsiCbcSolverInterface.hpp"
00062 #endif
00063
00064
00065 #include "CoinTime.hpp"
00066 #include "FP.hpp"
00067
00068 using namespace Bonmin;
00069
00070 OptParam params;
00071 static double BeginTimeGLOB;
00072
00077 #ifdef COIN_HAS_CPX
00078 int findGoodSolution(OsiCpxSolverInterface &mip, double maxTime,
00079 int& nTotalNodes)
00080 {
00081 double beginTime = CoinCpuTime();
00082 int saveIntSolLim;
00083 int saveNodLim;
00084 double saveTiLim;
00085 CPXCENVptr env2 = mip.getEnvironmentPtr();
00086 CPXCLPptr lp = mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
00087 CPXgetintparam(env2,CPX_PARAM_INTSOLLIM, &saveIntSolLim);
00088 CPXgetintparam(env2,CPX_PARAM_NODELIM, &saveNodLim);
00089 CPXgetdblparam(env2,CPX_PARAM_TILIM, &saveTiLim);
00090
00091 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 10000);
00092 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, params.minNodes_);
00093 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime);
00094
00095 mip.branchAndBound();
00096 nTotalNodes += CPXgetnodecnt(env2,lp);
00097 int mipstat = CPXgetstat(env2,lp);
00098 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
00099 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, params.nodeInterval_);
00100 while (mipstat == CPXMIP_NODE_LIM_INFEAS) {
00101 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime - CoinCpuTime() + beginTime);
00102 mip.branchAndBound();
00103 nTotalNodes += CPXgetnodecnt(env2,lp);
00104 mipstat = CPXgetstat(env2,lp);
00105 }
00106 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, saveIntSolLim);
00107 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, saveNodLim);
00108 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, saveTiLim);
00109 return mipstat;
00110 }
00111 #else
00112 int findGoodSolution(OsiSolverInterface &mip, int&nodeNumber)
00113 {
00114 std::cerr<<"Does not work without CPX"<<std::endl;
00115 throw -1;
00116 }
00117
00118 #endif
00119 struct ResolutionInformation
00120 {
00121 double time;
00122 double nlp_time;
00123 double mip_time;
00124 int n_iterations;
00125 int n_nlp_iterations;
00126 int n_mip_nodes;
00127 ResolutionInformation():time(0.), nlp_time(0.), mip_time(0.), n_iterations(0), n_nlp_iterations(0), n_mip_nodes(0)
00128 {}
00129 ResolutionInformation& operator +=(ResolutionInformation &other)
00130 {
00131 time +=other.time;
00132 nlp_time+=(other.nlp_time);
00133 mip_time+=(other.mip_time);
00134 n_iterations+=(other.n_iterations);
00135 n_nlp_iterations+=other.n_nlp_iterations;
00136 n_mip_nodes+=other.n_mip_nodes;
00137 return *this;
00138 }
00139 };
00140
00141 #if 0
00142 void
00143 writeBoundFiles(AmplInterface& nlp, const double * originalLower, const double * originalUpper)
00144 {
00145 const double * currentLower = nlp.getColLower();
00146 const double * currentUpper = nlp.getColUpper();
00147
00148 CoinRelFltEq eq;
00149 std::string fBoundsName;
00150 nlp.getStrParam(OsiProbName,fBoundsName);
00151 fBoundsName+="_bounds_FP";
00152 std::string fModName = fBoundsName + ".mod";
00153 std::ofstream fBounds;
00154 std::ofstream fMod;
00155 bool hasVarNames = 0;
00156 if(nlp.getVarNames()!=NULL )
00157 hasVarNames=1;
00158 if(hasVarNames)
00159 fMod.open(fModName.c_str());
00160 fBounds.open(fBoundsName.c_str());
00161
00162 for(int i = 0 ; i < nlp.getNumCols() ; i++) {
00163
00164 if(nlp.isInteger(i) && !eq(currentLower[i],originalLower[i])) {
00165 if(hasVarNames)
00166 fMod<<"bounds"<<i<<": "
00167 <<nlp.getVarNames()[i]<<" >= "
00168 <<currentLower[i]<<";\n";
00169 fBounds<<"LO"<<"\t"<<i<<"\t"<<currentLower[i]<<std::endl;
00170 }
00171 if(nlp.isInteger(i) && !eq(currentUpper[i],originalUpper[i])) {
00172 if(hasVarNames)
00173 fMod<<"bounds"<<i<<": "
00174 <<nlp.getVarNames()[i]<<" <= "
00175 <<currentUpper[i]<<";\n";
00176 fBounds<<"UP"<<"\t"<<i<<"\t"<<currentUpper[i]<<std::endl;
00177 }
00178 }
00179 }
00180 #endif
00181
00182 double FP(AmplInterface &nlp,
00183 OsiSolverInterface &linearModel,
00184 int numIntCols, int * inds, double * vals,
00185 double maxTime, int maxIter, ResolutionInformation& info,
00186 double ub, int &provenInfeas, double *&solution)
00187 {
00188 provenInfeas=0;
00189 #ifdef COIN_HAS_CPX
00190
00191 OsiCpxSolverInterface * cpxSi = dynamic_cast<OsiCpxSolverInterface *>
00192 (&linearModel);
00193 OsiCpxSolverInterface &mip = *cpxSi;
00194 #endif
00195
00196 int &nMajorIt = info.n_iterations;
00197 int &nTotalNodes = info.n_mip_nodes;
00198 int &nNlpIterations = info.n_nlp_iterations;
00199 double& bbTime = info.mip_time;
00200 double& nlpTime = info.nlp_time;
00201 info.time = -CoinCpuTime();
00202 double objValue = DBL_MAX;
00203 {
00204 OsiCuts cs;
00205
00206 nlpTime -= CoinCpuTime();
00207 double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
00208 nlpTime += CoinCpuTime();
00209 if(dist < 1e-06)
00210 {
00211
00212 std::cout<<"Feasibility subproblem has objective 0 while problem was claimed infeasible before"<<std::endl
00213 <<"I am confused, exiting with error"<<std::endl;
00214 throw -1;
00215 }
00216 OsiSolverInterface::ApplyCutsReturnCode acRc;
00217
00218
00219
00220
00221
00222
00223 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00224 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00225 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00226 <<" were inconsistent for this problem" <<std::endl;
00227 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00228 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00229 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00230 std::cout << std::endl << std::endl;
00231 }
00232
00233
00234 int numcols = linearModel.getNumCols();
00235 double * saveObj = new double[numcols];
00236 CoinCopyN(linearModel.getObjCoefficients(), numcols, saveObj);
00237
00238
00239
00240
00241 if(ub < 1e100) {
00242
00243 linearModel.setColUpper(numcols-1,ub);
00244
00245
00246
00247
00248
00249 }
00250 for(int i = 0;i< numcols ;i++)
00251 if(!linearModel.isInteger(i))
00252 linearModel.setObjCoeff(i,0.);
00253
00254
00255 int nAddedCuts = 0;
00256
00257 bool solved = 0;
00258 bool numericFailure = 0;
00259 int numScaleIncrease = 0;
00260 while((nMajorIt < maxIter) && (!solved && ( CoinCpuTime() + info.time < maxTime) && !numericFailure && numScaleIncrease < 5)) {
00261
00262 for(int i = 0; i < numIntCols; i++) {
00263 linearModel.setObjCoeff(inds[i],1 - 2* nlp.getColSolution()[inds[i]]);
00264 }
00265 #ifdef COIN_HAS_CPX
00266
00267
00268
00269
00270
00271 #else
00272 linearModel.initialSolve();
00273 CbcModel mip(linearModel);
00274
00275 CbcStrategyDefault defaultStrategy;
00276 mip.setStrategy(defaultStrategy);
00277 mip.solver()->messageHandler()->setLogLevel(0);
00278
00279
00280 mip.setMaximumSeconds(maxTime - CoinCpuTime() - info.time);
00281 mip.setMaximumSolutions(1);
00282 #endif
00283
00284 bbTime -= CoinCpuTime();
00285 #ifdef COIN_HAS_CPX
00286
00287 int mipstat = findGoodSolution(mip, maxTime - CoinCpuTime() - info.time, nTotalNodes);
00288 #else
00289
00290 #endif
00291
00292 bbTime += CoinCpuTime();
00293
00294 OsiCuts cs;
00295 nMajorIt++;
00296
00297 #ifdef COIN_HAS_CPX
00298
00299 const double * colsol = mip.getColSolution();
00300
00301 nTotalNodes += CPXgetnodecnt(cpxSi->getEnvironmentPtr(),cpxSi->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00302
00303
00304 if(mipstat == CPXMIP_INFEASIBLE) {
00305 provenInfeas=1;
00306 info.time += CoinCpuTime();
00307 for(int i = 0 ; i < numcols; i++)
00308 linearModel.setObjCoeff(i, saveObj[i]);
00309 return 1e75;
00310 }
00311 else if(mipstat == CPXMIP_TIME_LIM_FEAS || mipstat == CPXMIP_TIME_LIM_INFEAS) {
00312 provenInfeas= -1;
00313 info.time += CoinCpuTime();
00314 for(int i = 0 ; i < numcols; i++)
00315 linearModel.setObjCoeff(i, saveObj[i]);
00316 return 1e75;
00317 }
00318 #else
00319 if(mip.bestSolution() == NULL)
00320 break;
00321 if(mip.isNodeLimitReached()) {
00322 break;
00323 }
00324
00325
00326 nTotalNodes += mip.getNodeCount();
00327 const double * colsol = mip.bestSolution();
00328 #endif
00329
00330 for(int i = 0 ; i < numIntCols ; i++) {
00331 vals[i] = floor(colsol[inds[i]] + 0.5);
00332 vals[i] = max(mip.getColLower()[inds[i]],vals[i]);
00333 vals[i] = min(mip.getColUpper()[inds[i]],vals[i]);
00334 }
00335 nlpTime -= CoinCpuTime();
00336 double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
00337 nlpTime += CoinCpuTime();
00338 nNlpIterations += nlp.getIterationCount();
00339 if(!nlp.isProvenOptimal()) {
00340 std::cout<<"?????"<<std::endl;
00341 throw -1;
00342 }
00343 nlpTime += CoinCpuTime();
00344 if(dist < 1e-04)
00345 {
00346 solved = 1;
00347 double norm_inf = DBL_MAX;
00348 for(int i = 0 ; i < nlp.getNumCols() && solved; i++)
00349 {
00350 if(nlp.isInteger(i)) {
00351 const double &value = nlp.getColSolution()[i];
00352 double IIf = fabs( value- floor(value + 0.5));
00353 norm_inf = min(norm_inf, IIf);
00354 if(fabs( value- floor(value + 0.5)) > 1e-06) {
00355 std::cout<<"Variable "<<i<<" has integer infeasibility : "<<IIf<<std::endl;
00356 solved=0;
00357 numericFailure = 1;
00358 #if 0
00359
00360 numScaleIncrease++;
00361 nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
00362 for(int i = 0 ; i < numIntCols ; i++) {
00363 nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
00364 nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
00365 }
00366 #endif
00367
00368 }
00369 }
00370 }
00371 std::cout<<"Found a solution with maximal integer infeasibility of "<<norm_inf<<std::endl;
00372 }
00373
00374 OsiSolverInterface::ApplyCutsReturnCode acRc;
00375 std::cout<<"Cut violation :"<<cs.rowCut(0).violated(colsol)
00376 <<std::endl;
00377
00378
00379 linearModel.writeMps("test1");
00380 acRc = linearModel.applyCuts(cs);
00381 linearModel.writeMps("test2");
00382 nAddedCuts += cs.sizeRowCuts();
00383
00384 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00385 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00386 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00387 <<" were inconsistent for this problem" <<std::endl;
00388 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00389 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00390 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00391 std::cout << std::endl << std::endl;
00392
00393
00394 if(solved) {
00395
00396 nlp.setColSolution(nlp.getColSolution());
00397 nlp.setRowPrice(nlp.getRowPrice());
00398 nlp.solver()->enableWarmStart();
00399
00400 for(int i = 0; i < numIntCols; i++) {
00401 nlp.setColLower(inds[i], vals[i]);
00402 nlp.setColUpper(inds[i], vals[i]);
00403 }
00404
00405 nlp.initialSolve();
00406 if(nlp.isProvenOptimal()) {
00407 OsiCuts cs;
00408 nlp.getOuterApproximation(cs,1, NULL, true);
00409 linearModel.applyCuts(cs);
00410 nAddedCuts += cs.sizeRowCuts();
00411 std::cout<<" FP found easible solution of value "<<nlp.getObjValue()<<" in "<< CoinCpuTime() - BeginTimeGLOB<<" seconds, "<<nMajorIt<<" major iterations, took"
00412 <<nTotalNodes<<" nodes."<<std::endl;
00413 objValue = nlp.getObjValue();
00414 if (solution==NULL) solution = new double[numcols];
00415 CoinCopyN(nlp.getColSolution(), numcols, solution);
00416 }
00417 else {
00418 std::cout<<" FP converged in "<<CoinCpuTime()- BeginTimeGLOB<<" seconds, "<<nMajorIt<<" major iterations, took"
00419 <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
00420 std::cout<<"Increasing scaling of objective and restarting"<<std::endl;
00421 solved = 0;
00422 numScaleIncrease++;
00423 numericFailure = 1;
00424 #if 0
00425
00426 nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
00427 for(int i = 0 ; i < numIntCols ; i++) {
00428 nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
00429 nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
00430 }
00431 #endif
00432
00433 }
00434
00435 }
00436
00437 }
00438
00439 if(!solved) {
00440 if(numericFailure)
00441 std::cout<<"FP aborted because of a numberical failure : "<<CoinCpuTime() - BeginTimeGLOB<<", "<<nMajorIt<<" major iterations, took"
00442 <<nTotalNodes<<" nodes."<<std::endl;
00443 else
00444 std::cout<<"FP aborted on time limit elapsed time : "<<CoinCpuTime() - BeginTimeGLOB<<", "<<nMajorIt<<" major iterations, took"
00445 <<nTotalNodes<<" nodes."<<std::endl;
00446 provenInfeas = -1;
00447
00448 }
00449 std::cout<<"Nlp solve time : "<<nlpTime<<" B&B solve time : "<<bbTime<<std::endl;
00450
00451 for(int i = 0 ; i < numcols; i++)
00452 linearModel.setObjCoeff(i, saveObj[i]);
00453 delete [] saveObj;
00454 info.time += CoinCpuTime();
00455 return objValue;
00456 }
00457
00458
00459 #if 0 // Old enhanced OA code without the iFP at beginning
00460 int main3 (int argc, char *argv[])
00461 {
00462 std::cout.precision(11);
00463 using namespace Ipopt;
00464
00465
00466 char * pbName = new char[strlen(argv[1])+1];
00467 strcpy(pbName, argv[1]);
00468 BonminAmplInterface solver1(argv);
00469 solver1.turnOnSolverOutput();
00470 bool doFp=true;
00471
00472 #ifdef COIN_HAS_CPX
00473
00474 OsiCpxSolverInterface solver;
00475 OsiCpxSolverInterface &mip = solver;
00476 #else
00477
00478 OsiClpSolverInterface solver;
00479 #endif
00480
00481 solver.messageHandler()->setLogLevel(0);
00482
00483
00484
00485 BeginTimeGLOB= CoinCpuTime();
00486 double bbTime = 0;
00487 double nlpTime = - BeginTimeGLOB;
00488
00489 ResolutionInformation FP_infos;
00490 solver1.extractLinearRelaxation(solver, 1);
00491 nlpTime += CoinCpuTime();
00492
00493 #ifdef EXTRACT_LIN_REL
00494
00495 std::string linearName="Lin";
00496 linearName += pbName;
00497 solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00498 #endif
00499
00500 int numIntCols = 0;
00501 int * inds = new int[solver1.getNumCols()];
00502 double * x = new double[solver1.getNumCols()];
00503
00504
00505 for(int i = 0; i < solver1.getNumCols(); i++) {
00506 if(solver1.isInteger(i)) {
00507 inds[numIntCols++] = i;
00508 }
00509 }
00510
00511 double lb = solver1.getObjValue();
00512 double ub = DBL_MAX;
00513
00514
00515 params.maxTime_ = 2.*3600.;
00516
00517
00518 int nMajorIt = 0;
00519 int nTotalNodes = 0;
00520 int nTimesFPCalled =0;
00521 double precision = 1e-02;
00522
00523 double firstIterationTime;
00524 solver.messageHandler()->setLogLevel(0);
00525 int numNotFound = 0;
00526 bool minutePassed = 0;
00527 while( (ub - lb) > precision && (CoinCpuTime() - BeginTimeGLOB < params.maxTime_) ) {
00528 solver.resolve();
00529 #ifndef COIN_HAS_CPX
00530
00531 CbcModel mip(solver);
00532
00533 CbcStrategyDefault defaultStrategy(1,8,4);
00534 mip.setStrategy(defaultStrategy);
00535
00536 mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
00537 mip.setMaximumSolutions(3);
00538 mip.setCutoff(ub);
00539 mip.setLogLevel(1);
00540 mip.solver()->messageHandler()->setLogLevel(0);
00541
00542 #else
00543
00544
00545 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
00546 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, ub);
00547
00548 #endif
00549
00550 bbTime -= CoinCpuTime();
00551
00552
00553
00554
00555
00556
00557 mip.branchAndBound();
00558 bbTime += CoinCpuTime();
00559 if(nMajorIt==0) {
00560 firstIterationTime = bbTime;
00561 }
00562 #ifdef EXTRACT_LIN_REL
00563 std::string oa_name=linearName;
00564 std::ostringstream t;
00565 t<<nMajorIt;
00566 oa_name+=t.str();
00567 #endif
00568
00569 OsiCuts cs;
00570 #ifdef COIN_HAS_CPX
00571
00572 const double *colsol=mip.getColSolution();
00573 nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00574 double new_lb = 0;
00575 int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00576
00577 if((mipstat != CPXMIP_OPTIMAL && mipstat != CPXMIP_OPTIMAL_TOL)) {
00578 int status = CPXgetbestobjval(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL), &new_lb);
00579 lb=max(lb,new_lb);
00580 if(status)
00581 throw CoinError("Error in getting CPLEX best bound","IpCbcOACutGenerator2","siBestObj");
00582 }
00583 else
00584 lb = mip.getObjValue();
00585
00586 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, mip.getInfinity());
00587
00588 #else
00589
00590
00591 nTotalNodes += mip.getNodeCount();
00592 const double * colsol = mip.bestSolution();
00593 lb = mip.getBestPossibleObjValue();
00594 if(mip.bestSolution() == NULL)
00595 break;
00596 #endif
00597
00598 nMajorIt++;
00599
00600 std::cout<<"Found new lower bound of : "<<lb<<std::endl;
00601 if ((ub - lb) < precision)
00602 break;
00603 for(int i = 0 ; i < numIntCols ; i++) {
00604 x[i] = max(solver.getColLower()[inds[i]],colsol[inds[i]]);
00605 x[i] = min(solver.getColUpper()[inds[i]],x[i]);
00606 solver1.setColLower(inds[i], x[i]);
00607 solver1.setColUpper(inds[i], x[i]);
00608 }
00609 solver1.turnOnSolverOutput();
00610 solver1.initialSolve();
00611 if(solver1.isProvenOptimal()) {
00612 std::cout<<pbName<<" OA found easible solution of value "
00613 <<solver1.getObjValue()<<" in "
00614 <<CoinCpuTime() -BeginTimeGLOB<<" seconds, "
00615 <<nMajorIt<<" major iterations, took"
00616 <<nTotalNodes<<" nodes."<<std::endl;
00617 if(ub > 1e100)
00618 std::cout<<"First solution "<< solver1.getObjValue()<<" after "
00619 <<CoinCpuTime() -BeginTimeGLOB<<" secs."<<std::endl;
00620 if(!minutePassed && (CoinCpuTime() -BeginTimeGLOB) > 61.) {
00621 minutePassed = true;
00622 std::cout<<"Solution at minute "<< ub<<" after "
00623 <<CoinCpuTime() -BeginTimeGLOB<<" secs."<<std::endl;
00624 }
00625
00626 OsiCuts cs;
00627 solver1.getOuterApproximation(cs,1);
00628 OsiSolverInterface::ApplyCutsReturnCode acRc;
00629 ub = min (solver1.getObjValue(), ub);
00630 acRc = solver.applyCuts(cs);
00631
00632 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00633 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00634 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00635 <<" were inconsistent for this problem" <<std::endl;
00636 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00637 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00638 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00639 std::cout << std::endl << std::endl;
00640
00641 }
00642 else if(solver1.isAbandoned() || solver1.isIterationLimitReached()) {
00643 TNLPSolver::UnsolvedError * E=solver1.newUnsolvedError(0,argv[1],"not solved");
00644 E->writeDiffFiles();
00645 solver1.turnOnSolverOutput();
00646 solver1.initialSolve();
00647
00648 std::cerr<<"Error"<<std::endl;
00649 throw -1;
00650 }
00651
00652 else {
00653 if (doFp)
00654 {
00655
00656 for(int i = 0 ; i < numIntCols ; i++)
00657 {
00658 solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
00659 solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
00660 }
00661 ResolutionInformation infos;
00662 nTimesFPCalled++;
00663 int provenInfeas =0;
00664
00665
00666
00667 double fpTime = min( params.maxTime_ - CoinCpuTime() + BeginTimeGLOB,120.);
00668 int fpMaxIter = 30;
00669 ub = min(ub, FP(solver1, solver,numIntCols, inds, x,fpTime, fpMaxIter, infos, ub*(1-1e-03), provenInfeas, solution));
00670 if(provenInfeas == 1)
00671 lb = mip.getInfinity();
00672 if(provenInfeas == -1)
00673 {
00674 numNotFound++;
00675 if(numNotFound==2)
00676 doFp = 0;
00677 }
00678
00679
00680 FP_infos+=infos;
00681 }
00682 else {
00683 std::cout<<"Adding feasibility cuts based on 1-norm of"
00684 <<"constraint satisfaction"<<std::endl;
00685
00686 solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs);
00687 OsiSolverInterface::ApplyCutsReturnCode acRc;
00688 acRc = solver.applyCuts(cs);
00689
00690 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00691 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00692 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00693 <<" were inconsistent for this problem" <<std::endl;
00694 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00695 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00696 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00697 std::cout << std::endl << std::endl;
00698 }
00699 }
00700
00701 }
00702 double time = CoinCpuTime() - BeginTimeGLOB;
00703 std::cout<<pbName<<" FP driven OA found optimal solution of value "<<ub<<" (lower bound is "<<lb<<") in "<<time<<" seconds, "<<nMajorIt<<" major iterations, "
00704 <<nTimesFPCalled<<" calls to FP ( "
00705 <<FP_infos.n_iterations <<" it, and "<< FP_infos.time
00706 <<"sec) overall took "
00707 <<nTotalNodes + FP_infos.n_mip_nodes<<" nodes."<<std::endl;
00708 std::cout<<"Nlp solve time : "<<nlpTime + FP_infos.nlp_time<<" B&B solve time : "<<bbTime + FP_infos.mip_time<<std::endl;
00709 delete[] inds;
00710 delete[] x;
00711 delete [] pbName;
00712 return 0;
00713 }
00714 #endif
00715
00716
00717
00718
00720 int enhancedOA(AmplInterface & solver1, bool doFp,
00721 double *& solution)
00722 {
00723 bool nonConvex = 0;
00724 std::string pbName;
00725 solver1.getStrParam(OsiProbName, pbName);
00726 #ifdef COIN_HAS_CPX
00727
00728 OsiCpxSolverInterface solver;
00729 OsiCpxSolverInterface &mip = solver;
00730 #else
00731
00732 OsiClpSolverInterface solver;
00733 #endif
00734
00735 solver.messageHandler()->setLogLevel(0);
00736
00737
00738
00739 BeginTimeGLOB= CoinCpuTime();
00740 double bbTime = 0;
00741 double nlpTime = - BeginTimeGLOB;
00742
00743 ResolutionInformation FP_infos;
00744 solver1.extractLinearRelaxation(solver, !nonConvex);
00745
00746 nlpTime += CoinCpuTime();
00747 int numcols = solver.getNumCols();
00748 double * saveObj = new double[numcols];
00749 CoinCopyN(solver.getObjCoefficients(), numcols, saveObj);
00750
00751 #ifdef EXTRACT_LIN_REL
00752 {
00753 std::string linearName="Lin";
00754 linearName += basename(pbName);
00755 solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00756 }
00757 #endif
00758
00759 int numIntCols = 0;
00760 int * inds = new int[solver1.getNumCols()];
00761 double * x = new double[solver1.getNumCols()];
00762
00763
00764 for(int i = 0; i < solver1.getNumCols(); i++) {
00765 if(solver1.isInteger(i)) {
00766 inds[numIntCols++] = i;
00767 x[i]=0.5;
00768 }
00769 }
00770
00771 double lb = solver1.getObjValue();
00772 double ub = DBL_MAX;
00773
00774
00775 double maxFpTime = 60.;
00776
00777
00778
00779 int nMajorIt = 0;
00780 int nTotalNodes = 0;
00781 int nTimesFPCalled =0;
00782 int nTimesInitFPCalled =0;
00783 int nFPIterations =0;
00784 double precision = 1e-04;
00785 bool solved = 0;
00786 double firstIterationTime;
00787 solver.messageHandler()->setLogLevel(0);
00788 int numNotFound = 0;
00789
00790 bool feasible = 1;
00791
00792 double FPTime = 0.;
00793
00794 if(doFp)
00795 {
00796
00797 while(ub-lb> precision && feasible && ( CoinCpuTime() - BeginTimeGLOB < maxFpTime)) {
00798 nTimesInitFPCalled++;
00799 int nIt = 0;
00800 solved = 0;
00801 while(feasible & !solved &&
00802 (CoinCpuTime() - BeginTimeGLOB < maxFpTime)
00803
00804 ) {
00805 nIt++;
00806 nFPIterations++;
00807 for(int i = 0 ; i < numIntCols; i++) {
00808 solver.setObjCoeff(inds[i],1 - 2* solver1.getColSolution()[inds[i]]);
00809 }
00810
00811 if(!nonConvex) {
00812 int sign = ub > 0 ? -1 : 1;
00813 solver.setColUpper(solver.getNumCols()-1,ub*(1 + sign *1e-03));
00814 }
00815 #ifndef COIN_HAS_CPX
00816 solver.initialSolve();
00817 solver.messageHandler()->setLogLevel(0);
00818
00819 CbcStrategyDefault defaultStrategy;
00820 CbcModel mip(solver);
00821 mip.setStrategy(defaultStrategy);
00822 mip.solver()->messageHandler()->setLogLevel(0);
00823 mip.setLogLevel(0);
00824 mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() - BeginTimeGLOB);
00825 mip.setMaximumSolutions(1);
00826 #else
00827
00828
00829 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM,
00830 maxFpTime - CoinCpuTime() + BeginTimeGLOB);
00831
00832 #endif
00833
00834 bbTime -= CoinCpuTime();
00835
00836 #ifdef COIN_HAS_CPX
00837
00838 int mipstat = findGoodSolution(mip, maxFpTime - CoinCpuTime() - BeginTimeGLOB, nTotalNodes);
00839 #else
00840
00841 mip.branchAndBound();
00842 #endif
00843
00844 bbTime += CoinCpuTime();
00845
00846
00847 OsiCuts cs;
00848 #ifdef COIN_HAS_CPX
00849
00850 const double *colsol=mip.getColSolution();
00851
00852
00853 if(mipstat == CPXMIP_INFEASIBLE || mipstat == CPXMIP_INForUNBD) {
00854 std::cout<<"I found an Infeasible mip"<<std::endl;
00855 feasible = 0;
00856 break;
00857 }
00858
00859 #else
00860 if(mip.bestSolution() == NULL)
00861 break;
00862 if(mip.isNodeLimitReached()) {
00863 break;
00864 }
00865
00866
00867 nTotalNodes += mip.getNodeCount();
00868 const double * colsol = mip.bestSolution();
00869 #endif
00870
00871 nMajorIt++;
00872 for(int i = 0 ; i < numIntCols ; i++) {
00873 x[i] = max(solver1.getColLower()[inds[i]],colsol[inds[i]]);
00874 x[i] = min(solver1.getColUpper()[inds[i]],x[i]);
00875
00876 }
00877 nlpTime -= CoinCpuTime();
00878 double dist = solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
00879 nlpTime += CoinCpuTime();
00880 std::cout<<"Dist : "<<dist<<std::endl;
00881 if(dist < 1e-04)
00882 {
00883 solved = 1;
00884 for(int i = 0 ; i < solver1.getNumCols() && solved; i++)
00885 {
00886 if(solver1.isInteger(i)) {
00887 const double &value = solver1.getColSolution()[i];
00888 if(fabs( value- floor(value + 0.5)) > 1e-04) {
00889 solved = 0;
00890 std::cout<<"Variable "<<i<<" has integer infeasibility : "<<fabs( value- floor(value + 0.5))<<std::endl;
00891 }
00892 }
00893 }
00894
00895 }
00896
00897 OsiSolverInterface::ApplyCutsReturnCode acRc;
00898 acRc = solver.applyCuts(cs);
00899
00900 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00901 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00902 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00903 <<" were inconsistent for this problem" <<std::endl;
00904 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00905 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00906 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00907 std::cout << std::endl << std::endl;
00908 }
00909 double time = CoinCpuTime() - BeginTimeGLOB;
00910 if(solved) {
00911
00912 for(int i = 0; i < numIntCols; i++) {
00913 solver1.setColLower(inds[i], x[i]);
00914 solver1.setColUpper(inds[i], x[i]);
00915 }
00916 solver1.turnOnSolverOutput();
00917 solver1.initialSolve();
00918 if(solver1.isProvenOptimal()) {
00919 ub = min(solver1.getObjValue(), ub);
00920 if (solution==NULL) solution = new double[numcols];
00921 CoinCopyN(solver1.getColSolution(), numcols, solution);
00922 std::cout<<pbName<<" FP found easible solution of value "
00923 <<solver1.getObjValue()<<" in "
00924 <<CoinCpuTime() - BeginTimeGLOB<<" seconds, "
00925 <<nMajorIt<<" major iterations, took"
00926 <<nTotalNodes<<" nodes."<<std::endl;
00927 if(nTimesInitFPCalled == 1)
00928 std::cout<<"FIRST sol found in : "<<time<<" seconds value is "<<solver1.getObjValue()
00929 <<" iterations "<<nFPIterations<<std::endl;
00930
00931 solved = 1;
00932 OsiCuts cs;
00933 solver1.getOuterApproximation(cs,1, NULL, true);
00934 OsiSolverInterface::ApplyCutsReturnCode acRc;
00935 acRc = solver.applyCuts(cs);
00936 for(int i = 0 ; i < numIntCols ; i++) {
00937 solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
00938 solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
00939 }
00940 solver1.initialSolve();
00941 }
00942 else
00943 {
00944 std::cout<<pbName<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
00945 <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
00946 CoinPackedVector v;
00947 double rhs = 1.;
00948 for(int i = 0 ; i < numIntCols ; i++) {
00949 v.insert(inds[i], -(2*x[i] - 1));
00950 rhs -= x[i];
00951 }
00952 OsiCuts cs;
00953 OsiRowCut cut;
00954 cut.setRow(v);
00955 cut.setLb(rhs);
00956 cut.setUb(1e300);
00957 cut.setGloballyValid();
00958 cs.insert(cut);
00959
00960 OsiSolverInterface::ApplyCutsReturnCode acRc;
00961 acRc = solver.applyCuts(cs);
00962
00963 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00964 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00965 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
00966 <<" were inconsistent for this problem" <<std::endl;
00967 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00968 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00969 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00970 std::cout << std::endl << std::endl;
00971
00972 }
00973 }
00974 else if (feasible) {
00975
00976 solver1.initialSolve();
00977 nIt = 0;
00978 }
00979 if(nonConvex)
00980 break;
00981 }
00982
00983
00984 for(int i = 0 ; i < numcols; i++)
00985 solver.setObjCoeff(i, saveObj[i]);
00986
00987 solver.setColUpper(solver.getNumCols()-1,ub);
00988
00989 FPTime = CoinCpuTime() - BeginTimeGLOB;
00990 std::cout<<"Best known sol after : "<<FPTime<<" is "<<ub<<std::endl;
00991
00992 #ifdef EXTRACT_LIN_REL
00993
00994 {
00995 std::string linearName="Lin1Min";
00996 linearName += basename(pbName);
00997 solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00998 }
00999 #endif
01000
01001 solver.setColUpper(solver.getNumCols()-1,ub);
01002
01003 }
01004 while( (ub - lb) > precision && (CoinCpuTime() - BeginTimeGLOB < params.maxTime_) ) {
01005 #ifdef EXTRACT_LIN_REL
01006 {
01007 if(CoinCpuTime() - BeginTimeGLOB > 5. * 60.)
01008 {
01009 std::string linearName="Lin5min";
01010 linearName += basename(pbName);
01011 solver.writeMpsNative(linearName.c_str(),NULL,NULL);
01012 return 0;
01013 }
01014 }
01015 #endif
01016
01017 solver.resolve();
01018 #ifndef COIN_HAS_CPX
01019
01020 solver.messageHandler()->setLogLevel(0);
01021 CbcModel mip(solver);
01022
01023 CbcStrategyDefault defaultStrategy(1,8,4);
01024 mip.setStrategy(defaultStrategy);
01025
01026 mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
01027 mip.setMaximumSolutions(3);
01028 mip.setCutoff(ub);
01029 mip.setLogLevel(1);
01030 mip.solver()->messageHandler()->setLogLevel(0);
01031
01032 #else
01033
01034
01035 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
01036 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, ub);
01037
01038 #endif
01039
01040 bbTime -= CoinCpuTime();
01041
01042
01043
01044
01045
01046 mip.branchAndBound();
01047 bbTime += CoinCpuTime();
01048 if(nMajorIt==0) {
01049 firstIterationTime = bbTime;
01050 }
01051
01052 OsiCuts cs;
01053 #ifdef COIN_HAS_CPX
01054
01055 const double *colsol=mip.getColSolution();
01056 nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01057 double new_lb = 0;
01058 int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01059
01060 if((mipstat != CPXMIP_OPTIMAL && mipstat != CPXMIP_OPTIMAL_TOL)) {
01061 int status = CPXgetbestobjval(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL), &new_lb);
01062 lb=max(lb,new_lb);
01063 if(status)
01064 throw CoinError("Error in getting CPLEX best bound","IpCbcOACutGenerator2","siBestObj");
01065 }
01066 else
01067 lb = mip.getObjValue();
01068
01069 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, mip.getInfinity());
01070
01071 #else
01072
01073
01074 nTotalNodes += mip.getNodeCount();
01075 const double * colsol = mip.bestSolution();
01076 lb = mip.getBestPossibleObjValue();
01077 if(mip.bestSolution() == NULL)
01078 break;
01079 #endif
01080
01081 nMajorIt++;
01082
01083 std::cout<<"Found new lower bound of : "<<lb<<std::endl;
01084 if ((ub - lb) < precision)
01085 break;
01086 bool same = 1;
01087 for(int i = 0 ; i < numIntCols ; i++) {
01088 double value;
01089 value = max(solver.getColLower()[inds[i]],colsol[inds[i]]);
01090 value = min(solver.getColUpper()[inds[i]],value);
01091 if(fabs(value - x[i])>0.0001) {
01092 same = 0;
01093 }
01094 x[i] = value;
01095 solver1.setColLower(inds[i], x[i]);
01096 solver1.setColUpper(inds[i], x[i]);
01097 }
01098 if(same) {
01099 std::cout<<"Converged on same solution"<<std::endl;
01100 lb += 1e10;
01101 break;
01102 }
01103 solver1.turnOnSolverOutput();
01104 solver1.initialSolve();
01105 if(solver1.isProvenOptimal()) {
01106 if (solution==NULL) solution = new double[numcols];
01107 CoinCopyN(solver1.getColSolution(), numcols, solution);
01108 std::cout<<pbName<<" OA found easible solution of value "
01109 <<solver1.getObjValue()<<" in "
01110 <<CoinCpuTime() -BeginTimeGLOB<<" seconds, "
01111 <<nMajorIt<<" major iterations, took"
01112 <<nTotalNodes<<" nodes."<<std::endl;
01113
01114 OsiCuts cs;
01115 solver1.getOuterApproximation(cs,1, false, true);
01116 OsiSolverInterface::ApplyCutsReturnCode acRc;
01117 ub = min (solver1.getObjValue(), ub);
01118 acRc = solver.applyCuts(cs);
01119
01120 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01121 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01122 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
01123 <<" were inconsistent for this problem" <<std::endl;
01124 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01125 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01126 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01127 std::cout << std::endl << std::endl;
01128
01129 }
01130
01131
01132 else {
01133 if (doFp)
01134 {
01135
01136 for(int i = 0 ; i < numIntCols ; i++)
01137 {
01138 solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
01139 solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
01140 }
01141 ResolutionInformation infos;
01142 nTimesFPCalled++;
01143 int provenInfeas =0;
01144
01145
01146
01147 double fpTime = min( params.maxTime_ - CoinCpuTime() + BeginTimeGLOB,120.);
01148 int fpMaxIter = 5;
01149 ub = min(ub, FP(solver1, solver,numIntCols, inds, x,fpTime, fpMaxIter, infos, ub*(1-1e-03), provenInfeas, solution));
01150 if(provenInfeas == 1)
01151 lb = mip.getInfinity();
01152 if(provenInfeas == -1)
01153 {
01154 numNotFound++;
01155
01156
01157 }
01158
01159
01160 FP_infos+=infos;
01161 }
01162 else {
01163 std::cout<<"Adding feasibility cuts based on 1-norm of"
01164 <<"constraint satisfaction"<<std::endl;
01165
01166 solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
01167 OsiSolverInterface::ApplyCutsReturnCode acRc;
01168 acRc = solver.applyCuts(cs);
01169
01170 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01171 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01172 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
01173 <<" were inconsistent for this problem" <<std::endl;
01174 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01175 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01176 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01177 std::cout << std::endl << std::endl;
01178 }
01179 }
01180
01181 }
01182 std::string algoName=(doFp)?" enhanced OA":" OA";
01183 double time = CoinCpuTime() - BeginTimeGLOB;
01184 std::cout<<pbName<<algoName<<" found optimal solution of value "<<ub<<" (lower bound is "<<lb<<") in "<<time<<" seconds, "<<nMajorIt<<" major iterations, "<<std::endl;
01185 if(doFp)
01186 {
01187 std::cout<<nTimesInitFPCalled<<" calls to initial FP ( "
01188 <<nFPIterations <<" it, and "<< FPTime
01189 <<"sec)"
01190 <<std::endl
01191 <<nTimesFPCalled<<"subsequent calls to FP ( "
01192 <<FP_infos.n_iterations <<" it, and "<< FP_infos.time
01193 <<"sec)"<<std::endl;
01194 }
01195 std::cout<<"Nlp solve time : "<<nlpTime + FP_infos.nlp_time<<" B&B solve time : "<<bbTime + FP_infos.mip_time<<std::endl;
01196 delete[] inds;
01197 delete[] x;
01198 return 0;
01199 }
01200
01201
01202
01204 int iteratedFP (AmplInterface& solver1, bool standAlone,
01205 double * &solution)
01206 {
01207
01208
01209 using namespace Ipopt;
01210 std::string pbName;
01211 solver1.getStrParam(OsiProbName,pbName);
01212
01213 #ifdef COIN_HAS_CPX
01214
01215 OsiCpxSolverInterface solver;
01216 OsiCpxSolverInterface &mip = solver;
01217 #else
01218
01219 OsiClpSolverInterface solver;
01220 #endif
01221
01222 solver.messageHandler()->setLogLevel(0);
01223
01224
01225 double beginTime= CoinCpuTime();
01226 double bbTime = 0;
01227 double nlpTime = - beginTime;
01228
01229 solver1.extractLinearRelaxation(solver, 1);
01230 nlpTime += CoinCpuTime();
01231
01232
01233 int numIntCols = 0;
01234 int numcols = solver1.getNumCols();
01235 int * inds = new int[numcols];
01236 double * x = new double[numcols];
01237
01238 for(int i = 0; i < numcols; i++) {
01239 if(solver1.isInteger(i)) {
01240 inds[numIntCols++] = i;
01241 }
01242 }
01243
01244
01245
01246 for(int i = 0;i<solver.getNumCols();i++)
01247 solver.setObjCoeff(i,0.);
01248
01249
01250
01251
01252
01253 int nMajorIt = 0;
01254 int nTotalNodes = 0;
01255
01256
01257 double lb=-DBL_MAX;
01258 double ub=DBL_MAX;
01259 bool solved = 0;
01260 while(ub-lb>1e-05) {
01261 double ub = 1e100;
01262 bool feasible = 1;
01263 while(feasible && ( CoinCpuTime() - beginTime < params.maxTime_) ) {
01264 int nIt = 0;
01265 while(feasible & !solved && (CoinCpuTime() - beginTime < params.maxTime_)) {
01266 nIt++;
01267 for(int i = 0 ; i < numIntCols; i++) {
01268 solver.setObjCoeff(inds[i],1 - 2* solver1.getColSolution()[inds[i]]);
01269 }
01270
01271 solver.setColUpper(solver.getNumCols()-1,ub);
01272
01273 #ifndef COIN_HAS_CPX
01274
01275 solver.initialSolve();
01276 CbcModel mip(solver);
01277 CbcStrategyDefault defaultStrategy;
01278 mip.setStrategy(defaultStrategy);
01279 mip.solver()->messageHandler()->setLogLevel(0);
01280
01281
01282 mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() +beginTime);
01283 mip.setMaximumSolutions(1);
01284 #else
01285
01286 CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 3);
01287 CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM,
01288 params.maxTime_ - CoinCpuTime() + beginTime);
01289
01290 #endif
01291
01292 bbTime -= CoinCpuTime();
01293 mip.branchAndBound();
01294 bbTime += CoinCpuTime();
01295
01296
01297 OsiCuts cs;
01298 #ifdef COIN_HAS_CPX
01299
01300 const double *colsol=mip.getColSolution();
01301 nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01302 int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01303 if(mipstat == CPXMIP_INFEASIBLE || mipstat == CPXMIP_INForUNBD) {
01304 std::cout<<"Infeasible mip"<<std::endl;
01305 feasible = 0;
01306 break;
01307 }
01308
01309 #else
01310 if(mip.bestSolution() == NULL)
01311 break;
01312 if(mip.isNodeLimitReached()) {
01313 break;
01314 }
01315
01316
01317 nTotalNodes += mip.getNodeCount();
01318 const double * colsol = mip.bestSolution();
01319 #endif
01320
01321 nMajorIt++;
01322 for(int i = 0 ; i < numIntCols ; i++) {
01323 x[i] = max(solver1.getColLower()[inds[i]],colsol[inds[i]]);
01324 x[i] = min(solver1.getColUpper()[inds[i]],x[i]);
01325
01326 }
01327 nlpTime -= CoinCpuTime();
01328 double dist = solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
01329 nlpTime += CoinCpuTime();
01330 std::cout<<"Dist : "<<dist<<std::endl;
01331 if(dist < 1e-05)
01332 {
01333 solved = 1;
01334 for(int i = 0 ; i < numcols && solved; i++)
01335 {
01336 if(solver1.isInteger(i)) {
01337 const double &value = solver1.getColSolution()[i];
01338 if(fabs( value- floor(value + 0.5)) > 1e-04) {
01339 solved = 0;
01340 std::cout<<"Variable "<<i<<" has integer infeasibility : "<<fabs( value- floor(value + 0.5))<<std::endl;
01341 }
01342 }
01343 }
01344 feasible = !solved;
01345 }
01346
01347 OsiSolverInterface::ApplyCutsReturnCode acRc;
01348 acRc = solver.applyCuts(cs);
01349
01350 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01351 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01352 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
01353 <<" were inconsistent for this problem" <<std::endl;
01354 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01355 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01356 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01357 std::cout << std::endl << std::endl;
01358 }
01359 double time = CoinCpuTime() - beginTime;
01360 if(solved) {
01361
01362 for(int i = 0; i < numIntCols; i++) {
01363 solver1.setColLower(inds[i], x[i]);
01364 solver1.setColUpper(inds[i], x[i]);
01365 }
01366 solver1.turnOnSolverOutput();
01367 solver1.initialSolve();
01368 if(solver1.isProvenOptimal()) {
01369 ub = min(solver1.getObjValue() * (1-1e-4), ub);
01370 if (solution==NULL) solution = new double[numcols];
01371 CoinCopyN(solver1.getColSolution(), numcols, solution);
01372 std::cout<<pbName<<" FP found easible solution of value "
01373 <<solver1.getObjValue()<<" in "
01374 <<CoinCpuTime() - beginTime<<" seconds, "
01375 <<nMajorIt<<" major iterations, took"
01376 <<nTotalNodes<<" nodes."<<std::endl;
01377 solved = 1;
01378 OsiCuts cs;
01379 solver1.getOuterApproximation(cs,1, NULL, true);
01380 OsiSolverInterface::ApplyCutsReturnCode acRc;
01381 acRc = solver.applyCuts(cs);
01382 for(int i = 0 ; i < numIntCols ; i++) {
01383 solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
01384 solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
01385 }
01386
01387 }
01388 else
01389 std::cout<<pbName<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
01390 <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
01391
01392 }
01393 else if (feasible) {
01394
01395 solver1.initialSolve();
01396 nIt = 0;
01397 }
01398 }
01399 if(standAlone)
01400 break;
01401 }
01402 int returncode = 0;
01403 if(!standAlone)
01404 {
01405
01406 if(solved) {
01407 std::cout<<"iterated FP finished in : "<<CoinCpuTime() - beginTime
01408 <<", value of optimum "<<ub
01409 <<", "<<nMajorIt<<" major iterations, took"
01410 <<nTotalNodes<<" nodes."<<std::endl;
01411 returncode = 1;
01412 }
01413 else
01414 std::cout<<"iterated FP aborted on time limit elapsed time : "
01415 <<CoinCpuTime() - beginTime<<", "<<nMajorIt
01416 <<" major iterations, took"
01417 <<nTotalNodes<<" nodes."<<std::endl;
01418 returncode = 0;
01419 }
01420 else returncode = 0;
01421 std::cout<<"Nlp solve time : "<<nlpTime
01422 <<" B&B solve time : "<<bbTime<<std::endl;
01423 delete[] inds;
01424 delete[] x;
01425 return returncode;
01426 }
01427
01428
01429
01430 double FPGeneralIntegers(AmplInterface &nlp, OsiSolverInterface &linearModel,
01431 int numIntCols, int * inds, double * vals, double maxTime,
01432 ResolutionInformation& info, double ub, bool &provenInfeas)
01433
01434 {
01435 provenInfeas=0;
01436 #ifdef COIN_HAS_CPX
01437
01438 OsiCpxSolverInterface * cpxSi = dynamic_cast<OsiCpxSolverInterface *>
01439 (&linearModel);
01440 OsiCpxSolverInterface &mip = *cpxSi;
01441 #endif
01442
01443 int &nMajorIt = info.n_iterations;
01444 int &nTotalNodes = info.n_mip_nodes;
01445 int &nNlpIterations = info.n_nlp_iterations;
01446 double& bbTime = info.mip_time;
01447 double& nlpTime = info.nlp_time;
01448 double &time = info.time = -CoinCpuTime();
01449 double objValue = DBL_MAX;
01450 {
01451 OsiCuts cs;
01452
01453 nlpTime -= CoinCpuTime();
01454 double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
01455 nlpTime += CoinCpuTime();
01456 if(dist < 1e-06)
01457 {
01458
01459 std::cout<<"Feasibility subproblem has objective 0 while problem was claimed infeasible before"<<std::endl
01460 <<"I am confused, exiting with error"<<std::endl;
01461 throw -1;
01462 }
01463
01464 OsiSolverInterface::ApplyCutsReturnCode acRc;
01465
01466
01467
01468
01469
01470
01471 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01472 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01473 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
01474 <<" were inconsistent for this problem" <<std::endl;
01475 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01476 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01477 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01478 std::cout << std::endl << std::endl;
01479 }
01480
01481
01482 double * saveObj = new double[linearModel.getNumCols()];
01483 CoinCopyN(linearModel.getObjCoefficients(), linearModel.getNumCols(), saveObj);
01484
01485
01486
01487 int origNumCols = linearModel.getNumCols();
01488 int origNumRows = linearModel.getNumRows();
01489
01490 double * colLb = new double[2*numIntCols+1];
01491 double * colUb = new double[2*numIntCols+1];
01492 double * obj = new double[2*numIntCols];
01493 int k = 0;
01494 CoinPackedVectorBase ** emptyCols = new CoinPackedVectorBase*[2*numIntCols+1];
01495 for(int i = 0; i < nlp.getNumCols(); i++) {
01496 if(nlp.isInteger(i)) {
01497 colLb[k] = colLb[k+1] = 0.;
01498 colUb[k] = colUb[k+1] = nlp.getColUpper()[i] - nlp.getColLower()[i];
01499 obj[k] = obj[k+1] = 1.;
01500 emptyCols[k] = new CoinPackedVector;
01501 emptyCols[k + 1] = new CoinPackedVector;
01502 k+=2;
01503 }
01504 }
01505
01506 linearModel.addCols(2*numIntCols, emptyCols, colLb, colUb, obj);
01507 k = 0;
01508
01509 for(int i = 0; i < nlp.getNumCols(); i++) {
01510 if(nlp.isInteger(i)) {
01511 CoinPackedVector * rowToAdd = (CoinPackedVector *) emptyCols[k];
01512 rowToAdd->insert(i,1.);
01513 rowToAdd->insert(origNumCols + 2*k ,-1.);
01514 rowToAdd->insert(origNumCols + 2*k + 1,1.);
01515 colUb[k] = 0.;
01516 k++;
01517 }
01518 }
01519
01520
01521 if(ub < 1e100) {
01522
01523 linearModel.setColUpper(origNumCols-1,ub);
01524
01525
01526
01527
01528
01529 }
01530 for(int i = 0;i<origNumCols;i++)
01531 linearModel.setObjCoeff(i,0.);
01532 linearModel.addRows(numIntCols, emptyCols, colLb, colUb);
01533
01534 int nAddedCuts = 0;
01535
01536 bool solved = 0;
01537 bool numericFailure = 0;
01538 int numScaleIncrease = 0;
01539 while(!solved && (time + CoinCpuTime() < maxTime) && !numericFailure && numScaleIncrease < 5 ) {
01540 for(int i = 0; i < numIntCols; i++) {
01541 linearModel.setRowBounds(origNumRows + i, nlp.getColSolution()[inds[i]], nlp.getColSolution()[inds[i]]);
01542 }
01543 #ifdef COIN_HAS_CPX
01544
01545 CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
01546 CPXsetdblparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime - CoinCpuTime() + time);
01547
01548 #else
01549
01550 linearModel.initialSolve();
01551 CbcModel mip(linearModel);
01552
01553 CbcStrategyDefault defaultStrategy;
01554 mip.setStrategy(defaultStrategy);
01555 mip.solver()->messageHandler()->setLogLevel(0);
01556
01557 mip.setMaximumSeconds(maxTime - CoinCpuTime() + time);
01558 mip.setMaximumSolutions(1);
01559 #endif
01560
01561 bbTime -= CoinCpuTime();
01562 mip.branchAndBound();
01563 bbTime += CoinCpuTime();
01564
01565 OsiCuts cs;
01566 nMajorIt++;
01567
01568 #ifdef COIN_HAS_CPX
01569
01570 const double * colsol = mip.getColSolution();
01571
01572 nTotalNodes += CPXgetnodecnt(cpxSi->getEnvironmentPtr(),cpxSi->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01573 int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01574 if(mipstat == CPXMIP_INFEASIBLE) {
01575 provenInfeas=1;
01576 return 1e75;
01577 }
01578 #else
01579 if(mip.bestSolution() == NULL)
01580 break;
01581 if(mip.isNodeLimitReached()) {
01582 break;
01583 }
01584
01585
01586 nTotalNodes += mip.getNodeCount();
01587 const double * colsol = mip.bestSolution();
01588 #endif
01589
01590 for(int i = 0 ; i < numIntCols ; i++) {
01591 vals[i] = floor(colsol[inds[i]] + 0.5);
01592 vals[i] = max(mip.getColLower()[inds[i]],vals[i]);
01593 vals[i] = min(mip.getColUpper()[inds[i]],vals[i]);
01594 }
01595 nlpTime -= CoinCpuTime();
01596 double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
01597 nlpTime += CoinCpuTime();
01598 nNlpIterations += nlp.getIterationCount();
01599 if(!nlp.isProvenOptimal()) {
01600 std::cout<<"?????"<<std::endl;
01601 throw -1;
01602 }
01603 nlpTime += CoinCpuTime();
01604 if(dist < 1e-04)
01605 {
01606 solved = 1;
01607 double norm_inf = DBL_MAX;
01608 for(int i = 0 ; i < nlp.getNumCols() && solved; i++)
01609 {
01610 if(nlp.isInteger(i)) {
01611 const double &value = nlp.getColSolution()[i];
01612 double IIf = fabs( value- floor(value + 0.5));
01613 norm_inf = min(norm_inf, IIf);
01614 if(fabs( value- floor(value + 0.5)) > 1e-04) {
01615 std::cout<<"Variable "<<i<<" has integer infeasibility : "<<IIf<<std::endl;
01616 solved=0;
01617 numericFailure = 1;
01618 #if 0
01619
01620 numScaleIncrease++;
01621 nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
01622 for(int i = 0 ; i < numIntCols ; i++) {
01623 nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
01624 nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
01625 }
01626 #endif
01627
01628 }
01629 }
01630 }
01631 std::cout<<"Found a solution with maximal integer infeasibility of "<<norm_inf<<std::endl;
01632 }
01633
01634 OsiSolverInterface::ApplyCutsReturnCode acRc;
01635 std::cout<<"Cut violation :"<<cs.rowCut(0).violated(colsol)
01636 <<std::endl;
01637
01638
01639 linearModel.writeMps("test1");
01640 acRc = linearModel.applyCuts(cs);
01641 linearModel.writeMps("test2");
01642 nAddedCuts += cs.sizeRowCuts();
01643
01644 std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01645 std::cout <<" " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01646 std::cout <<" " <<acRc.getNumInconsistentWrtIntegerModel()
01647 <<" were inconsistent for this problem" <<std::endl;
01648 std::cout <<" " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01649 std::cout <<" " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01650 std::cout <<" " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01651 std::cout << std::endl << std::endl;
01652
01653
01654 if(solved) {
01655
01656 nlp.setColSolution(nlp.getColSolution());
01657 nlp.setRowPrice(nlp.getRowPrice());
01658 nlp.solver()->enableWarmStart();
01659
01660 for(int i = 0; i < numIntCols; i++) {
01661 nlp.setColLower(inds[i], vals[i]);
01662 nlp.setColUpper(inds[i], vals[i]);
01663 }
01664
01665 nlp.initialSolve();
01666 if(nlp.isProvenOptimal()) {
01667 OsiCuts cs;
01668 nlp.getOuterApproximation(cs,1, NULL, true);
01669 linearModel.applyCuts(cs);
01670 nAddedCuts += cs.sizeRowCuts();
01671 std::cout<<" FP found easible solution of value "<<nlp.getObjValue()<<" in "<<time + CoinCpuTime()<<" seconds, "<<nMajorIt<<" major iterations, took"
01672 <<nTotalNodes<<" nodes."<<std::endl;
01673 objValue = nlp.getObjValue();
01674 }
01675 else {
01676 std::cout<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
01677 <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
01678 std::cout<<"Increasing scaling of objective and restarting"<<std::endl;
01679 solved = 0;
01680 numScaleIncrease++;
01681 numericFailure = 1;
01682 #if 0
01683
01684 nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
01685 for(int i = 0 ; i < numIntCols ; i++) {
01686 nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
01687 nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
01688 }
01689 #endif
01690
01691 }
01692
01693 }
01694
01695 }
01696 time+=CoinCpuTime();
01697
01698 if(!solved) {
01699 if(numericFailure)
01700 std::cout<<"FP aborted because of a numberical failure : "<<time<<", "<<nMajorIt<<" major iterations, took"
01701 <<nTotalNodes<<" nodes."<<std::endl;
01702 else
01703 std::cout<<"FP aborted on time limit elapsed time : "<<time<<", "<<nMajorIt<<" major iterations, took"
01704 <<nTotalNodes<<" nodes."<<std::endl;
01705
01706 }
01707 std::cout<<"Nlp solve time : "<<nlpTime<<" B&B solve time : "<<bbTime<<std::endl;
01708
01709
01710
01711 int * colToDelete = new int[2*numIntCols];
01712 int * rowToDelete = new int[numIntCols];
01713 for(int i = 0 ; i < numIntCols ; i++) {
01714 colToDelete[2*i] = origNumCols + 2*i;
01715 colToDelete[2*i +1] = origNumCols + 2*i + 1;
01716 rowToDelete[i] = origNumRows + i;
01717 }
01718
01719
01720 linearModel.deleteCols(2*numIntCols, colToDelete);
01721 linearModel.deleteRows(numIntCols, rowToDelete);
01722
01723
01724 assert(linearModel.getNumCols()==origNumCols);
01725 assert(linearModel.getNumRows()==origNumRows + nAddedCuts);
01726
01727 for(int i = 0 ; i < origNumCols; i++)
01728 linearModel.setObjCoeff(i, saveObj[i]);
01729 delete [] saveObj;
01730 return objValue;
01731 }