00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "CbcModel.hpp"
00013 #include "CoinTime.hpp"
00014 #include "CoinHelperFunctions.hpp"
00015
00016 #include "CouenneExprAux.hpp"
00017 #include "CouenneFeasPump.hpp"
00018 #include "CouenneProblem.hpp"
00019 #include "CouenneProblemElem.hpp"
00020 #include "CouenneCutGenerator.hpp"
00021 #include "CouenneTNLP.hpp"
00022 #include "CouenneFPpool.hpp"
00023 #include "CouenneRecordBestSol.hpp"
00024
00025 #ifdef COIN_HAS_SCIP
00026
00027 #include "scip/scip.h"
00028 #include "scip/cons_linear.h"
00029 #include "scip/scipdefplugins.h"
00030 #endif
00031
00032 using namespace Ipopt;
00033 using namespace Couenne;
00034
00035 void printDist (CouenneProblem *p, double *iSol, double *nSol);
00036 void printCmpSol (int n, double *iSol, double *nSol, int direction);
00037
00038
00039 int CouenneFeasPump::solution (double &objVal, double *newSolution) {
00040
00041 const int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
00042
00043 if (
00044 (CoinCpuTime () > problem_ -> getMaxCpuTime ()) ||
00045 ((numberSolvePerLevel_ >= 0) &&
00046 (CoinDrand48 () > 1. / CoinMax
00047 (1., (double) ((CoinMax (0, depth - numberSolvePerLevel_ + 1)) *
00048 (depth - numberSolvePerLevel_ + 1))))))
00049 return 0;
00050
00051 problem_ -> Jnlst () -> Printf (J_ERROR, J_NLPHEURISTIC, "FP: BEGIN\n");
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 if (!nlp_)
00062 nlp_ = new CouenneTNLP (problem_);
00063
00064 problem_ -> domain () -> push (problem_ -> nVars (),
00065 problem_ -> domain () -> x (),
00066 problem_ -> domain () -> lb (),
00067 problem_ -> domain () -> ub ());
00068
00069
00070
00071
00072 nlp_ -> setInitSol (problem_ -> domain () -> x ());
00073
00075
00076 if ((multHessNLP_ > 0.) ||
00077 (multHessMILP_ > 0.))
00078 nlp_ -> getSaveOptHessian () = true;
00079
00080 if (problem_ -> Jnlst () -> ProduceOutput
00081 (J_WARNING, J_NLPHEURISTIC)) {
00082
00083 printf ("initial NLP:\n");
00084 problem_ -> print ();
00085 printf ("---------------------\n");
00086 }
00087
00088
00089 ApplicationReturnStatus status = app_ -> OptimizeTNLP (nlp_);
00090
00092
00093 problem_ -> domain () -> pop ();
00094
00095 if ((status != Solve_Succeeded) &&
00096 (status != Solved_To_Acceptable_Level))
00097
00098 problem_ -> Jnlst () -> Printf
00099 (J_ERROR, J_NLPHEURISTIC, "Feasibility Pump: error in initial NLP problem\n");
00100
00101 if ((multHessNLP_ > 0.) ||
00102 (multHessMILP_ > 0.))
00103 nlp_ -> getSaveOptHessian () = false;
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 CouNumber
00135 *nSol = NULL,
00136 *iSol = NULL,
00137 *best = NULL;
00138
00139 int
00140 niter = 0,
00141 nsuciter = 0,
00142 retval = 0,
00143 nSep = 0,
00144 objInd = problem_ -> Obj (0) -> Body () -> Index ();
00145
00147
00148
00149
00150
00151
00152
00153
00155
00156
00157
00158
00159
00160 problem_ -> domain () -> push (model_ -> solver ());
00161
00162 expression *originalObjective = problem_ -> Obj (0) -> Body ();
00163
00164 do {
00165
00166
00167
00168
00169
00170
00171
00172 double z = solveMILP (nSol, iSol, niter, &nsuciter);
00173
00174
00175
00176 if (!iSol || z >= COIN_DBL_MAX/2) {
00177
00178 problem_ -> Jnlst () -> Printf
00179 (Ipopt::J_ERROR, J_NLPHEURISTIC,
00180 "FP: breaking out of loop upon %p==NULL or %.3f large\n", iSol, z);
00181
00182 break;
00183 }
00184
00185 bool isChecked = false;
00186
00187
00188
00189
00190
00191
00192
00193
00194 CouenneFPsolution checkedSol (problem_, iSol, false);
00195
00196 if (tabuPool_. find (checkedSol) != tabuPool_ . end ()) {
00197
00198
00199
00200
00201 if (tabuMgt_ == FP_TABU_NONE) break;
00202
00203 else if ((tabuMgt_ == FP_TABU_POOL) && !(pool_ -> Set (). empty ())) {
00204
00205
00206 do {
00207
00208
00209 pool_ -> findClosestAndReplace (iSol, nSol, problem_ -> nVars ());
00210
00211 CouenneFPsolution newSol (problem_, iSol);
00212
00213
00214 if (tabuPool_. find(newSol) == tabuPool_ . end ())
00215 break;
00216
00217
00218 if (pool_ -> Set ().empty())
00219 {
00220 delete[] iSol;
00221 iSol = NULL;
00222 }
00223
00224 } while( !pool_ -> Set ().empty() );
00225
00226 } else if (((tabuMgt_ == FP_TABU_CUT) ||
00227 ((pool_ -> Set (). empty ()) && iSol))) {
00228
00229 OsiCuts cs;
00230
00231 problem_ -> domain () -> push (problem_ -> nVars (), iSol, NULL, NULL);
00232 couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL);
00233 problem_ -> domain () -> pop ();
00234
00235 if (cs.sizeRowCuts ()) {
00236
00237 milp_ -> applyCuts (cs);
00238
00239 if (nSep++ < nSepRounds_)
00240 continue;
00241
00242 } else break;
00243
00244 } else if ((tabuMgt_ == FP_TABU_PERTURB) && iSol) {
00245
00246
00247
00248 const CouNumber
00249 *lb = milp_ -> getColLower (),
00250 *ub = milp_ -> getColUpper ();
00251
00252 double
00253 downMoves = 0.,
00254 upMoves = 0.;
00255
00256 int startIndex = (int) floor (CoinDrand48 () * problem_ -> nOrigVars ());
00257
00258 for (int ii=problem_ -> nOrigVars (); ii--; lb++, ub++) {
00259
00260
00261
00262 int i = (startIndex + ii) % problem_ -> nOrigVars ();
00263
00264 if (problem_ -> Var (i) -> isInteger ()) {
00265
00266 double
00267 rnd = CoinDrand48 (),
00268 down = 0.,
00269 up = 1.;
00270
00271
00272
00273
00274
00275 #define RND_DECR_EXPONENT .5
00276
00277 if (iSol [i] >= lb [i] - 1.) down = 1. / pow (1. + (downMoves += 1.), RND_DECR_EXPONENT);
00278 if (iSol [i] <= ub [i] + 1.) up = 1. - 1. / pow (1. + (upMoves += 1.), RND_DECR_EXPONENT);
00279
00280 if (rnd < down) iSol [i] -= 1.;
00281 else if (rnd > up) iSol [i] += 1.;
00282 }
00283 }
00284 }
00285
00286 } else tabuPool_. insert (CouenneFPsolution (problem_, iSol));
00287
00288 #ifdef FM_CHECKNLP2
00289 isChecked = problem_ -> checkNLP2 (iSol, 0, false,
00290 true,
00291 true,
00292 problem_->getFeasTol());
00293 if(isChecked) {
00294 z = problem_->getRecordBestSol()->getModSolVal();
00295 }
00296
00297 #else
00298 isChecked = problem_ -> checkNLP (iSol, z, true);
00299 #endif
00300
00301 if (isChecked) {
00302
00303 problem_ -> Jnlst () -> Printf
00304 (Ipopt::J_ERROR, J_NLPHEURISTIC,
00305 "FP: IP solution is MINLP feasible\n");
00306
00307
00308
00309 retval = 1;
00310 objVal = z;
00311
00312 #ifdef FM_CHECKNLP2
00313 #ifdef FM_TRACE_OPTSOL
00314 problem_->getRecordBestSol()->update();
00315 best = problem_->getRecordBestSol()->getSol();
00316 objVal = problem_->getRecordBestSol()->getVal();
00317 #else
00318 best = problem_->getRecordBestSol()->getModSol(problem_->nVars());
00319 objVal = z;
00320 #endif
00321 #else
00322 #ifdef FM_TRACE_OPTSOL
00323 problem_->getRecordBestSol()->update(iSol, problem_->nVars(),
00324 z, problem_->getFeasTol());
00325 best = problem_->getRecordBestSol()->getSol();
00326 objVal = problem_->getRecordBestSol()->getVal();
00327 #else
00328 best = iSol;
00329 objVal = z;
00330 #endif
00331 #endif
00332
00333 if (z < problem_ -> getCutOff ()) {
00334
00335 problem_ -> setCutOff (objVal);
00336
00337 t_chg_bounds *chg_bds = NULL;
00338
00339 if (objInd >= 0) {
00340
00341 chg_bds = new t_chg_bounds [problem_ -> nVars ()];
00342 chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
00343 }
00344
00345
00346 bool is_still_feas = problem_ -> btCore (chg_bds);
00347
00348
00349
00350 if (chg_bds)
00351 delete [] chg_bds;
00352
00353 if (!is_still_feas)
00354 break;
00355
00356
00357 const CouNumber
00358 *plb = problem_ -> Lb (),
00359 *pub = problem_ -> Ub (),
00360 *mlb = milp_ -> getColLower (),
00361 *mub = milp_ -> getColUpper ();
00362
00363 for (int i=problem_ -> nVars (), j=0; i--; ++j, ++plb, ++pub) {
00364
00365 if (*plb > *mlb++) milp_ -> setColLower (j, *plb);
00366 if (*pub < *mub++) milp_ -> setColUpper (j, *pub);
00367 }
00368 }
00369
00370 break;
00371
00372 } else if (milpCuttingPlane_ == FP_CUT_EXTERNAL ||
00373 milpCuttingPlane_ == FP_CUT_POST) {
00374
00375
00376
00377
00378 OsiCuts cs;
00379
00380 problem_ -> domain () -> push (milp_);
00381 couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL);
00382 problem_ -> domain () -> pop ();
00383
00384 if (cs.sizeRowCuts ()) {
00385
00386
00387
00388 milp_ -> applyCuts (cs);
00389
00390
00391 if (milpCuttingPlane_ == FP_CUT_EXTERNAL &&
00392 nSep++ < nSepRounds_)
00393 continue;
00394 }
00395 }
00396
00397
00398
00399
00400
00401 nSep = 0;
00402
00403
00404
00405
00406
00407
00408 z = solveNLP (iSol, nSol);
00409
00410 if ((nSol && iSol) &&
00411 (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC))) {
00412
00413 double dist = 0.;
00414 int nNonint = 0;
00415
00416 for (int i = 0; i < problem_ -> nVars (); ++i) {
00417
00418 if (problem_ -> Var (i) -> isInteger () &&
00419 (fabs (iSol [i] - ceil (iSol [i] - .5)) > 1e-4))
00420 ++nNonint;
00421
00422 dist +=
00423 (iSol [i] - nSol [i]) *
00424 (iSol [i] - nSol [i]);
00425 }
00426
00427 printf ("FP: after NLP, distance %g, %d nonintegers\n", sqrt (dist), nNonint);
00428 }
00429
00430 if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_NLPHEURISTIC)) {
00431
00432 printDist (problem_, iSol, nSol);
00433 printCmpSol (problem_ -> nVars (), iSol, nSol, 0);
00434 }
00435
00436 if (z > COIN_DBL_MAX/2)
00437 break;
00438
00439
00440
00441 isChecked = false;
00442
00443 if (nSol) {
00444
00445 #ifdef FM_CHECKNLP2
00446 isChecked = problem_->checkNLP2(nSol, 0, false,
00447 true,
00448 true,
00449 problem_->getFeasTol());
00450 if(isChecked) {
00451 z = problem_->getRecordBestSol()->getModSolVal();
00452 }
00453 #else
00454 isChecked = problem_ -> checkNLP (nSol, z, true);
00455 #endif
00456 }
00457
00458 if (nSol && isChecked) {
00459
00460 retval = 1;
00461 objVal = z;
00462
00463 #ifdef FM_CHECKNLP2
00464 #ifdef FM_TRACE_OPTSOL
00465 problem_->getRecordBestSol()->update();
00466 best = problem_->getRecordBestSol()->getSol();
00467 objVal = problem_->getRecordBestSol()->getVal();
00468 #else
00469 best = problem_->getRecordBestSol()->getModSol(problem_ -> nVars ());
00470 objVal = z;
00471 #endif
00472 #else
00473 #ifdef FM_TRACE_OPTSOL
00474 problem_->getRecordBestSol()->update(nSol, problem_->nVars(),
00475 z, problem_->getFeasTol());
00476 best = problem_->getRecordBestSol()->getSol();
00477 objVal = problem_->getRecordBestSol()->getVal();
00478 #else
00479 best = nSol;
00480 objVal = z;
00481 #endif
00482 #endif
00483
00484 if (z < problem_ -> getCutOff ()) {
00485
00486 problem_ -> setCutOff (objVal);
00487
00488 t_chg_bounds *chg_bds = NULL;
00489
00490 if (objInd >= 0) {
00491
00492 chg_bds = new t_chg_bounds [problem_ -> nVars ()];
00493 chg_bds [objInd].setUpper (t_chg_bounds::CHANGED);
00494 }
00495
00496
00497 bool is_still_feas = problem_ -> btCore (chg_bds);
00498
00499 if (chg_bds)
00500 delete [] chg_bds;
00501
00502 if (!is_still_feas)
00503 break;
00504 }
00505 }
00506
00507 } while ((niter++ < maxIter_) &&
00508 (retval == 0));
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518 if (nlp_)
00519 nlp_ -> setObjective (originalObjective);
00520
00521 if (retval > 0) {
00522
00523 if (!nlp_)
00524 nlp_ = new CouenneTNLP (problem_);
00525
00526 problem_ -> domain () -> push (problem_ -> nVars (),
00527 problem_ -> domain () -> x (),
00528 problem_ -> domain () -> lb (),
00529 problem_ -> domain () -> ub ());
00530
00531
00532
00533
00534 fixIntVariables (best);
00535 nlp_ -> setInitSol (best);
00536
00538
00539
00540
00541
00542 status = app_ -> OptimizeTNLP (nlp_);
00543
00545
00546 problem_ -> domain () -> pop ();
00547
00548 if ((status != Solve_Succeeded) &&
00549 (status != Solved_To_Acceptable_Level))
00550
00551 problem_ -> Jnlst () -> Printf (J_ERROR, J_NLPHEURISTIC,
00552 "Feasibility Pump: error in final NLP problem (due to fixing integer variables?)\n");
00553
00554
00555
00556 double z = nlp_ -> getSolValue ();
00557
00558
00559 bool isChecked = false;
00560
00561 if (nSol) {
00562
00563 #ifdef FM_CHECKNLP2
00564 isChecked = problem_->checkNLP2(nSol, 0, false,
00565 true,
00566 true,
00567 problem_->getFeasTol());
00568 if (isChecked) {
00569 z = problem_->getRecordBestSol()->getModSolVal();
00570 }
00571
00572 #else
00573 isChecked = problem_ -> checkNLP (nSol, z, true);
00574 #endif
00575 }
00576
00577 if (nSol &&
00578 isChecked &&
00579 (z < problem_ -> getCutOff ())) {
00580
00581 #ifdef FM_CHECKNLP2
00582 #ifdef FM_TRACE_OPTSOL
00583 problem_->getRecordBestSol()->update();
00584 best = problem_->getRecordBestSol()->getSol();
00585 objVal = problem_->getRecordBestSol()->getVal();
00586 #else
00587 best = problem_->getRecordBestSol()->getModSol(problem_ -> nVars ());
00588 objVal = z;
00589 #endif
00590 #else
00591 #ifdef FM_TRACE_OPTSOL
00592 problem_->getRecordBestSol()->update(nSol, problem_->nVars(),
00593 z, problem_->getFeasTol());
00594 best = problem_->getRecordBestSol()->getSol();
00595 objVal = problem_->getRecordBestSol()->getVal();
00596 #else
00597 best = nSol;
00598 objVal = z;
00599 #endif
00600 #endif
00601
00602 problem_ -> setCutOff (objVal);
00603 }
00604 }
00605
00606 if (retval > 0) {
00607
00608 if (problem_ -> Jnlst () -> ProduceOutput (J_ERROR, J_NLPHEURISTIC)) {
00609
00610 printf ("FP: returning MINLP feasible solution:\n");
00611 printDist (problem_, best, nSol);
00612 }
00613
00614 CoinCopyN (best, problem_ -> nVars (), newSolution);
00615 }
00616
00617 if (iSol) delete [] iSol;
00618 if (nSol) delete [] nSol;
00619
00620
00621 problem_ -> domain () -> pop ();
00622
00623
00624
00625
00626
00627 delete milp_;
00628 delete postlp_;
00629 milp_ = postlp_ = NULL;
00630
00631 problem_ -> Jnlst () -> Printf
00632 (J_ERROR, J_NLPHEURISTIC, "FP: done ===================\n");
00633
00634 return retval;
00635 }
00636
00637
00638 void compDistSingle (CouenneProblem *p,
00639 int n, double *v,
00640 double &norm,
00641 int &nInfI,
00642 int &nInfN,
00643 double &infI,
00644 double &infN) {
00645
00646 p -> domain () -> push (n, v, NULL, NULL);
00647
00648 norm = infI = infN = 0.;
00649 nInfI = nInfN = 0;
00650
00651 while (n--) {
00652 norm += (*v * *v);
00653 ++v;
00654 }
00655 v -= n;
00656
00657 norm = sqrt (norm);
00658
00659 for (std::vector <exprVar *>::iterator i = p -> Variables (). begin ();
00660 i != p -> Variables (). end ();
00661 ++i) {
00662
00663 CouNumber
00664 vval = (**i) ();
00665
00666 if ((*i) -> Multiplicity () <= 0)
00667 continue;
00668
00669 if ((*i) -> isInteger ()) {
00670
00671 double inf = CoinMax (vval - floor (vval + COUENNE_EPS),
00672 ceil (vval - COUENNE_EPS) - vval);
00673
00674 if (inf > COUENNE_EPS) {
00675 ++nInfI;
00676 infI += inf;
00677 }
00678 }
00679
00680 if ((*i) -> Type () == AUX) {
00681
00682 double
00683 diff = 0.,
00684 fval = (*((*i) -> Image ())) ();
00685
00686 if ((*i) -> sign () != expression::AUX_GEQ) diff = CoinMax (diff, vval - fval);
00687 else if ((*i) -> sign () != expression::AUX_LEQ) diff = CoinMax (diff, fval - vval);
00688
00689 if (diff > COUENNE_EPS) {
00690 ++nInfN;
00691 infN += diff;
00692 }
00693 }
00694 }
00695
00696 p -> domain () -> pop ();
00697 }
00698
00699
00700
00701 void printDist (CouenneProblem *p, double *iSol, double *nSol) {
00702
00703 int nInfII = -1, nInfNI = -1, nInfIN = -1, nInfNN = -1;
00704
00705 double
00706 dist = -1.,
00707 normI = -1., normN = -1.,
00708 infII = -1., infNI = -1.,
00709 infIN = -1., infNN = -1.;
00710
00711 if (iSol) compDistSingle (p, p -> nVars (), iSol, normI, nInfII, nInfNI, infII, infNI);
00712 if (nSol) compDistSingle (p, p -> nVars (), nSol, normN, nInfIN, nInfNN, infIN, infNN);
00713
00714 if (iSol && nSol) {
00715
00716 dist = 0.;
00717
00718 for (int i = p -> nVars (); i--;)
00719 dist +=
00720 (iSol [i] - nSol [i]) *
00721 (iSol [i] - nSol [i]);
00722
00723 dist = sqrt (dist);
00724 }
00725
00726 printf ("FP: ");
00727
00728 printf ("IP norm i:%e n:%e dist %e #inf i:%4d n:%4d max inf i:%e n:%e ",
00729 normI, normN, dist, nInfII, nInfNI, infII, infNI);
00730
00731 printf ("NL #inf i:%4d n:%4d max inf i:%e n:%e\n",
00732 nInfIN, nInfNN, infIN, infNN);
00733 }
00734
00735
00736 #define WRAP 3
00737
00738 void printCmpSol (int n, double *iSol, double *nSol, int direction) {
00739
00740 printf ("i:%p n:%p\nFP # ",
00741 (void *) iSol, (void *) nSol);
00742
00743 double
00744 distance = 0.,
00745 diff;
00746
00747 char c =
00748 direction < 0 ? '<' :
00749 direction > 0 ? '>' : '-';
00750
00751 for (int i=0; i<n; i++) {
00752
00753 if (i && !(i % WRAP))
00754 printf ("\nFP # ");
00755
00756 double
00757 iS = iSol ? iSol [i] : 12345.,
00758 nS = nSol ? nSol [i] : 54321.;
00759
00760 printf ("[%4d %+e -%c- %+e (%e)] ",
00761 i, iS, c, nS,
00762 (iSol && nSol) ? fabs (iS - nS) : 0.);
00763
00764 if (iSol && nSol) {
00765 diff = iS - nS;
00766 distance += (diff*diff);
00767 }
00768 }
00769
00770 if (iSol && nSol) {
00771
00772 distance = sqrt (distance);
00773 printf ("\n### distance: %e\n", distance);
00774 }
00775 }