/home/coin/SVN-release/OS-2.4.0/Couenne/src/heuristics/CouenneFPFindSolution.cpp

Go to the documentation of this file.
00001 /* $Id: CouenneFPFindSolution.cpp 733 2011-07-03 22:32:52Z pbelotti $
00002  *
00003  * Name:    CouenneFPFindSolution.cpp
00004  * Authors: Pietro Belotti
00005  *          Timo Berthold, ZIB Berlin
00006  * Purpose: Find solution by looping through MILP solvers/heuristics
00007  * 
00008  * This file is licensed under the Eclipse Public License (EPL)
00009  */
00010 
00011 #include "CoinTime.hpp"
00012 
00013 #include "CouenneFeasPump.hpp"
00014 #include "CouenneFPpool.hpp"
00015 #include "CouenneProblem.hpp"
00016 #include "cons_rowcuts.h"
00017 
00018 #ifdef COIN_HAS_SCIP
00019 /* general SCIP includes */
00020 #include "scip/scip.h"
00021 #include "scip/cons_linear.h"
00022 #include "scip/scipdefplugins.h"
00023 #endif
00024 
00025 using namespace Couenne;
00026 
00028 double CouenneFeasPump::findSolution (double* &sol, int niter, int* nsuciter) {
00029 
00040 
00041   // What order should we use? I suggest we use priorities, assigned
00042   // at the beginning but changeable in the event of multiple failures
00043   // (or successes) of a given method.
00044   //
00045   // Rule of thumb: 
00046   //
00047   // 1) Assign all methods i a number p[i] (for instance those in the
00048   //    list above)
00049   //
00050   // 2) Call each in the order defined by p[i], return a solution if
00051   //    found, otherwise proceed to next method
00052   //
00053   // 3) If K consecutive successes at finding new solution (not
00054   //    necessarily new best feasible), --p[i]
00055   //
00056   // 4) if H consecutive failutes, ++p[i]
00057 
00058   double obj;
00059 
00061 
00062 #ifdef COIN_HAS_SCIP
00063 
00064   static int currentmilpmethod = 0;
00065 
00066   int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
00067 
00068   if (useSCIP_ && problem_ -> nIntVars () > 0) { // if LP, use Clp below
00069 
00070      SCIP* scip;
00071 
00072      SCIP_VAR** vars;
00073      const SCIP_Real* lbs;
00074      const SCIP_Real* ubs;
00075      const SCIP_Real* objs;
00076      const char* vartypes;
00077      const CoinPackedMatrix * matrix;
00078      const CoinBigIndex* rowstarts;
00079      const int* rowlengths;
00080      const SCIP_Real* coeffs;
00081      const SCIP_Real* lhss;
00082      const SCIP_Real* rhss;
00083      const int* indices;
00084 
00085      SCIP_Real timelimit;
00086 
00087      double infinity;
00088      int nvars;
00089      int nconss;
00090      int nscipsols;
00091 
00092      bool solveagain;
00093 
00094      // COUENNE_INFINITY , getInfinity()
00095 
00096      // get problem data
00097      nvars    = milp_ -> getNumCols ();
00098      nconss   = milp_ -> getNumRows ();
00099      infinity = milp_ -> getInfinity ();
00100 
00101      // get variable data
00102      lbs =      milp_ -> getColLower ();
00103      ubs =      milp_ -> getColUpper ();
00104      objs =     milp_ -> getObjCoefficients ();
00105      vartypes = milp_ -> getColType ();
00106 
00107      // get row data
00108      lhss = milp_ -> getRowLower ();
00109      rhss = milp_ -> getRowUpper ();
00110 
00111      // get matrix data
00112      matrix     = milp_ -> getMatrixByRow();
00113      rowstarts  = matrix -> getVectorStarts();
00114      rowlengths = matrix -> getVectorLengths();
00115      coeffs     = matrix -> getElements();
00116      indices    = matrix -> getIndices();
00117      
00118      if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
00119         SCIPdebugMessage("create SCIP problem instance with %d variables and %d constraints.\n", nvars, nconss);
00120      }
00121 
00122      // initialize SCIP
00123      SCIP_CALL_ABORT( SCIPcreate(&scip) );
00124      assert(scip != NULL);
00125    
00126      // include default SCIP plugins
00127      SCIP_CALL_ABORT( SCIPincludeDefaultPlugins(scip) );
00128 
00129      // include row cut constraint hanlder
00130      if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
00131      { 
00132         SCIP_CALL_ABORT( SCIPincludeConshdlrRowcuts(scip, couenneCG_, milp_) );
00133      }
00134 
00135      // create problem instance in SCIP
00136      SCIP_CALL_ABORT( SCIPcreateProb(scip, "auxiliary FeasPump MILP", NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
00137 
00138      // allocate local memory for SCIP variable array
00139      SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &vars, nvars) );
00140    
00141      // one variable for objective !!!!!!!!!
00142 
00143      // create variables 
00144      for (int i=0; i<nvars; i++) {
00145         char varname[SCIP_MAXSTRLEN];  
00146 
00147         // check that all data is in valid ranges
00148         assert( 0 <= vartypes[i] && vartypes[i] <= 2);
00149         checkInfinity(scip, lbs[i], infinity);
00150         checkInfinity(scip, ubs[i], infinity);
00151         
00152         // all variables are named x_i
00153         (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "x_%d", i);
00154         SCIP_CALL_ABORT( SCIPcreateVar(scip, &vars[i], varname, lbs[i], ubs[i], objs[i], 
00155               vartypes[i] == 0 ? SCIP_VARTYPE_CONTINUOUS : (vartypes[i] == 1 ? SCIP_VARTYPE_BINARY : SCIP_VARTYPE_INTEGER),
00156               TRUE, FALSE, NULL, NULL, NULL, NULL, NULL) );
00157 
00158         // add the variable to SCIP
00159         SCIP_CALL_ABORT( SCIPaddVar(scip, vars[i]) );
00160      }
00161 
00162      // create constraints
00163      for (int i=0; i<nconss; i++) {
00164 
00165         SCIP_CONS* cons;
00166         
00167         char consname[SCIP_MAXSTRLEN];  
00168         (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN, "row_%d", i);
00169 
00170         // check that all data is in valid ranges
00171         checkInfinity(scip, lhss[i], infinity);
00172         checkInfinity(scip, rhss[i], infinity);
00173 
00174         // create an empty linear constraint
00175         SCIP_CALL_ABORT( SCIPcreateConsLinear(scip, &cons, consname, 0, NULL, NULL, lhss[i], rhss[i], 
00176               TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
00177              
00178         // add variables to constraint
00179         for(int j=rowstarts[i]; j<rowstarts[i]+rowlengths[i]; j++)        
00180         {
00181            checkInfinity(scip, coeffs[j], infinity);
00182            SCIP_CALL_ABORT( SCIPaddCoefLinear(scip, cons, vars[indices[j]], coeffs[j]) );
00183         }
00184 
00185         // add constraint to SCIP
00186         SCIP_CALL_ABORT( SCIPaddCons(scip, cons) );
00187         SCIP_CALL_ABORT( SCIPreleaseCons(scip, &cons) );        
00188      }
00189   
00190      // determine the method to solve the MILP
00191      if (milpMethod_ == 0 && niter == 0)
00192      {
00193         // initialize currentmilpmethod: at the root node we solve the MILP with SCIP default
00194         // deeper in the tree, we will solve the MILP by a FP heuristic
00195         if( depth == 0 )
00196            currentmilpmethod = 2;
00197         else 
00198            currentmilpmethod = 4;
00199      }
00200      else if (milpMethod_ != 0)
00201         currentmilpmethod = milpMethod_; // use a fixed method to solve the MILP
00202      
00203      
00204      // MILP solving loop. If the MILP terminates without a solution, it might get resolved with a more expensive atrategy
00205      do {
00206         solveagain = false;
00207         
00208         // reset parameters if MILP is solved agian
00209         SCIP_CALL( SCIPresetParams(scip) );
00210         
00211         // deactivate SCIP output
00212         if (!(problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC))) {
00213            SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 0) );
00214         }
00215         
00216         // do not abort subproblem on CTRL-C
00217         SCIP_CALL( SCIPsetBoolParam(scip, "misc/catchctrlc", FALSE) );
00218         
00219         // set time limit
00220         timelimit = problem_ -> getMaxCpuTime () - CoinCpuTime ();
00221         SCIP_CALL( SCIPsetRealParam(scip, "limits/time", timelimit) );        
00222 
00223 
00224         if (problem_ -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_NLPHEURISTIC)) {
00225            SCIPinfoMessage(scip, NULL, "using MILP method: %d\n",currentmilpmethod);
00226         }
00227         
00228         // tune SCIP differently, depending on the chosen method to solve the MILP
00235         switch(currentmilpmethod)
00236         {
00237         case -1: // solve the MILP completely. SCIP's default setting should be best for this
00238            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
00239            { 
00240               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
00241            }
00242            break;
00243 
00244         case 1: // Be aggressive in finding feasible solutions, but lazy about the dual bound. 
00245            // Set limits on overall nodes and stall nodes (nodes without incumbent improvement).
00246            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/stallnodes", 1000) );
00247            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 10000) );
00248 
00249            // disable cutting plane separation 
00250            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
00251         
00252            // disable expensive presolving 
00253            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
00254 
00255            // use aggressive primal heuristics 
00256            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_AGGRESSIVE, TRUE) );
00257         
00258            // use best estimate node selection 
00259            if( SCIPfindNodesel(scip, "estimate") != NULL )
00260            {
00261               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "nodeselection/estimate/stdpriority", INT_MAX/4) ); 
00262            }
00263         
00264            // use inference branching 
00265            if( SCIPfindBranchrule(scip, "inference") != NULL )
00266            {
00267               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
00268            }
00269         
00270            // disable conflict analysis 
00271            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useprop", FALSE) );
00272            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useinflp", FALSE) );
00273            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/useboundlp", FALSE) );
00274            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/usesb", FALSE) );
00275            SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "conflict/usepseudo", FALSE) );
00276 
00277            break;
00278 
00279         case 2: // default SCIP with node limits
00280            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/stallnodes", 500) );
00281            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 5000) );
00282            break;
00283         case 3: // solve the MILP with RENS. Disable most other features, enable RENS
00284            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 1) );
00285 
00286            // disable cutting plane separation 
00287            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
00288         
00289            // disable expensive presolving 
00290            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
00291 
00292            // besides RENS, only use cheap heuristics 
00293            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
00294 
00295            // use inference branching 
00296            if( SCIPfindBranchrule(scip, "inference") != NULL )
00297            {
00298               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
00299            }
00300 
00301            // ensure that RENS is called
00302            if( SCIPfindHeur(scip, "rens") != NULL )
00303            {
00304               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/rens/freq", 0) );
00305               SCIP_CALL_ABORT( SCIPsetRealParam(scip, "heuristics/rens/minfixingrate", 0.0) );
00306            }
00307            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
00308            { 
00309               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
00310            }
00311            break;
00312 
00313         case 4: // solve the MILP with Feasibility Pump. Disable most other features, enable stage 3 for feaspump
00314            SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "limits/nodes", 1) );
00315 
00316            // disable cutting plane separation 
00317            SCIP_CALL_ABORT( SCIPsetSeparating(scip, SCIP_PARAMSETTING_OFF, TRUE) );
00318         
00319            // disable expensive presolving 
00320            SCIP_CALL_ABORT( SCIPsetPresolving(scip, SCIP_PARAMSETTING_FAST, TRUE) );
00321 
00322            // besides feaspump, only use cheap heuristics 
00323            SCIP_CALL_ABORT( SCIPsetHeuristics(scip, SCIP_PARAMSETTING_FAST, TRUE) );
00324 
00325            // use inference branching 
00326            if( SCIPfindBranchrule(scip, "inference") != NULL )
00327            {
00328               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "branching/inference/priority", INT_MAX/4) );
00329            }
00330 
00331            // ensure that feasibility pump is called
00332            if( SCIPfindHeur(scip, "feaspump") != NULL )
00333            {
00334               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/feaspump/freq", 0) );
00335               SCIP_CALL_ABORT( SCIPsetIntParam(scip, "heuristics/feaspump/maxsols", -1) );
00336 
00337               if( SCIPversion() <= 2.01 )
00338               {
00339                  SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "heuristics/feaspump2/stage3", TRUE) );
00340               }
00341               else
00342               {
00343                  SCIP_CALL_ABORT( SCIPsetBoolParam(scip, "heuristics/feaspump/stage3", TRUE) );
00344               }
00345            }
00346            if( milpCuttingPlane_ == FP_CUT_INTEGRATED )
00347            { 
00348               SCIP_CALL_ABORT( SCIPsetLongintParam(scip, "constraints/rowcuts/maxcuttingrounds", 0) );
00349            }
00350            break;
00351 
00352         default:
00353            printf("invalid MILP method: %d\n", currentmilpmethod);
00354            assert(false);
00355            break;
00356         }
00357 
00358 #if 0
00359      // writes MILP problem and SCIP settings into a file
00360      SCIP_CALL_ABORT( SCIPwriteOrigProblem(scip, "debug.lp", NULL, FALSE) );
00361      SCIP_CALL_ABORT( SCIPwriteParams(scip, "debug.set", FALSE,TRUE) );
00362 #endif
00363           
00364         // solve the MILP
00365         SCIP_CALL_ABORT( SCIPsolve(scip) );
00366 
00367         nscipsols =  SCIPgetNSols(scip);
00368      
00369         // copy the solution
00370         if( nscipsols)
00371         {
00372            SCIP_SOL** scipsols;
00373            SCIP_SOL* bestsol;
00374            SCIP_Real cutoffbound;
00375 
00376            int nstoredsols;
00377 
00378            /* get incumbent solution */
00379            bestsol = SCIPgetBestSol(scip);
00380            assert(bestsol != NULL);
00381 
00382            /* get SCIP solution pool */
00383            scipsols = SCIPgetSols(scip);
00384            assert(scipsols != NULL);
00385 
00386            if (!sol)
00387              sol = new CouNumber [nvars];
00388 
00389            // get solution values and objective of incumbent
00390            SCIP_CALL_ABORT( SCIPgetSolVals(scip, bestsol, nvars, vars, sol) );
00391            obj = SCIPgetSolOrigObj(scip, bestsol);
00392 
00393            nstoredsols = 0;
00394 
00395            // if we use an objective feaspump, the objective function might be negative
00396            cutoffbound = obj > 0 ? 2*obj : obj / 2.0;
00397 
00398            // insert other SCIP solutions into solution pool
00399            // do not store too many or too poor solutions 
00400            for(int i=1; i<nscipsols && nstoredsols < 10 && 
00401                   SCIPgetSolOrigObj(scip,scipsols[i]) <= cutoffbound; i++){
00402               double* tmpsol;
00403 
00404               tmpsol = new CouNumber [nvars];
00405            
00406               // get solution values
00407               SCIP_CALL_ABORT( SCIPgetSolVals(scip, scipsols[i], nvars, vars, tmpsol) );
00408               CouenneFPsolution couennesol = CouenneFPsolution (problem_, tmpsol);
00409 
00410               // add solutions to the pool if they are not in the tabu list
00411               if (tabuPool_. find (couennesol) == tabuPool_ . end () 
00412                   && pool_ -> Set (). find (couennesol) == pool_ -> Set(). end ()
00413                  ){
00414                  pool_ -> Set (). insert (couennesol);
00415 
00416                  ++nstoredsols;
00417               }
00418            }
00419 
00420            ++(*nsuciter);
00421 
00422            // if we succeeded five times in a row, try a cheaper MILP_ solving method next time
00423            // TODO: if we want to use time limits, hitting the time limit would be another good reason to switch
00424            if( *nsuciter >= 3 && currentmilpmethod < 4 )
00425            {
00426               ++currentmilpmethod;
00427               *nsuciter = 0;
00428            }          
00429         }
00430         //try to use a more aggressive, more expensive way to solve the MILP
00431         else if( milpMethod_ == 0 && currentmilpmethod > 1 )
00432         {
00433            --currentmilpmethod;
00434            solveagain = true;
00435            *nsuciter = 0;
00436 
00437            // throw away the current solution process
00438            SCIP_CALL( SCIPfreeTransform(scip) );
00439         }
00440         else
00441            obj = COIN_DBL_MAX;  
00442 
00443      } while (solveagain);
00444    
00445      // release variables before freeing them
00446      for (int i=0; i<nvars; i++) {
00447         SCIP_CALL_ABORT( SCIPreleaseVar(scip, &vars[i]) );
00448      }
00449 
00450      // free memory
00451      SCIPfreeMemoryArray(scip, &vars);
00452      SCIP_CALL_ABORT( SCIPfree(&scip) );
00453    
00454      BMScheckEmptyMemory();     
00455   }
00456   else
00457 #endif      
00458   {
00459 
00460      if (problem_ -> nIntVars () > 0) milp_ -> branchAndBound ();
00461      else                             milp_ -> initialSolve ();
00462 
00463      if (!sol)
00464        sol = new CouNumber [problem_ -> nVars ()];
00465 
00466      if (milp_ -> getColSolution ())
00467        CoinCopyN (milp_ -> getColSolution (), problem_ -> nVars (), sol);
00468      else {
00469 
00470        if (sol)
00471          delete [] sol;
00472        sol = NULL;
00473      }
00474 
00475      obj = milp_ -> getObjValue ();
00476   }
00477 
00478   return obj;
00479 }
00480 
00482 void CouenneFeasPump::init_MILP () {
00483 
00484 }

Generated on Thu Sep 22 03:05:58 2011 by  doxygen 1.4.7