1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file heur_feaspump.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief Objective Feasibility Pump 2.0 28 * @author Timo Berthold 29 * @author Domenico Salvagnin 30 */ 31 32 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 33 34 #include "blockmemshell/memory.h" 35 #include "scip/cons_linear.h" 36 #include "scip/heur_feaspump.h" 37 #include "scip/pub_heur.h" 38 #include "scip/pub_message.h" 39 #include "scip/pub_misc.h" 40 #include "scip/pub_misc_sort.h" 41 #include "scip/pub_var.h" 42 #include "scip/scip_branch.h" 43 #include "scip/scip_cons.h" 44 #include "scip/scip_copy.h" 45 #include "scip/scip_general.h" 46 #include "scip/scip_heur.h" 47 #include "scip/scip_lp.h" 48 #include "scip/scip_mem.h" 49 #include "scip/scip_message.h" 50 #include "scip/scip_nodesel.h" 51 #include "scip/scip_numerics.h" 52 #include "scip/scip_param.h" 53 #include "scip/scip_pricer.h" 54 #include "scip/scip_prob.h" 55 #include "scip/scip_probing.h" 56 #include "scip/scip_randnumgen.h" 57 #include "scip/scip_sol.h" 58 #include "scip/scip_solve.h" 59 #include "scip/scip_solvingstats.h" 60 #include "scip/scip_tree.h" 61 #include "scip/scip_var.h" 62 #include <string.h> 63 64 #define HEUR_NAME "feaspump" 65 #define HEUR_DESC "objective feasibility pump 2.0" 66 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_OBJDIVING 67 #define HEUR_PRIORITY -1000000 68 #define HEUR_FREQ 20 69 #define HEUR_FREQOFS 0 70 #define HEUR_MAXDEPTH -1 71 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPPLUNGE 72 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */ 73 74 #define DEFAULT_MAXLPITERQUOT 0.01 /**< maximal fraction of diving LP iterations compared to node LP iterations */ 75 #define DEFAULT_MAXLPITEROFS 1000 /**< additional number of allowed LP iterations */ 76 #define DEFAULT_MAXSOLS 10 /**< total number of feasible solutions found up to which heuristic is called 77 * (-1: no limit) */ 78 #define DEFAULT_MAXLOOPS 10000 /**< maximal number of pumping rounds (-1: no limit) */ 79 #define DEFAULT_MAXSTALLLOOPS 10 /**< maximal number of pumping rounds without fractionality improvement (-1: no limit) */ 80 #define DEFAULT_MINFLIPS 10 /**< minimum number of random variables to flip, if a 1-cycle is encountered */ 81 #define DEFAULT_CYCLELENGTH 3 /**< maximum length of cycles to be checked explicitly in each round */ 82 #define DEFAULT_PERTURBFREQ 100 /**< number of iterations until a random perturbation is forced */ 83 #define DEFAULT_OBJFACTOR 0.1 /**< factor by which the regard of the objective is decreased in each round, 84 * 1.0 for dynamic, depending on solutions already found */ 85 #define DEFAULT_ALPHA 1.0 /**< initial weight of the objective function in the convex combination */ 86 #define DEFAULT_ALPHADIFF 1.0 /**< threshold difference for the convex parameter to perform perturbation */ 87 #define DEFAULT_BEFORECUTS TRUE /**< should the feasibility pump be called at root node before cut separation? */ 88 #define DEFAULT_USEFP20 FALSE /**< should an iterative round-and-propagate scheme be used to find the integral points? */ 89 #define DEFAULT_PERTSOLFOUND TRUE /**< should a random perturbation be performed if a feasible solution was found? */ 90 #define DEFAULT_STAGE3 FALSE /**< should we solve a local branching sub-MIP if no solution could be found? */ 91 #define DEFAULT_NEIGHBORHOODSIZE 18 /**< radius of the neighborhood to be searched in stage 3 */ 92 #define DEFAULT_COPYCUTS TRUE /**< should all active cuts from the cutpool of the original SCIP be copied to 93 * constraints of the subscip 94 */ 95 96 #define MINLPITER 5000 /**< minimal number of LP iterations allowed in each LP solving call */ 97 98 #define DEFAULT_RANDSEED 13 /**< initial random seed */ 99 100 /** primal heuristic data */ 101 struct SCIP_HeurData 102 { 103 SCIP_SOL* sol; /**< working solution */ 104 SCIP_SOL* roundedsol; /**< rounded solution */ 105 SCIP_Longint nlpiterations; /**< number of LP iterations used in this heuristic */ 106 SCIP_Real maxlpiterquot; /**< maximal fraction of diving LP iterations compared to node LP iterations */ 107 SCIP_Real objfactor; /**< factor by which the regard of the objective is decreased in each round, 108 * 1.0 for dynamic, depending on solutions already found */ 109 SCIP_Real alpha; /**< initial weight of the objective function in the convex combination */ 110 SCIP_Real alphadiff; /**< threshold difference for the convex parameter to perform perturbation */ 111 112 int maxlpiterofs; /**< additional number of allowed LP iterations */ 113 int maxsols; /**< total number of feasible solutions found up to which heuristic is called 114 * (-1: no limit) */ 115 int maxloops; /**< maximum number of loops (-1: no limit) */ 116 int maxstallloops; /**< maximal number of pumping rounds without fractionality improvement (-1: no limit) */ 117 int minflips; /**< minimum number of random variables to flip, if a 1-cycle is encountered */ 118 int cyclelength; /**< maximum length of cycles to be checked explicitly in each round */ 119 int perturbfreq; /**< number of iterations until a random perturbation is forced */ 120 int nsuccess; /**< number of runs that produced at least one feasible solution */ 121 int neighborhoodsize; /**< radius of the neighborhood to be searched in stage 3 */ 122 123 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */ 124 SCIP_Bool beforecuts; /**< should the feasibility pump be called at root node before cut separation? */ 125 SCIP_Bool usefp20; /**< should an iterative round-and-propagate scheme be used to find the integral points? */ 126 SCIP_Bool pertsolfound; /**< should a random perturbation be performed if a feasible solution was found? */ 127 SCIP_Bool stage3; /**< should we solve a local branching sub-MIP if no solution could be found? */ 128 SCIP_Bool copycuts; /**< should all active cuts from cutpool be copied to constraints in 129 * subproblem? 130 */ 131 }; 132 133 /* copies SCIP to probing SCIP and creates variable hashmap */ 134 static 135 SCIP_RETCODE setupProbingSCIP( 136 SCIP* scip, /**< SCIP data structure */ 137 SCIP** probingscip, /**< sub-SCIP data structure */ 138 SCIP_HASHMAP** varmapfw, /**< mapping of SCIP variables to sub-SCIP variables */ 139 SCIP_Bool copycuts, /**< should all active cuts from cutpool of scip copied to constraints in subscip */ 140 SCIP_Bool* success /**< was copying successful? */ 141 ) 142 { 143 /* check if we are already at the maximal tree depth */ 144 if( SCIP_MAXTREEDEPTH <= SCIPgetDepth(scip) ) 145 { 146 *success = FALSE; 147 return SCIP_OKAY; 148 } 149 150 /* initializing the subproblem */ 151 SCIP_CALL( SCIPcreate(probingscip) ); 152 153 /* create the variable mapping hash map */ 154 SCIP_CALL( SCIPhashmapCreate(varmapfw, SCIPblkmem(*probingscip), SCIPgetNVars(scip)) ); 155 *success = FALSE; 156 157 /* copy SCIP instance */ 158 SCIP_CALL( SCIPcopyConsCompression(scip, *probingscip, *varmapfw, NULL, "feaspump", NULL, NULL, 0, FALSE, FALSE, 159 FALSE, TRUE, success) ); 160 161 if( copycuts ) 162 { 163 /* copies all active cuts from cutpool of sourcescip to linear constraints in targetscip */ 164 SCIP_CALL( SCIPcopyCuts(scip, *probingscip, *varmapfw, NULL, FALSE, NULL) ); 165 } 166 167 return SCIP_OKAY; 168 } 169 170 /** set appropriate parameters for probing SCIP in FP2 */ 171 static 172 SCIP_RETCODE setupSCIPparamsFP2( 173 SCIP* scip, /**< SCIP data structure */ 174 SCIP* probingscip /**< sub-SCIP data structure */ 175 ) 176 { 177 if( SCIPisParamFixed(probingscip, "heuristics/" HEUR_NAME "/freq") ) 178 { 179 SCIPwarningMessage(scip, "unfixing parameter heuristics/" HEUR_NAME "/freq in probingscip of " HEUR_NAME " heuristic to avoid recursive calls\n"); 180 SCIP_CALL( SCIPunfixParam(probingscip, "heuristics/" HEUR_NAME "/freq") ); 181 } 182 SCIP_CALL( SCIPsetIntParam(probingscip, "heuristics/" HEUR_NAME "/freq", -1) ); 183 184 /* do not abort subproblem on CTRL-C */ 185 SCIP_CALL( SCIPsetBoolParam(probingscip, "misc/catchctrlc", FALSE) ); 186 187 #ifndef SCIP_DEBUG 188 /* disable output to console */ 189 SCIP_CALL( SCIPsetIntParam(probingscip, "display/verblevel", 0) ); 190 #endif 191 192 /* do not multiaggregate variables, because otherwise they have to be skipped in the fix-and-propagate loop */ 193 SCIP_CALL( SCIPsetBoolParam(probingscip, "presolving/donotmultaggr", TRUE) ); 194 195 /* limit to root node solving */ 196 SCIP_CALL( SCIPsetLongintParam(probingscip, "limits/nodes", 1LL) ); 197 198 /* disable LP solving and expensive techniques */ 199 if( SCIPisParamFixed(probingscip, "lp/solvefreq") ) 200 { 201 SCIPwarningMessage(scip, "unfixing parameter lp/solvefreq in probingscip of " HEUR_NAME " heuristic to speed up propagation\n"); 202 SCIP_CALL( SCIPunfixParam(probingscip, "lp/solvefreq") ); 203 } 204 SCIP_CALL( SCIPsetIntParam(probingscip, "lp/solvefreq", -1) ); 205 SCIP_CALL( SCIPsetBoolParam(probingscip, "conflict/enable", FALSE) ); 206 SCIP_CALL( SCIPsetBoolParam(probingscip, "constraints/disableenfops", TRUE) ); 207 SCIP_CALL( SCIPsetBoolParam(probingscip, "constraints/knapsack/negatedclique", FALSE) ); 208 SCIP_CALL( SCIPsetPresolving(probingscip, SCIP_PARAMSETTING_FAST, TRUE) ); 209 SCIP_CALL( SCIPsetHeuristics(probingscip, SCIP_PARAMSETTING_OFF, TRUE) ); 210 211 return SCIP_OKAY; 212 } 213 214 /** set appropriate parameters for probing SCIP in Stage 3 */ 215 static 216 SCIP_RETCODE setupSCIPparamsStage3( 217 SCIP* scip, /**< SCIP data structure */ 218 SCIP* probingscip /**< sub-SCIP data structure */ 219 ) 220 { 221 SCIP_CALL( SCIPcopyParamSettings(scip, probingscip) ); 222 /* do not abort subproblem on CTRL-C */ 223 SCIP_CALL( SCIPsetBoolParam(probingscip, "misc/catchctrlc", FALSE) ); 224 225 #ifndef SCIP_DEBUG 226 /* disable output to console */ 227 SCIP_CALL( SCIPsetIntParam(probingscip, "display/verblevel", 0) ); 228 #endif 229 /* set limits for the subproblem */ 230 SCIP_CALL( SCIPcopyLimits(scip, probingscip) ); 231 SCIP_CALL( SCIPsetLongintParam(probingscip, "limits/nodes", 1000LL) ); 232 SCIP_CALL( SCIPsetLongintParam(probingscip, "limits/stallnodes", 100LL) ); 233 234 /* forbid recursive call of heuristics and separators solving sub-SCIPs */ 235 SCIP_CALL( SCIPsetSubscipsOff(probingscip, TRUE) ); 236 if( SCIPisParamFixed(probingscip, "heuristics/" HEUR_NAME "/freq") ) 237 { 238 SCIPwarningMessage(scip,"unfixing parameter heuristics/" HEUR_NAME "/freq in probingscip of " HEUR_NAME " heuristic to avoid recursive calls\n"); 239 SCIP_CALL( SCIPunfixParam(probingscip, "heuristics/" HEUR_NAME "/freq") ); 240 } 241 SCIP_CALL( SCIPsetIntParam(probingscip, "heuristics/feaspump/freq", -1) ); 242 243 /* disable heuristics which aim to feasibility instead of optimality */ 244 if( !SCIPisParamFixed(probingscip, "heuristics/octane/freq") ) 245 { 246 SCIP_CALL( SCIPsetIntParam(probingscip, "heuristics/octane/freq", -1) ); 247 } 248 if( !SCIPisParamFixed(probingscip, "heuristics/objpscostdiving/freq") ) 249 { 250 SCIP_CALL( SCIPsetIntParam(probingscip, "heuristics/objpscostdiving/freq", -1) ); 251 } 252 if( !SCIPisParamFixed(probingscip, "heuristics/rootsoldiving/freq") ) 253 { 254 SCIP_CALL( SCIPsetIntParam(probingscip, "heuristics/rootsoldiving/freq", -1) ); 255 } 256 257 /* disable cutting plane separation */ 258 SCIP_CALL( SCIPsetSeparating(probingscip, SCIP_PARAMSETTING_OFF, TRUE) ); 259 260 /* disable expensive presolving */ 261 SCIP_CALL( SCIPsetPresolving(probingscip, SCIP_PARAMSETTING_FAST, TRUE) ); 262 263 /* use best estimate node selection */ 264 if( SCIPfindNodesel(probingscip, "estimate") != NULL && !SCIPisParamFixed(probingscip, "nodeselection/estimate/stdpriority") ) 265 { 266 SCIP_CALL( SCIPsetIntParam(probingscip, "nodeselection/estimate/stdpriority", INT_MAX/4) ); 267 } 268 269 /* use inference branching */ 270 if( SCIPfindBranchrule(probingscip, "inference") != NULL && !SCIPisParamFixed(probingscip, "branching/inference/priority") ) 271 { 272 SCIP_CALL( SCIPsetIntParam(probingscip, "branching/inference/priority", INT_MAX/4) ); 273 } 274 275 /* disable conflict analysis */ 276 if( !SCIPisParamFixed(probingscip, "conflict/enable") ) 277 { 278 SCIP_CALL( SCIPsetBoolParam(probingscip, "conflict/enable", FALSE) ); 279 } 280 281 return SCIP_OKAY; 282 } 283 284 /** checks whether a variable is one of the currently most fractional ones */ 285 static 286 void insertFlipCand( 287 SCIP_VAR** mostfracvars, /**< sorted array of the currently most fractional variables */ 288 SCIP_Real* mostfracvals, /**< array of their fractionality, decreasingly sorted */ 289 int* nflipcands, /**< number of fractional variables already labeled to be flipped*/ 290 int maxnflipcands, /**< typically randomized number of maximum amount of variables to flip */ 291 SCIP_VAR* var, /**< variable to be checked */ 292 SCIP_Real frac /**< fractional value of the variable */ 293 ) 294 { 295 int i; 296 297 assert(mostfracvars != NULL); 298 assert(mostfracvals != NULL); 299 assert(nflipcands != NULL); 300 301 /* instead of the fractional value use the fractionality */ 302 if( frac > 0.5 ) 303 frac = 1 - frac; 304 305 /* if there are already enough candidates and the variable is less fractional, return, else reserve the last entry */ 306 if( *nflipcands >= maxnflipcands ) 307 { 308 if( frac <= mostfracvals[*nflipcands-1] ) 309 return; 310 else 311 (*nflipcands)--; 312 } 313 314 /* shifting var and frac through the (sorted) arrays */ 315 for( i = *nflipcands; i > 0 && mostfracvals[i-1] < frac; i-- ) 316 { 317 mostfracvars[i] = mostfracvars[i-1]; 318 mostfracvals[i] = mostfracvals[i-1]; 319 } 320 assert(0 <= i && i <= *nflipcands && *nflipcands < maxnflipcands); 321 322 /* insert the variable and its fractionality */ 323 mostfracvars[i] = var; 324 mostfracvals[i] = frac; 325 326 /* we've found another candidate */ 327 (*nflipcands)++; 328 } 329 330 /** set solution value in rounded solution and update objective coefficient accordingly */ 331 static 332 SCIP_RETCODE updateVariableRounding( 333 SCIP* scip, /**< SCIP data structure */ 334 SCIP_HEURDATA* heurdata, /**< heuristic data structure */ 335 SCIP_VAR* var, /**< variable */ 336 SCIP_Real solval, /**< solution value for this variable */ 337 SCIP_Real alpha, /**< weight of original objective function */ 338 SCIP_Real scalingfactor /**< factor to scale the original objective function with */ 339 ) 340 { 341 SCIP_Real lb; 342 SCIP_Real ub; 343 SCIP_Real newobjcoeff; 344 SCIP_Real orgobjcoeff; 345 346 assert(heurdata != NULL); 347 assert(var != NULL); 348 349 lb = SCIPvarGetLbLocal(var); 350 ub = SCIPvarGetUbLocal(var); 351 352 /* update rounded solution */ 353 assert(SCIPisFeasLE(scip, lb, solval) && SCIPisFeasLE(scip, solval, ub)); 354 assert(SCIPisIntegral(scip,solval)); 355 SCIP_CALL( SCIPsetSolVal(scip, heurdata->roundedsol, var, solval) ); 356 357 /* modify objective towards rounded solution value if it is at one of the variable bounds */ 358 orgobjcoeff = SCIPvarGetObj(var); 359 if( SCIPisEQ(scip, solval, lb) ) 360 newobjcoeff = (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff; 361 else if( SCIPisEQ(scip, solval, ub) ) 362 newobjcoeff = - (1.0 - alpha)/scalingfactor + alpha * orgobjcoeff; 363 else 364 newobjcoeff = alpha * orgobjcoeff; 365 366 SCIP_CALL( SCIPchgVarObjDive(scip, var, newobjcoeff) ); 367 368 return SCIP_OKAY; 369 } 370 371 372 /** flips the roundings of the most fractional variables, if a 1-cycle was found */ 373 static 374 SCIP_RETCODE handle1Cycle( 375 SCIP* scip, /**< SCIP data structure */ 376 SCIP_HEURDATA* heurdata, /**< data of this special heuristic */ 377 SCIP_VAR** mostfracvars, /**< sorted array of the currently most fractional variables */ 378 int nflipcands, /**< number of variables to flip */ 379 SCIP_Real alpha, /**< factor how much the original objective is regarded */ 380 SCIP_Real scalingfactor /**< factor to scale the original objective function with */ 381 ) 382 { 383 int i; 384 385 /* flip rounding to the opposite side */ 386 for( i = 0; i < nflipcands; i++ ) 387 { 388 SCIP_VAR* var; 389 SCIP_Real solval; 390 SCIP_Real roundedsolval; 391 392 var = mostfracvars[i]; 393 solval = SCIPvarGetLPSol(var); 394 roundedsolval = SCIPgetSolVal(scip, heurdata->roundedsol, var); 395 396 assert(! SCIPisFeasIntegral(scip, solval)); 397 assert(SCIPisFeasIntegral(scip, roundedsolval)); 398 399 /* flip to the opposite rounded solution value */ 400 if( roundedsolval > solval ) 401 solval = SCIPfeasFloor(scip, solval); 402 else 403 { 404 solval = SCIPfeasCeil(scip, solval); 405 } 406 407 SCIPdebugMsg(scip, "1-cycle flip: variable <%s> [%g,%g] LP sol %.15g sol %.15g -> %.15g\n", 408 SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), 409 SCIPvarGetLPSol(var), SCIPgetSolVal(scip, heurdata->roundedsol, var), solval); 410 411 SCIP_CALL( updateVariableRounding(scip, heurdata, var, solval, alpha, scalingfactor) ); 412 } 413 return SCIP_OKAY; 414 } 415 416 /** flips the roundings of randomly chosen fractional variables, preferring highly fractional ones, 417 * if a longer cycle was found 418 */ 419 static 420 SCIP_RETCODE handleCycle( 421 SCIP* scip, /**< SCIP data structure */ 422 SCIP_HEURDATA* heurdata, /**< data of this special heuristic */ 423 SCIP_VAR** vars, /**< array of all variables */ 424 int nbinandintvars, /**< number of general integer and 0-1 variables */ 425 SCIP_Real alpha, /**< factor how much the original objective is regarded */ 426 SCIP_Real scalingfactor /**< factor to scale the original objective function with */ 427 ) 428 { 429 int i; 430 431 /* flip variables randomized biased on their fractionality */ 432 for( i = 0; i < nbinandintvars; i++ ) 433 { 434 SCIP_VAR* var; 435 SCIP_Real solval; 436 SCIP_Real frac; 437 SCIP_Real flipprob; 438 SCIP_Real roundedsolval; 439 440 var = vars[i]; 441 solval = SCIPvarGetLPSol(var); 442 443 /* skip variables with integral solution values */ 444 if( SCIPisFeasIntegral(scip, solval) ) 445 continue; 446 447 frac = SCIPfeasFrac(scip, solval); 448 flipprob = SCIPrandomGetReal(heurdata->randnumgen, -0.3, 0.7); 449 450 /* flip, iff the sum of the randomized number and the fractionality is big enough */ 451 if( MIN(frac, 1.0 - frac) + MAX(flipprob, 0.0) > 0.5 ) 452 { 453 roundedsolval = SCIPgetSolVal(scip, heurdata->roundedsol, var); 454 assert(SCIPisFeasIntegral(scip, roundedsolval)); 455 456 /* flip the solution to the opposite side */ 457 if( roundedsolval > solval ) 458 solval = SCIPfloor(scip, solval); 459 else 460 solval = SCIPceil(scip, solval); 461 462 /* update rounded solution value and objective coefficient */ 463 SCIP_CALL( updateVariableRounding(scip, heurdata, var, solval, alpha, scalingfactor) ); 464 } 465 } 466 467 return SCIP_OKAY; 468 } 469 470 /** create the extra constraint of local branching and add it to subscip */ 471 static 472 SCIP_RETCODE addLocalBranchingConstraint( 473 SCIP* scip, /**< SCIP data structure of the original problem */ 474 SCIP* probingscip, /**< SCIP data structure of the subproblem */ 475 SCIP_HASHMAP* varmapfw, /**< mapping of SCIP variables to sub-SCIP variables */ 476 SCIP_SOL* bestsol, /**< SCIP solution */ 477 SCIP_Real neighborhoodsize /**< rhs for LB constraint */ 478 ) 479 { 480 SCIP_CONS* cons; /* local branching constraint to create */ 481 SCIP_VAR** consvars; 482 SCIP_VAR** vars; 483 484 int nbinvars; 485 int nconsvars; 486 int i; 487 SCIP_Real lhs; 488 SCIP_Real rhs; 489 SCIP_Real* consvals; 490 char consname[SCIP_MAXSTRLEN]; 491 492 (void) SCIPsnprintf(consname, SCIP_MAXSTRLEN, "%s_localbranchcons", SCIPgetProbName(scip)); 493 494 /* get vars data */ 495 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, NULL, NULL, NULL) ); 496 /* memory allocation */ 497 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nbinvars) ); 498 SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nbinvars) ); 499 nconsvars = 0; 500 501 /* set initial left and right hand sides of local branching constraint */ 502 lhs = 0.0; 503 rhs = neighborhoodsize; 504 505 /* create the distance (to incumbent) function of the binary variables */ 506 for( i = 0; i < nbinvars; i++ ) 507 { 508 SCIP_Real solval; 509 510 solval = SCIPgetSolVal(scip, bestsol, vars[i]); 511 assert( SCIPisFeasIntegral(scip, solval) ); 512 513 /* is variable i part of the binary support of closest sol? */ 514 if( SCIPisFeasEQ(scip,solval,1.0) ) 515 { 516 consvals[nconsvars] = -1.0; 517 rhs -= 1.0; 518 lhs -= 1.0; 519 } 520 else 521 consvals[nconsvars] = 1.0; 522 consvars[nconsvars] = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, vars[i]); 523 if( consvars[nconsvars] == NULL ) 524 continue; 525 SCIP_CALL( SCIPchgVarObj(probingscip, consvars[nconsvars], consvals[nconsvars]) ); 526 assert( SCIPvarGetType(consvars[nconsvars]) == SCIP_VARTYPE_BINARY ); 527 ++nconsvars; 528 } 529 530 /* creates localbranching constraint and adds it to subscip */ 531 SCIP_CALL( SCIPcreateConsLinear(probingscip, &cons, consname, nconsvars, consvars, consvals, 532 lhs, rhs, FALSE, FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) ); 533 SCIP_CALL( SCIPaddCons(probingscip, cons) ); 534 SCIP_CALL( SCIPreleaseCons(probingscip, &cons) ); 535 536 /* free local memory */ 537 SCIPfreeBufferArray(scip, &consvals); 538 SCIPfreeBufferArray(scip, &consvars); 539 540 return SCIP_OKAY; 541 } 542 543 /** creates new solutions for the original problem by copying the solutions of the subproblem */ 544 static 545 SCIP_RETCODE createNewSols( 546 SCIP* scip, /**< original SCIP data structure */ 547 SCIP* subscip, /**< SCIP structure of the subproblem */ 548 SCIP_HASHMAP* varmapfw, /**< mapping of SCIP variables to sub-SCIP variables */ 549 SCIP_HEUR* heur, /**< heuristic structure */ 550 SCIP_Bool* success /**< used to store whether new solution was found or not */ 551 ) 552 { 553 SCIP_VAR** vars; /* the original problem's variables */ 554 int nvars; 555 SCIP_VAR** subvars; 556 int i; 557 558 assert(scip != NULL); 559 assert(subscip != NULL); 560 561 /* get variables' data */ 562 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) ); 563 564 /* for copying a solution we need an explicit mapping */ 565 SCIP_CALL( SCIPallocBufferArray(scip, &subvars, nvars) ); 566 for( i = 0; i < nvars; i++ ) 567 subvars[i] = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, vars[i]); 568 569 SCIP_CALL( SCIPtranslateSubSols(scip, subscip, heur, subvars, success, NULL) ); 570 571 SCIPfreeBufferArray(scip, &subvars); 572 573 return SCIP_OKAY; 574 } 575 576 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 577 static 578 SCIP_DECL_HEURCOPY(heurCopyFeaspump) 579 { 580 assert(scip != NULL); 581 assert(heur != NULL); 582 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 583 584 /* call inclusion method of primal heuristic */ 585 SCIP_CALL( SCIPincludeHeurFeaspump(scip) ); 586 587 return SCIP_OKAY; 588 } 589 590 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 591 static 592 SCIP_DECL_HEURFREE(heurFreeFeaspump) 593 { 594 SCIP_HEURDATA* heurdata; 595 596 assert(heur != NULL); 597 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 598 assert(scip != NULL); 599 600 /* free heuristic data */ 601 heurdata = SCIPheurGetData(heur); 602 assert(heurdata != NULL); 603 SCIPfreeBlockMemory(scip, &heurdata); 604 SCIPheurSetData(heur, NULL); 605 606 return SCIP_OKAY; 607 } 608 609 610 /** initialization method of primal heuristic (called after problem was transformed) */ 611 static 612 SCIP_DECL_HEURINIT(heurInitFeaspump) 613 { 614 SCIP_HEURDATA* heurdata; 615 616 assert(heur != NULL); 617 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 618 619 /* get heuristic data */ 620 heurdata = SCIPheurGetData(heur); 621 assert(heurdata != NULL); 622 623 /* create working solution */ 624 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) ); 625 SCIP_CALL( SCIPcreateSol(scip, &heurdata->roundedsol, heur) ); 626 627 /* initialize data */ 628 heurdata->nlpiterations = 0; 629 heurdata->nsuccess = 0; 630 631 /* create random number generator */ 632 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, 633 DEFAULT_RANDSEED, TRUE) ); 634 635 return SCIP_OKAY; 636 } 637 638 /** deinitialization method of primal heuristic (called before transformed problem is freed) */ 639 static 640 SCIP_DECL_HEUREXIT(heurExitFeaspump) 641 { 642 SCIP_HEURDATA* heurdata; 643 644 assert(heur != NULL); 645 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 646 647 /* get heuristic data */ 648 heurdata = SCIPheurGetData(heur); 649 assert(heurdata != NULL); 650 651 /* free working solution */ 652 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) ); 653 SCIP_CALL( SCIPfreeSol(scip, &heurdata->roundedsol) ); 654 655 /* free random number generator */ 656 SCIPfreeRandom(scip, &heurdata->randnumgen); 657 658 return SCIP_OKAY; 659 } 660 661 662 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */ 663 static 664 SCIP_DECL_HEURINITSOL(heurInitsolFeaspump) 665 { 666 SCIP_HEURDATA* heurdata; 667 668 heurdata = SCIPheurGetData(heur); 669 assert(heurdata != NULL); 670 671 /* if the heuristic is called at the root node, we may want to be called directly after the initial root LP solve */ 672 if( heurdata->beforecuts && SCIPheurGetFreqofs(heur) == 0 ) 673 SCIPheurSetTimingmask(heur, SCIP_HEURTIMING_DURINGLPLOOP); 674 675 return SCIP_OKAY; 676 } 677 678 679 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */ 680 static 681 SCIP_DECL_HEUREXITSOL(heurExitsolFeaspump) 682 { 683 /* reset the timing mask to its default value */ 684 SCIPheurSetTimingmask(heur, HEUR_TIMING); 685 686 return SCIP_OKAY; 687 } 688 689 /** calculates an adjusted maximal number of LP iterations */ 690 static 691 SCIP_Longint adjustedMaxNLPIterations( 692 SCIP_Longint maxnlpiterations, /**< regular maximal number of LP iterations */ 693 SCIP_Longint nsolsfound, /**< total number of solutions found so far by SCIP */ 694 int nstallloops /**< current number of stalling rounds */ 695 ) 696 { 697 if( nstallloops <= 1 ) 698 { 699 if( nsolsfound == 0 ) 700 return 4*maxnlpiterations; 701 else 702 return 2*maxnlpiterations; 703 } 704 else 705 return maxnlpiterations; 706 } 707 708 /** execution method of primal heuristic */ 709 static 710 SCIP_DECL_HEUREXEC(heurExecFeaspump) 711 { 712 SCIP_HEURDATA* heurdata; 713 SCIP_SOL* tmpsol; /* only used for swapping */ 714 SCIP_SOL** lastroundedsols;/* solutions of the last pumping rounds (depending on heurdata->cyclelength) */ 715 SCIP_SOL* closestsol; /* rounded solution closest to the LP relaxation: used for stage3 */ 716 SCIP_Real* lastalphas; /* alpha values associated to solutions in lastroundedsols */ 717 718 SCIP* probingscip; /* copied SCIP structure, used for round-and-propagate loop of feasibility pump 2.0 */ 719 SCIP_HASHMAP* varmapfw; /* mapping of SCIP variables to sub-SCIP variables */ 720 721 SCIP_VAR** vars; 722 SCIP_VAR** pseudocands; 723 SCIP_VAR** tmppseudocands; 724 SCIP_VAR** mostfracvars; /* the 30 most fractional variables, needed to avoid 1-cycles */ 725 SCIP_VAR* var; 726 727 SCIP_Real* mostfracvals; /* the values of the variables above */ 728 SCIP_Real oldsolval; /* one value of the last solution */ 729 SCIP_Real solval; /* one value of the actual solution */ 730 SCIP_Real frac; /* the fractional part of the value above */ 731 SCIP_Real objfactor; /* factor by which the regard of the objective is decreased in each round, in [0,0.99] */ 732 SCIP_Real alpha; /* factor how the original objective is regarded, used for convex combination of two functions */ 733 SCIP_Real objnorm; /* Euclidean norm of the objective function, used for scaling */ 734 SCIP_Real scalingfactor; /* factor to scale the original objective function with */ 735 SCIP_Real mindistance; /* distance of the closest rounded solution from the LP relaxation: used for stage3 */ 736 737 SCIP_Longint nlpiterations; /* number of LP iterations done during one pumping round */ 738 SCIP_Longint maxnlpiterations; /* maximum number of LP iterations for this heuristic */ 739 SCIP_Longint nsolsfound; /* number of solutions found by this heuristic */ 740 SCIP_Longint ncalls; /* number of calls of this heuristic */ 741 SCIP_Longint nbestsolsfound; /* current total number of best solution updates in SCIP */ 742 743 SCIP_LPSOLSTAT lpsolstat; /* status of the LP solution */ 744 745 int nvars; /* number of variables */ 746 int nbinvars; /* number of 0-1-variables */ 747 int nintvars; /* number of integer variables */ 748 int nfracs; /* number of fractional variables updated after each pumping round*/ 749 int nflipcands; /* how many flipcands (most frac. var.) have been found */ 750 int npseudocands; 751 int maxnflipcands; /* maximal number of candidates to flip in the current pumping round */ 752 int nloops; /* how many pumping rounds have been made */ 753 int maxflips; /* maximum number of flips, if a 1-cycle is found (depending on heurdata->minflips) */ 754 int maxloops; /* maximum number of pumping rounds */ 755 int nstallloops; /* number of loops without reducing the current best number of factional variables */ 756 int maxstallloops; /* maximal number of allowed stalling loops */ 757 int bestnfracs; /* best number of fractional variables */ 758 int i; 759 int j; 760 761 SCIP_Bool success; 762 SCIP_Bool lperror; 763 SCIP_Bool* cycles; /* are there short cycles */ 764 765 SCIP_RETCODE retcode; 766 767 assert(heur != NULL); 768 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 769 assert(scip != NULL); 770 assert(result != NULL); 771 assert(SCIPhasCurrentNodeLP(scip)); 772 773 *result = SCIP_DELAYED; 774 775 /* do not call heuristic of node was already detected to be infeasible */ 776 if( nodeinfeasible ) 777 return SCIP_OKAY; 778 779 /* only call heuristic, if an optimal LP solution is at hand */ 780 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 781 return SCIP_OKAY; 782 783 /* only call heuristic, if the LP objective value is smaller than the cutoff bound */ 784 if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) ) 785 return SCIP_OKAY; 786 787 /* only call heuristic, if the LP solution is basic (which allows fast resolve in diving) */ 788 if( !SCIPisLPSolBasic(scip) ) 789 return SCIP_OKAY; 790 791 *result = SCIP_DIDNOTRUN; 792 793 /* don't dive two times at the same node */ 794 if( SCIPgetLastDivenode(scip) == SCIPgetNNodes(scip) && SCIPgetDepth(scip) > 0 ) 795 return SCIP_OKAY; 796 797 /* only call feaspump once at the root */ 798 if( SCIPgetDepth(scip) == 0 && SCIPheurGetNCalls(heur) > 0 ) 799 return SCIP_OKAY; 800 801 /* reset the timing mask to its default value (at the root node it could be different) */ 802 SCIPheurSetTimingmask(heur, HEUR_TIMING); 803 804 /* only call the heuristic before cutting planes if we do not have an incumbent and no pricer exists */ 805 if( heurtiming == SCIP_HEURTIMING_DURINGLPLOOP && SCIPgetNSolsFound(scip) > 0 && SCIPgetNPricers(scip) == 0 ) 806 return SCIP_OKAY; 807 808 /* get heuristic's data */ 809 heurdata = SCIPheurGetData(heur); 810 assert(heurdata != NULL); 811 812 /* only apply heuristic, if only a few solutions have been found and no pricer exists */ 813 if( heurdata->maxsols >= 0 && SCIPgetNSolsFound(scip) > heurdata->maxsols && SCIPgetNPricers(scip) == 0 ) 814 return SCIP_OKAY; 815 816 /* get all variables of LP and number of fractional variables in LP solution that should be integral */ 817 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) ); 818 nfracs = SCIPgetNLPBranchCands(scip); 819 assert(0 <= nfracs && nfracs <= nbinvars + nintvars); 820 if( nfracs == 0 ) 821 return SCIP_OKAY; 822 823 /* calculate the maximal number of LP iterations until heuristic is aborted */ 824 nlpiterations = SCIPgetNLPIterations(scip); 825 ncalls = SCIPheurGetNCalls(heur); 826 nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + heurdata->nsuccess; 827 maxnlpiterations = (SCIP_Longint)((1.0 + 10.0*(nsolsfound+1.0)/(ncalls+1.0)) * heurdata->maxlpiterquot * nlpiterations); 828 maxnlpiterations += heurdata->maxlpiterofs; 829 830 /* don't try to dive, if we took too many LP iterations during diving */ 831 if( heurdata->nlpiterations >= maxnlpiterations ) 832 return SCIP_OKAY; 833 834 /* at the first root call, allow more iterations if there is no feasible solution yet */ 835 if( SCIPheurGetNCalls(heur) == 0 && SCIPgetNSolsFound(scip) == 0 && SCIPgetDepth(scip) == 0 ) 836 maxnlpiterations += nlpiterations; 837 838 /* allow at least a certain number of LP iterations in this dive */ 839 maxnlpiterations = MAX(maxnlpiterations, heurdata->nlpiterations + MINLPITER); 840 841 /* calculate maximal number of flips and loops */ 842 maxflips = 3*heurdata->minflips; 843 maxloops = (heurdata->maxloops == -1 ? INT_MAX : heurdata->maxloops); 844 maxstallloops = (heurdata->maxstallloops == -1 ? INT_MAX : heurdata->maxstallloops); 845 846 SCIPdebugMsg(scip, "executing feasibility pump heuristic, nlpiters=%" SCIP_LONGINT_FORMAT ", maxnlpit:%" SCIP_LONGINT_FORMAT ", maxflips:%d \n", 847 nlpiterations, maxnlpiterations, maxflips); 848 849 *result = SCIP_DIDNOTFIND; 850 851 probingscip = NULL; 852 varmapfw = NULL; 853 854 if( heurdata->usefp20 ) 855 { 856 SCIP_Bool valid; 857 858 /* ignore value of valid */ 859 SCIP_CALL( setupProbingSCIP(scip, &probingscip, &varmapfw, heurdata->copycuts, &valid) ); 860 861 if( probingscip != NULL ) 862 { 863 SCIP_CALL( setupSCIPparamsFP2(scip, probingscip) ); 864 865 retcode = SCIPsolve(probingscip); 866 867 /* errors in solving the subproblem should not kill the overall solving process; 868 * hence, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. */ 869 if( retcode != SCIP_OKAY ) 870 { 871 #ifndef NDEBUG 872 SCIP_CALL( retcode ); 873 #endif 874 SCIPwarningMessage(scip, "Error while solving subproblem in feaspump heuristic; sub-SCIP terminated with code <%d>\n", retcode); 875 876 /* free hash map and copied SCIP */ 877 SCIPhashmapFree(&varmapfw); 878 SCIP_CALL( SCIPfree(&probingscip) ); 879 return SCIP_OKAY; 880 } 881 882 if( SCIPgetStage(probingscip) != SCIP_STAGE_SOLVING) 883 { 884 SCIP_STATUS probingstatus = SCIPgetStatus(probingscip); 885 886 if( probingstatus == SCIP_STATUS_OPTIMAL ) 887 { 888 assert( SCIPgetNSols(probingscip) > 0 ); 889 SCIP_CALL( createNewSols(scip, probingscip, varmapfw, heur, &success) ); 890 if( success ) 891 *result = SCIP_FOUNDSOL; 892 } 893 894 /* free hash map and copied SCIP */ 895 SCIPhashmapFree(&varmapfw); 896 SCIP_CALL( SCIPfree(&probingscip) ); 897 return SCIP_OKAY; 898 } 899 SCIP_CALL( SCIPsetLongintParam(probingscip, "limits/nodes", 2LL) ); 900 901 /* set SCIP into probing mode and create root node of the probing tree */ 902 SCIP_CALL( SCIPstartProbing(probingscip) ); 903 904 /* this should always be fulfilled */ 905 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(probingscip)); 906 907 SCIP_CALL( SCIPnewProbingNode(probingscip) ); 908 909 SCIPdebugMsg(scip, "successfully copied SCIP instance -> feasibility pump 2.0 can be used.\n"); 910 } 911 else 912 { 913 assert(varmapfw == NULL); 914 915 SCIPdebugMsg(scip, "SCIP reached the depth limit -> skip heuristic\n"); 916 return SCIP_OKAY; 917 } 918 } /*lint !e438*/ 919 920 /* memory allocation */ 921 SCIP_CALL( SCIPallocBufferArray(scip, &mostfracvars, maxflips) ); 922 SCIP_CALL( SCIPallocBufferArray(scip, &mostfracvals, maxflips) ); 923 SCIP_CALL( SCIPallocBufferArray(scip, &lastroundedsols, heurdata->cyclelength) ); 924 SCIP_CALL( SCIPallocBufferArray(scip, &lastalphas, heurdata->cyclelength) ); 925 SCIP_CALL( SCIPallocBufferArray(scip, &cycles, heurdata->cyclelength) ); 926 927 for( j = 0; j < heurdata->cyclelength; j++ ) 928 { 929 SCIP_CALL( SCIPcreateSol(scip, &lastroundedsols[j], heur) ); 930 } 931 932 closestsol = NULL; 933 if( heurdata->stage3 ) 934 { 935 SCIP_CALL( SCIPcreateSol(scip, &closestsol, heur) ); 936 } 937 938 /* start diving */ 939 SCIP_CALL( SCIPstartDive(scip) ); 940 941 /* lp was solved optimal */ 942 lperror = FALSE; 943 lpsolstat = SCIP_LPSOLSTAT_OPTIMAL; 944 945 /* pumping rounds */ 946 nsolsfound = SCIPgetNBestSolsFound(scip); 947 if( heurdata->objfactor == 1.0 ) 948 objfactor = MIN(1.0 - 0.1 / (SCIP_Real)(1 + nsolsfound), 0.999); 949 else 950 objfactor = heurdata->objfactor; 951 952 /* scale distance function and original objective to the same norm */ 953 objnorm = SCIPgetObjNorm(scip); 954 objnorm = MAX(objnorm, 1.0); 955 scalingfactor = SQRT((SCIP_Real)(nbinvars + nintvars)) / objnorm; 956 957 /* data initialization */ 958 alpha = heurdata->alpha; 959 nloops = 0; 960 nstallloops = 0; 961 nbestsolsfound = SCIPgetNBestSolsFound(scip); 962 bestnfracs = INT_MAX; 963 mindistance = SCIPinfinity(scip); 964 965 /* pumping loop */ 966 while( nfracs > 0 967 && heurdata->nlpiterations < adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops) 968 && nloops < maxloops && nstallloops < maxstallloops 969 && !SCIPisStopped(scip) ) 970 { 971 int minimum; 972 SCIP_Real* pseudocandsfrac; 973 SCIP_Longint nlpiterationsleft; 974 int iterlimit; 975 976 /* decrease convex combination scalar */ 977 nloops++; 978 alpha *= objfactor; 979 980 SCIPdebugMsg(scip, "feasibility pump loop %d: %d fractional variables (alpha: %.4f, stall: %d/%d)\n", 981 nloops, nfracs, alpha, nstallloops, maxstallloops); 982 983 success = FALSE; 984 985 SCIP_CALL( SCIPlinkLPSol(scip, heurdata->roundedsol) ); 986 987 /* randomly choose maximum number of variables to flip in current pumping round in case of a 1-cycle */ 988 maxnflipcands = SCIPrandomGetInt(heurdata->randnumgen, MIN(nfracs/2+1, heurdata->minflips), MIN(nfracs, maxflips)); 989 nflipcands = 0; 990 991 /* get all unfixed integer variables */ 992 SCIP_CALL( SCIPgetPseudoBranchCands(scip, &tmppseudocands, &npseudocands, NULL) ); 993 SCIP_CALL( SCIPduplicateBufferArray(scip, &pseudocands, tmppseudocands, npseudocands) ); 994 995 /* get array of all fractional variables and sort it w.r.t. their fractionalities */ 996 if( heurdata->usefp20 ) 997 { 998 SCIP_CALL( SCIPallocBufferArray(scip, &pseudocandsfrac, npseudocands) ); 999 1000 for( i = 0; i < npseudocands; i++ ) 1001 { 1002 frac = SCIPfeasFrac(scip, SCIPvarGetLPSol(pseudocands[i])); 1003 pseudocandsfrac[i] = MIN(frac, 1.0-frac); /* always a number between 0 and 0.5 */ 1004 if( SCIPvarGetType(pseudocands[i]) == SCIP_VARTYPE_BINARY ) 1005 pseudocandsfrac[i] -= 10.0; /* binaries always come first */ 1006 } 1007 SCIPsortRealPtr(pseudocandsfrac, (void**)pseudocands, npseudocands); 1008 SCIPfreeBufferArray(scip, &pseudocandsfrac); 1009 1010 SCIPdebugMsg(scip, "iteratively fix and propagate variables\n"); 1011 } 1012 1013 for( i = 0; i < npseudocands; i++ ) 1014 { 1015 SCIP_VAR* probingvar; 1016 SCIP_Bool infeasible; 1017 SCIP_Longint ndomreds; 1018 1019 var = pseudocands[i]; 1020 1021 /* round the LP solution */ 1022 solval = SCIPvarGetLPSol(var); 1023 frac = SCIPfeasFrac(scip, solval); 1024 1025 /* round randomly if the value is close to 0.5 */ 1026 if( SCIPisEQ(scip, frac, 0.5) ) 1027 { 1028 if( SCIPrandomGetReal(heurdata->randnumgen, 0.0, 1.0) <= 0.5 ) 1029 solval = SCIPfloor(scip, solval); 1030 else 1031 solval = SCIPceil(scip, solval); 1032 } 1033 else 1034 solval = SCIPfloor(scip, solval + 0.5); 1035 1036 /* ensure, that the fixing value is inside the local domains */ 1037 if( heurdata->usefp20 ) 1038 { 1039 SCIP_Real lbprobing; 1040 SCIP_Real ubprobing; 1041 1042 probingvar = (SCIP_VAR*) SCIPhashmapGetImage(varmapfw, var); 1043 /* skip if variable does not exist in probingscip */ 1044 if( probingvar != NULL ) 1045 { 1046 lbprobing = SCIPvarGetLbLocal(probingvar); 1047 ubprobing = SCIPvarGetUbLocal(probingvar); 1048 1049 solval = MAX(solval, lbprobing); 1050 solval = MIN(solval, ubprobing); 1051 1052 /* fix the variable and propagate the domain change */ 1053 if( !SCIPisFeasEQ(probingscip, lbprobing, ubprobing) && SCIPvarIsActive(SCIPvarGetTransVar(probingvar)) ) 1054 { 1055 assert(SCIPisFeasLE(probingscip, lbprobing, ubprobing)); 1056 SCIPdebugMsg(scip, "try to fix variable <%s> (domain [%f,%f] to %f\n", SCIPvarGetName(probingvar), lbprobing, ubprobing, 1057 solval); 1058 SCIP_CALL( SCIPfixVarProbing(probingscip, probingvar, solval) ); 1059 SCIP_CALL( SCIPpropagateProbing(probingscip, -1, &infeasible, &ndomreds) ); 1060 SCIPdebugMsg(scip, " -> reduced %" SCIP_LONGINT_FORMAT " domains\n", ndomreds); 1061 1062 if( infeasible ) 1063 { 1064 SCIPdebugMsg(scip, " -> infeasible!\n"); 1065 SCIP_CALL( SCIPbacktrackProbing(probingscip, 0) ); 1066 } 1067 } 1068 else 1069 { 1070 SCIPdebugMsg(scip, "variable <%s> is already fixed to %f\n", SCIPvarGetName(probingvar), solval); 1071 } 1072 } 1073 } 1074 1075 SCIP_CALL( updateVariableRounding(scip, heurdata, var, solval, alpha, scalingfactor) ); 1076 1077 /* check whether the variable is one of the most fractionals and label if so */ 1078 if( SCIPisFeasPositive(scip, frac) ) 1079 insertFlipCand(mostfracvars, mostfracvals, &nflipcands, maxnflipcands, var, frac); 1080 } 1081 1082 if( heurdata->usefp20 ) 1083 { 1084 SCIP_CALL( SCIPbacktrackProbing(probingscip, 0) ); 1085 } 1086 1087 /* change objective coefficients for continuous variables */ 1088 for( i = nbinvars+nintvars; i < nvars; i++ ) 1089 { 1090 SCIP_CALL( SCIPchgVarObjDive(scip, vars[i], alpha * SCIPvarGetObj(vars[i])) ); 1091 } 1092 1093 SCIPfreeBufferArray(scip, &pseudocands); 1094 1095 /* initialize cycle check */ 1096 minimum = MIN(heurdata->cyclelength, nloops-1); 1097 for( j = 0; j < heurdata->cyclelength; j++ ) 1098 cycles[j] = (nloops > j+1) && (REALABS(lastalphas[j] - alpha) < heurdata->alphadiff); 1099 1100 /* check for j-cycles */ 1101 for( i = 0; i < nbinvars+nintvars; i++ ) 1102 { 1103 solval = SCIPgetSolVal(scip, heurdata->roundedsol, vars[i]); 1104 1105 /* cycles exist, iff all solution values are equal */ 1106 for( j = 0; j < minimum; j++ ) 1107 { 1108 oldsolval = SCIPgetSolVal(scip, lastroundedsols[j], vars[i]); 1109 cycles[j] = cycles[j] && SCIPisFeasEQ(scip, solval, oldsolval); 1110 } 1111 } 1112 1113 /* force to flip variables at random after a couple of pumping rounds, 1114 * or if a new best solution in the current region has been found 1115 */ 1116 assert(heurdata->perturbfreq > 0); 1117 if( nloops % heurdata->perturbfreq == 0 || (heurdata->pertsolfound && SCIPgetNBestSolsFound(scip) > nbestsolsfound) ) 1118 { 1119 SCIPdebugMsg(scip, " -> random perturbation\n"); 1120 SCIP_CALL( handleCycle(scip, heurdata, vars, nintvars+nbinvars, alpha, scalingfactor) ); 1121 nbestsolsfound = SCIPgetNBestSolsFound(scip); 1122 } 1123 else 1124 { 1125 minimum = MIN(heurdata->cyclelength, nloops-1); 1126 1127 for( j = 0; j < minimum; j++ ) 1128 { 1129 /* if we got the same rounded solution as in some step before, we have to flip some variables */ 1130 if( cycles[j] ) 1131 { 1132 /* 1-cycles have a special flipping rule (flip most fractional variables) */ 1133 if( j == 0 ) 1134 { 1135 SCIPdebugMsg(scip, " -> avoiding 1-cycle: flipping %d candidates\n", nflipcands); 1136 SCIP_CALL( handle1Cycle(scip, heurdata, mostfracvars, nflipcands, alpha, scalingfactor) ); 1137 } 1138 else 1139 { 1140 SCIPdebugMsg(scip, " -> avoiding %d-cycle by random flip\n", j+1); 1141 SCIP_CALL( handleCycle(scip, heurdata, vars, nintvars+nbinvars, alpha, scalingfactor) ); 1142 } 1143 break; 1144 } 1145 } 1146 } 1147 1148 /* the LP with the new (distance) objective is solved */ 1149 nlpiterations = SCIPgetNLPIterations(scip); 1150 nlpiterationsleft = adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops) - heurdata->nlpiterations; 1151 iterlimit = MAX((int)nlpiterationsleft, MINLPITER); 1152 SCIPdebugMsg(scip, " -> solve LP with iteration limit %d\n", iterlimit); 1153 1154 if( heurdata->stage3 ) 1155 { 1156 SCIP_CALL( SCIPunlinkSol(scip, heurdata->roundedsol) ); 1157 } 1158 1159 retcode = SCIPsolveDiveLP(scip, iterlimit, &lperror, NULL); 1160 lpsolstat = SCIPgetLPSolstat(scip); 1161 1162 /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic. 1163 * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. 1164 */ 1165 if( retcode != SCIP_OKAY ) 1166 { 1167 #ifndef NDEBUG 1168 if( lpsolstat != SCIP_LPSOLSTAT_UNBOUNDEDRAY ) 1169 { 1170 SCIP_CALL( retcode ); 1171 } 1172 #endif 1173 SCIPwarningMessage(scip, "Error while solving LP in Feaspump heuristic; LP solve terminated with code <%d>\n", retcode); 1174 SCIPwarningMessage(scip, "This does not affect the remaining solution procedure --> continue\n"); 1175 } 1176 1177 /* update iteration count */ 1178 heurdata->nlpiterations += SCIPgetNLPIterations(scip) - nlpiterations; 1179 SCIPdebugMsg(scip, " -> number of iterations: %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ", lperror=%u, lpsolstat=%d\n", 1180 heurdata->nlpiterations, adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops), lperror, lpsolstat); 1181 1182 /* check whether LP was solved optimal */ 1183 if( lperror || lpsolstat != SCIP_LPSOLSTAT_OPTIMAL ) 1184 break; 1185 1186 if( heurdata->stage3 ) 1187 { 1188 SCIP_Real distance; /* distance of the current rounded solution from the LP solution */ 1189 1190 assert(closestsol != NULL); 1191 1192 /* calculate distance */ 1193 distance = 0.0; 1194 for( i = 0; i < nbinvars+nintvars; i++ ) 1195 { 1196 SCIP_Real roundedval; 1197 SCIP_Real lpval; 1198 1199 roundedval = SCIPgetSolVal(scip, heurdata->roundedsol, vars[i]); 1200 lpval = SCIPvarGetLPSol(vars[i]); 1201 distance += REALABS(roundedval - lpval); 1202 } 1203 1204 /* copy solution and update minimum distance */ 1205 if( SCIPisLT(scip, distance, mindistance) ) 1206 { 1207 for( i = 0; i < nbinvars+nintvars; i++ ) 1208 { 1209 assert(SCIPisIntegral(scip,SCIPgetSolVal(scip, heurdata->roundedsol, vars[i]))); 1210 SCIP_CALL( SCIPsetSolVal(scip, closestsol, vars[i], SCIPgetSolVal(scip, heurdata->roundedsol, vars[i])) ); 1211 } 1212 mindistance = distance; 1213 } 1214 } 1215 1216 /* swap the last solutions */ 1217 SCIP_CALL( SCIPunlinkSol(scip, heurdata->roundedsol) ); 1218 tmpsol = lastroundedsols[heurdata->cyclelength-1]; 1219 for( j = heurdata->cyclelength-1; j > 0; j-- ) 1220 { 1221 lastroundedsols[j] = lastroundedsols[j-1]; 1222 lastalphas[j] = lastalphas[j-1]; 1223 } 1224 lastroundedsols[0] = heurdata->roundedsol; 1225 lastalphas[0] = alpha; 1226 heurdata->roundedsol = tmpsol; 1227 1228 /* check for improvement in number of fractionals */ 1229 nfracs = SCIPgetNLPBranchCands(scip); 1230 if( nfracs < bestnfracs ) 1231 { 1232 bestnfracs = nfracs; 1233 nstallloops = 0; 1234 } 1235 else 1236 nstallloops++; 1237 1238 SCIPdebugMsg(scip, " -> loop finished: %d fractional variables (stall: %d/%d, iterations: %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ")\n", 1239 nfracs, nstallloops, maxstallloops, heurdata->nlpiterations, adjustedMaxNLPIterations(maxnlpiterations, nsolsfound, nstallloops)); 1240 } 1241 1242 /* try final solution, if no more fractional variables are left */ 1243 if( nfracs == 0 && !lperror && lpsolstat == SCIP_LPSOLSTAT_OPTIMAL ) 1244 { 1245 success = FALSE; 1246 1247 SCIP_CALL( SCIPlinkLPSol(scip, heurdata->sol) ); 1248 SCIPdebugMsg(scip, "feasibility pump found solution (%d fractional variables)\n", nfracs); 1249 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) ); 1250 if( success ) 1251 *result = SCIP_FOUNDSOL; 1252 } 1253 1254 /* end diving */ 1255 SCIP_CALL( SCIPendDive(scip) ); 1256 1257 /* end probing in order to be able to apply stage 3 */ 1258 if( heurdata->usefp20 ) 1259 { 1260 SCIP_CALL( SCIPendProbing(probingscip) ); 1261 } 1262 1263 /* only do stage 3 if we have not found a solution yet */ 1264 /* only do stage 3 if the distance of the closest infeasible solution to the polyhedron is below a certain threshold */ 1265 if( heurdata->stage3 && (*result != SCIP_FOUNDSOL) && SCIPisLE(scip, mindistance, (SCIP_Real) heurdata->neighborhoodsize) ) 1266 { 1267 SCIP_Bool cancopy; 1268 assert(closestsol != NULL); 1269 assert(!SCIPisInfinity(scip, mindistance) || nloops == 0); 1270 1271 /* if we do not use feasibility pump 2.0, we have not created a copied SCIP instance yet */ 1272 if( heurdata->usefp20 ) 1273 { 1274 assert(probingscip != NULL); 1275 SCIP_CALL( SCIPfreeTransform(probingscip) ); 1276 } 1277 else 1278 { 1279 assert(probingscip == NULL); 1280 SCIP_CALL( setupProbingSCIP(scip, &probingscip, &varmapfw, heurdata->copycuts, &success) ); 1281 } 1282 1283 /* check whether there is enough time and memory left */ 1284 SCIP_CALL( SCIPcheckCopyLimits(scip, &cancopy) ); 1285 1286 if( cancopy ) 1287 { 1288 SCIP_CALL( setupSCIPparamsStage3(scip, probingscip) ); 1289 1290 /* the neighborhood size is double the distance plus another ten percent */ 1291 mindistance = SCIPceil(scip, 2.2*mindistance); 1292 1293 SCIP_CALL( addLocalBranchingConstraint(scip, probingscip, varmapfw, closestsol, mindistance) ); 1294 1295 /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic. 1296 * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. 1297 */ 1298 SCIP_CALL_ABORT( SCIPsolve(probingscip) ); 1299 1300 /* check, whether a solution was found */ 1301 if( SCIPgetNSols(probingscip) > 0 ) 1302 { 1303 success = FALSE; 1304 SCIP_CALL( createNewSols(scip, probingscip, varmapfw, heur, &success) ); 1305 if( success ) 1306 *result = SCIP_FOUNDSOL; 1307 } 1308 } 1309 } 1310 1311 if( *result == SCIP_FOUNDSOL ) 1312 heurdata->nsuccess++; 1313 1314 /* free hash map and copied SCIP */ 1315 if( varmapfw != NULL ) 1316 SCIPhashmapFree(&varmapfw); 1317 1318 if( probingscip != NULL ) 1319 { 1320 SCIP_CALL( SCIPfree(&probingscip) ); 1321 } 1322 1323 if( heurdata->stage3 ) 1324 { 1325 SCIP_CALL( SCIPfreeSol(scip, &closestsol) ); 1326 } 1327 1328 /* free memory */ 1329 for( j = 0; j < heurdata->cyclelength; j++ ) 1330 { 1331 SCIP_CALL( SCIPfreeSol(scip, &lastroundedsols[j]) ); 1332 } 1333 1334 SCIPfreeBufferArray(scip, &cycles); 1335 SCIPfreeBufferArray(scip, &lastalphas); 1336 SCIPfreeBufferArray(scip, &lastroundedsols); 1337 SCIPfreeBufferArray(scip, &mostfracvals); 1338 SCIPfreeBufferArray(scip, &mostfracvars); 1339 1340 SCIPdebugMsg(scip, "feasibility pump finished [%d iterations done].\n", nloops); 1341 1342 #ifdef SCIP_STATISTIC 1343 if( nfracs == 0 ) 1344 { 1345 double objval; 1346 double primalBound; 1347 1348 objval = SCIPgetSolOrigObj(scip, heurdata->sol); 1349 primalBound = SCIPgetPrimalbound(scip); 1350 SCIPstatisticMessage("feasibility pump found: 1, objval: %f, iterations: %d, primal bound: %f\n", objval, nloops, primalBound); 1351 } 1352 else 1353 { 1354 double primalBound; 1355 1356 primalBound = SCIPgetPrimalbound(scip); 1357 SCIPstatisticMessage("feasibility pump found: 0, objval: +inf, iterations: %d, primal bound: %f\n", nloops, primalBound); 1358 } 1359 1360 #endif /* SCIP_STATISTIC */ 1361 1362 return SCIP_OKAY; 1363 } 1364 1365 1366 /* 1367 * primal heuristic specific interface methods 1368 */ 1369 1370 /** creates the feaspump primal heuristic and includes it in SCIP */ 1371 SCIP_RETCODE SCIPincludeHeurFeaspump( 1372 SCIP* scip /**< SCIP data structure */ 1373 ) 1374 { 1375 SCIP_HEURDATA* heurdata; 1376 SCIP_HEUR* heur; 1377 1378 /* create Feaspump primal heuristic data */ 1379 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 1380 1381 /* include primal heuristic */ 1382 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 1383 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 1384 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecFeaspump, heurdata) ); 1385 1386 assert(heur != NULL); 1387 1388 /* set non-NULL pointers to callback methods */ 1389 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyFeaspump) ); 1390 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeFeaspump) ); 1391 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitFeaspump) ); 1392 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitFeaspump) ); 1393 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolFeaspump) ); 1394 SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolFeaspump) ); 1395 1396 /* add feaspump primal heuristic parameters */ 1397 SCIP_CALL( SCIPaddRealParam(scip, 1398 "heuristics/" HEUR_NAME "/maxlpiterquot", 1399 "maximal fraction of diving LP iterations compared to node LP iterations", 1400 &heurdata->maxlpiterquot, FALSE, DEFAULT_MAXLPITERQUOT, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 1401 SCIP_CALL( SCIPaddRealParam(scip, 1402 "heuristics/" HEUR_NAME "/objfactor", 1403 "factor by which the regard of the objective is decreased in each round, 1.0 for dynamic", 1404 &heurdata->objfactor, FALSE, DEFAULT_OBJFACTOR, 0.0, 1.0, NULL, NULL) ); 1405 SCIP_CALL( SCIPaddRealParam(scip, 1406 "heuristics/" HEUR_NAME "/alpha", 1407 "initial weight of the objective function in the convex combination", 1408 &heurdata->alpha, FALSE, DEFAULT_ALPHA, 0.0, 1.0, NULL, NULL) ); 1409 SCIP_CALL( SCIPaddRealParam(scip, 1410 "heuristics/" HEUR_NAME "/alphadiff", 1411 "threshold difference for the convex parameter to perform perturbation", 1412 &heurdata->alphadiff, FALSE, DEFAULT_ALPHADIFF, 0.0, 1.0, NULL, NULL) ); 1413 1414 SCIP_CALL( SCIPaddIntParam(scip, 1415 "heuristics/" HEUR_NAME "/maxlpiterofs", 1416 "additional number of allowed LP iterations", 1417 &heurdata->maxlpiterofs, FALSE, DEFAULT_MAXLPITEROFS, 0, INT_MAX, NULL, NULL) ); 1418 SCIP_CALL( SCIPaddIntParam(scip, 1419 "heuristics/" HEUR_NAME "/maxsols", 1420 "total number of feasible solutions found up to which heuristic is called (-1: no limit)", 1421 &heurdata->maxsols, TRUE, DEFAULT_MAXSOLS, -1, INT_MAX, NULL, NULL) ); 1422 SCIP_CALL( SCIPaddIntParam(scip, 1423 "heuristics/" HEUR_NAME "/maxloops", 1424 "maximal number of pumping loops (-1: no limit)", 1425 &heurdata->maxloops, TRUE, DEFAULT_MAXLOOPS, -1, INT_MAX, NULL, NULL) ); 1426 SCIP_CALL( SCIPaddIntParam(scip, 1427 "heuristics/" HEUR_NAME "/maxstallloops", 1428 "maximal number of pumping rounds without fractionality improvement (-1: no limit)", 1429 &heurdata->maxstallloops, TRUE, DEFAULT_MAXSTALLLOOPS, -1, INT_MAX, NULL, NULL) ); 1430 SCIP_CALL( SCIPaddIntParam(scip, 1431 "heuristics/" HEUR_NAME "/minflips", 1432 "minimum number of random variables to flip, if a 1-cycle is encountered", 1433 &heurdata->minflips, TRUE, DEFAULT_MINFLIPS, 1, INT_MAX, NULL, NULL) ); 1434 SCIP_CALL( SCIPaddIntParam(scip, 1435 "heuristics/" HEUR_NAME "/cyclelength", 1436 "maximum length of cycles to be checked explicitly in each round", 1437 &heurdata->cyclelength, TRUE, DEFAULT_CYCLELENGTH, 1, 100, NULL, NULL) ); 1438 SCIP_CALL( SCIPaddIntParam(scip, 1439 "heuristics/" HEUR_NAME "/perturbfreq", 1440 "number of iterations until a random perturbation is forced", 1441 &heurdata->perturbfreq, TRUE, DEFAULT_PERTURBFREQ, 1, INT_MAX, NULL, NULL) ); 1442 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/neighborhoodsize", 1443 "radius (using Manhattan metric) of the neighborhood to be searched in stage 3", 1444 &heurdata->neighborhoodsize, FALSE, DEFAULT_NEIGHBORHOODSIZE, 1, INT_MAX, NULL, NULL) ); 1445 1446 SCIP_CALL( SCIPaddBoolParam(scip, 1447 "heuristics/" HEUR_NAME "/beforecuts", 1448 "should the feasibility pump be called at root node before cut separation?", 1449 &heurdata->beforecuts, FALSE, DEFAULT_BEFORECUTS, NULL, NULL) ); 1450 SCIP_CALL( SCIPaddBoolParam(scip, 1451 "heuristics/" HEUR_NAME "/usefp20", 1452 "should an iterative round-and-propagate scheme be used to find the integral points?", 1453 &heurdata->usefp20, FALSE, DEFAULT_USEFP20, NULL, NULL) ); 1454 SCIP_CALL( SCIPaddBoolParam(scip, 1455 "heuristics/" HEUR_NAME "/pertsolfound", 1456 "should a random perturbation be performed if a feasible solution was found?", 1457 &heurdata->pertsolfound, FALSE, DEFAULT_PERTSOLFOUND, NULL, NULL) ); 1458 SCIP_CALL( SCIPaddBoolParam(scip, 1459 "heuristics/" HEUR_NAME "/stage3", 1460 "should we solve a local branching sub-MIP if no solution could be found?", 1461 &heurdata->stage3, FALSE, DEFAULT_STAGE3, NULL, NULL) ); 1462 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/copycuts", 1463 "should all active cuts from cutpool be copied to constraints in subproblem?", 1464 &heurdata->copycuts, TRUE, DEFAULT_COPYCUTS, NULL, NULL) ); 1465 1466 return SCIP_OKAY; 1467 } 1468