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 heuristics.c 26 * @ingroup OTHER_CFILES 27 * @brief methods commonly used by primal heuristics 28 * @author Gregor Hendel 29 */ 30 31 #include "scip/heuristics.h" 32 #include "scip/cons_linear.h" 33 #include "scip/scipdefplugins.h" 34 #include "scip/stat.h" 35 #include "scip/struct_scip.h" 36 37 #include "scip/pub_heur.h" 38 39 /* the indicator and SOS1 constraint handlers are included for the diving algorithm SCIPperformGenericDivingAlgorithm() */ 40 #include "scip/cons_indicator.h" 41 #include "scip/cons_sos1.h" 42 43 #define MINLPITER 10000 /**< minimal number of LP iterations allowed in each LP solving call */ 44 45 46 /** solve probing LP */ 47 static 48 SCIP_RETCODE solveLP( 49 SCIP* scip, /**< SCIP data structure */ 50 SCIP_DIVESET* diveset, /**< diving settings */ 51 SCIP_Longint maxnlpiterations, /**< maximum number of allowed LP iterations */ 52 SCIP_DIVECONTEXT divecontext, /**< context for diving statistics */ 53 SCIP_Bool* lperror, /**< pointer to store if an internal LP error occurred */ 54 SCIP_Bool* cutoff /**< pointer to store whether the LP was infeasible */ 55 ) 56 { 57 int lpiterationlimit; 58 SCIP_RETCODE retstat; 59 SCIP_Longint nlpiterations; 60 61 assert(lperror != NULL); 62 assert(cutoff != NULL); 63 64 nlpiterations = SCIPgetNLPIterations(scip); 65 66 /* allow at least MINLPITER more iterations so as not to run out of LP iterations during this solve */ 67 lpiterationlimit = (int)(maxnlpiterations - SCIPdivesetGetNLPIterations(diveset, divecontext)); 68 lpiterationlimit = MAX(lpiterationlimit, MINLPITER); 69 70 retstat = SCIPsolveProbingLP(scip, lpiterationlimit, lperror, cutoff); 71 72 /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic. 73 * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. 74 */ 75 #ifdef NDEBUG 76 if( retstat != SCIP_OKAY ) 77 { 78 SCIPwarningMessage(scip, "Error while solving LP in %s diving heuristic; LP solve terminated with code <%d>.\n", SCIPdivesetGetName(diveset), retstat); 79 } 80 #else 81 SCIP_CALL( retstat ); 82 #endif 83 84 /* update iteration count */ 85 SCIPupdateDivesetLPStats(scip, diveset, SCIPgetNLPIterations(scip) - nlpiterations, divecontext); 86 87 return SCIP_OKAY; 88 } 89 90 /** select the next variable and type of diving */ 91 static 92 SCIP_RETCODE selectNextDiving( 93 SCIP* scip, /**< SCIP data structure */ 94 SCIP_DIVESET* diveset, /**< dive set */ 95 SCIP_SOL* worksol, /**< current working solution */ 96 SCIP_Bool onlylpbranchcands, /**< should only LP branching candidates be considered? */ 97 SCIP_Bool storelpcandscores, /**< should the scores of the LP candidates be updated? */ 98 SCIP_VAR** lpcands, /**< LP branching candidates, or NULL if not needed */ 99 SCIP_Real * lpcandssol, /**< solution values LP branching candidates, or NULL if not needed */ 100 SCIP_Real* lpcandsfrac, /**< fractionalities of LP branching candidates, or NULL if not needed*/ 101 SCIP_Real* lpcandsscores, /**< array with LP branching candidate scores, or NULL */ 102 SCIP_Bool* lpcandroundup, /**< array to remember whether the preferred branching direction is upwards */ 103 int* nviollpcands, /**< pointer to store the number of LP candidates whose solution value already violates local bounds */ 104 int nlpcands, /**< number of current LP cands */ 105 SCIP_Bool* enfosuccess, /**< pointer to store whether a candidate was sucessfully found */ 106 SCIP_Bool* infeasible /**< pointer to store whether the diving can be immediately aborted because it is infeasible */ 107 ) 108 { 109 assert(scip != NULL); 110 assert(worksol != NULL); 111 assert(!onlylpbranchcands || lpcandsscores != NULL); 112 assert(!onlylpbranchcands || lpcandroundup != NULL); 113 assert(enfosuccess != NULL); 114 assert(infeasible != NULL); 115 assert(nviollpcands != NULL); 116 117 *nviollpcands = 0; 118 /* we use diving solution enforcement provided by the constraint handlers */ 119 if( !onlylpbranchcands ) 120 { 121 SCIP_CALL( SCIPgetDiveBoundChanges(scip, diveset, worksol, enfosuccess, infeasible) ); 122 } 123 else 124 { 125 int c; 126 int bestcandidx; 127 SCIP_Real bestscore; 128 SCIP_Real score; 129 130 bestscore = SCIP_REAL_MIN; 131 bestcandidx = -1; 132 133 SCIPclearDiveBoundChanges(scip); 134 135 /* search for the candidate that maximizes the dive set score function and whose solution value is still feasible */ 136 for( c = 0; c < nlpcands; ++c ) 137 { 138 assert(SCIPgetSolVal(scip, worksol, lpcands[c]) == lpcandssol[c]); /*lint !e777 doesn't like comparing floats for equality */ 139 140 /* scores are kept in arrays for faster reuse */ 141 if( storelpcandscores ) 142 { 143 SCIP_CALL( SCIPgetDivesetScore(scip, diveset, SCIP_DIVETYPE_INTEGRALITY, lpcands[c], lpcandssol[c], 144 lpcandsfrac[c], &lpcandsscores[c], &lpcandroundup[c]) ); 145 } 146 147 score = lpcandsscores[c]; 148 /* update the best candidate if it has a higher score and a solution value which does not violate one of the local bounds */ 149 if( SCIPisLE(scip, SCIPvarGetLbLocal(lpcands[c]), lpcandssol[c]) && SCIPisGE(scip, SCIPvarGetUbLocal(lpcands[c]), lpcandssol[c]) ) 150 { 151 if( score > bestscore ) 152 { 153 bestcandidx = c; 154 bestscore = score; 155 } 156 } 157 else 158 ++(*nviollpcands); 159 } 160 161 /* there is no guarantee that a candidate is found since local bounds might render all solution values infeasible */ 162 *enfosuccess = (bestcandidx >= 0); 163 if( *enfosuccess ) 164 { 165 /* if we want to round up the best candidate, it is added as the preferred bound change */ 166 SCIP_CALL( SCIPaddDiveBoundChange(scip, lpcands[bestcandidx], SCIP_BRANCHDIR_UPWARDS, 167 SCIPceil(scip, lpcandssol[bestcandidx]), lpcandroundup[bestcandidx]) ); 168 SCIP_CALL( SCIPaddDiveBoundChange(scip, lpcands[bestcandidx], SCIP_BRANCHDIR_DOWNWARDS, 169 SCIPfloor(scip, lpcandssol[bestcandidx]), ! lpcandroundup[bestcandidx]) ); 170 } 171 } 172 return SCIP_OKAY; 173 } 174 175 /** return the LP iteration budget suggestion for this dive set */ 176 static 177 SCIP_Longint getDivesetIterLimit( 178 SCIP* scip, /**< SCIP data structure */ 179 SCIP_DIVESET* diveset, /**< dive set data structure */ 180 SCIP_DIVECONTEXT divecontext /**< context for diving statistics */ 181 ) 182 { 183 SCIP_Longint iterlimit; 184 /*todo another factor of 10, REALLY? */ 185 iterlimit = (SCIP_Longint)((1.0 + 10*(SCIPdivesetGetNSols(diveset, divecontext)+1.0)/(SCIPdivesetGetNCalls(diveset, divecontext)+1.0)) * SCIPdivesetGetMaxLPIterQuot(diveset) * SCIPgetNNodeLPIterations(scip)); 186 iterlimit += SCIPdivesetGetMaxLPIterOffset(diveset); 187 iterlimit -= SCIPdivesetGetNLPIterations(diveset, divecontext); 188 189 return iterlimit; 190 } 191 192 /** performs a diving within the limits of the @p diveset parameters 193 * 194 * This method performs a diving according to the settings defined by the diving settings @p diveset; Contrary to the 195 * name, SCIP enters probing mode (not diving mode) and dives along a path into the tree. Domain propagation 196 * is applied at every node in the tree, whereas probing LPs might be solved less frequently. 197 * 198 * Starting from the current LP solution, the algorithm selects candidates which maximize the 199 * score defined by the @p diveset and whose solution value has not yet been rendered infeasible by propagation, 200 * and propagates the bound change on this candidate. 201 * 202 * The algorithm iteratively selects the the next (unfixed) candidate in the list, until either enough domain changes 203 * or the resolve frequency of the LP trigger an LP resolve (and hence, the set of potential candidates changes), 204 * or the last node is proven to be infeasible. It optionally backtracks and tries the 205 * other branching direction. 206 * 207 * After the set of remaining candidates is empty or the targeted depth is reached, the node LP is 208 * solved, and the old candidates are replaced by the new LP candidates. 209 * 210 * @see heur_guideddiving.c for an example implementation of a dive set controlling the diving algorithm. 211 * 212 * @note the node from where the algorithm is called is checked for a basic LP solution. If the solution 213 * is non-basic, e.g., when barrier without crossover is used, the method returns without performing a dive. 214 * 215 * @note currently, when multiple diving heuristics call this method and solve an LP at the same node, only the first 216 * call will be executed, see SCIPgetLastDiveNode() 217 * 218 * @todo generalize method to work correctly with pseudo or external branching/diving candidates 219 */ 220 SCIP_RETCODE SCIPperformGenericDivingAlgorithm( 221 SCIP* scip, /**< SCIP data structure */ 222 SCIP_DIVESET* diveset, /**< settings for diving */ 223 SCIP_SOL* worksol, /**< non-NULL working solution */ 224 SCIP_HEUR* heur, /**< the calling primal heuristic */ 225 SCIP_RESULT* result, /**< SCIP result pointer */ 226 SCIP_Bool nodeinfeasible, /**< is the current node known to be infeasible? */ 227 SCIP_Longint iterlim, /**< nonnegative iteration limit for the LP solves, or -1 for dynamic setting */ 228 int nodelimit, /**< nonnegative probing node limit or -1 if no limit should be used */ 229 SCIP_Real lpresolvedomchgquot, /**< percentage of immediate domain changes during probing to trigger LP resolve or -1 230 if diveset specific default should be used */ 231 SCIP_DIVECONTEXT divecontext /**< context for diving statistics */ 232 ) 233 { 234 SCIP_CONSHDLR* indconshdlr; /* constraint handler for indicator constraints */ 235 SCIP_CONSHDLR* sos1conshdlr; /* constraint handler for SOS1 constraints */ 236 SCIP_VAR** lpcands; 237 SCIP_Real* lpcandssol; 238 239 SCIP_VAR** previouscands; 240 SCIP_Real* lpcandsscores; 241 SCIP_Real* previousvals; 242 SCIP_Real* lpcandsfrac; 243 SCIP_Bool* lpcandroundup; 244 SCIP_Real searchubbound; 245 SCIP_Real searchavgbound; 246 SCIP_Real searchbound; 247 SCIP_Real ubquot; 248 SCIP_Real avgquot; 249 SCIP_Longint maxnlpiterations; 250 SCIP_Longint domreds; 251 int startndivecands; 252 int depth; 253 int maxdepth; 254 int maxdivedepth; 255 int totalnbacktracks; 256 int totalnprobingnodes; 257 int lastlpdepth; 258 int previouscandssize; 259 int lpcandsscoressize; 260 int nviollpcands; 261 SCIP_Longint oldnsolsfound; 262 SCIP_Longint oldnbestsolsfound; 263 SCIP_Longint oldnconflictsfound; 264 265 SCIP_Bool success; 266 SCIP_Bool leafsol; 267 SCIP_Bool enfosuccess; 268 SCIP_Bool lperror; 269 SCIP_Bool cutoff; 270 SCIP_Bool backtracked; 271 SCIP_Bool backtrack; 272 SCIP_Bool onlylpbranchcands; 273 274 int nlpcands; 275 int lpsolvefreq; 276 277 assert(scip != NULL); 278 assert(heur != NULL); 279 assert(result != NULL); 280 assert(worksol != NULL); 281 282 *result = SCIP_DELAYED; 283 284 /* do not call heuristic in node that was already detected to be infeasible */ 285 if( nodeinfeasible ) 286 return SCIP_OKAY; 287 288 /* only call heuristic, if an optimal LP solution is at hand */ 289 if( !SCIPhasCurrentNodeLP(scip) || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 290 return SCIP_OKAY; 291 292 /* only call heuristic, if the LP solution is basic (which allows fast resolve in diving) */ 293 if( !SCIPisLPSolBasic(scip) ) 294 return SCIP_OKAY; 295 296 /* when heuristic was called by scheduler, be less restrictive */ 297 if( divecontext != SCIP_DIVECONTEXT_SCHEDULER ) 298 { 299 /* don't dive two times at the same node */ 300 if( SCIPgetLastDivenode(scip) == SCIPgetNNodes(scip) && SCIPgetDepth(scip) > 0 ) 301 return SCIP_OKAY; 302 303 /* only call heuristic, if the LP objective value is smaller than the cutoff bound */ 304 if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) ) 305 return SCIP_OKAY; 306 } 307 308 *result = SCIP_DIDNOTRUN; 309 310 /* only try to dive, if we are in the correct part of the tree, given by minreldepth and maxreldepth */ 311 depth = SCIPgetDepth(scip); 312 maxdepth = SCIPgetMaxDepth(scip); 313 maxdepth = MAX(maxdepth, 30); 314 if( divecontext != SCIP_DIVECONTEXT_SCHEDULER && 315 (depth < SCIPdivesetGetMinRelDepth(diveset) * maxdepth || depth > SCIPdivesetGetMaxRelDepth(diveset) * maxdepth) ) 316 { 317 return SCIP_OKAY; 318 } 319 320 /* calculate the maximal number of LP iterations */ 321 if( iterlim < 0 ) 322 { 323 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + getDivesetIterLimit(scip, diveset, divecontext); 324 } 325 else 326 { 327 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + iterlim; 328 } 329 330 /* don't try to dive, if we took too many LP iterations during diving */ 331 if( SCIPdivesetGetNLPIterations(diveset, divecontext) >= maxnlpiterations ) 332 return SCIP_OKAY; 333 334 /* allow at least a certain number of LP iterations in this dive */ 335 if( SCIPdivesetGetNLPIterations(diveset, divecontext) + MINLPITER > maxnlpiterations ) 336 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + MINLPITER; 337 338 /* these constraint handlers are required for polishing an LP relaxation solution beyond rounding */ 339 indconshdlr = SCIPfindConshdlr(scip, "indicator"); 340 sos1conshdlr = SCIPfindConshdlr(scip, "SOS1"); 341 342 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, &lpcandsfrac, &nlpcands, NULL, NULL) ); 343 344 onlylpbranchcands = SCIPdivesetUseOnlyLPBranchcands(diveset); 345 /* don't try to dive, if there are no diving candidates */ 346 if( onlylpbranchcands && nlpcands == 0 ) 347 return SCIP_OKAY; 348 349 /* calculate the objective search bound */ 350 if( SCIPgetNSolsFound(scip) == 0 ) 351 { 352 ubquot = SCIPdivesetGetUbQuotNoSol(diveset); 353 avgquot = SCIPdivesetGetAvgQuotNoSol(diveset); 354 } 355 else 356 { 357 ubquot = SCIPdivesetGetUbQuot(diveset); 358 avgquot = SCIPdivesetGetAvgQuot(diveset); 359 } 360 361 if( ubquot > 0.0 ) 362 searchubbound = SCIPgetLowerbound(scip) + ubquot * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip)); 363 else 364 searchubbound = SCIPinfinity(scip); 365 366 if( avgquot > 0.0 ) 367 searchavgbound = SCIPgetLowerbound(scip) + avgquot * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip)); 368 else 369 searchavgbound = SCIPinfinity(scip); 370 371 searchbound = MIN(searchubbound, searchavgbound); 372 373 if( SCIPisObjIntegral(scip) ) 374 searchbound = SCIPceil(scip, searchbound); 375 376 /* calculate the maximal diving depth: 10 * min{number of integer variables, max depth}*/ 377 maxdivedepth = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip); 378 if ( sos1conshdlr != NULL ) 379 maxdivedepth += SCIPgetNSOS1Vars(sos1conshdlr); 380 maxdivedepth = MIN(maxdivedepth, maxdepth); 381 maxdivedepth *= 10; 382 383 /* if lpresolvedomchgquot is not explicitly given (so -1), then use the default for the current diveset */ 384 if( lpresolvedomchgquot < 0 ) 385 lpresolvedomchgquot = SCIPdivesetGetLPResolveDomChgQuot(diveset); 386 assert(lpresolvedomchgquot > 0.0 && lpresolvedomchgquot <= 1.0); 387 388 *result = SCIP_DIDNOTFIND; 389 390 /* start probing mode */ 391 SCIP_CALL( SCIPstartProbing(scip) ); 392 393 /* enables collection of variable statistics during probing */ 394 SCIPenableVarHistory(scip); 395 396 SCIPdebugMsg(scip, "(node %" SCIP_LONGINT_FORMAT ") executing %s heuristic: depth=%d, %d fractionals, dualbound=%g, avgbound=%g, cutoffbound=%g, searchbound=%g\n", 397 SCIPgetNNodes(scip), SCIPheurGetName(heur), SCIPgetDepth(scip), nlpcands, SCIPgetDualbound(scip), SCIPgetAvgDualbound(scip), 398 SCIPretransformObj(scip, SCIPgetCutoffbound(scip)), SCIPretransformObj(scip, searchbound)); 399 400 /* allocate buffer storage for previous candidates and their branching values for pseudo cost updates */ 401 lpsolvefreq = SCIPdivesetGetLPSolveFreq(diveset); 402 previouscandssize = MAX(1, lpsolvefreq); 403 SCIP_CALL( SCIPallocBufferArray(scip, &previouscands, previouscandssize) ); 404 SCIP_CALL( SCIPallocBufferArray(scip, &previousvals, previouscandssize) ); 405 406 /* keep some statistics */ 407 lperror = FALSE; 408 cutoff = FALSE; 409 lastlpdepth = -1; 410 startndivecands = nlpcands; 411 totalnbacktracks = 0; 412 totalnprobingnodes = 0; 413 oldnsolsfound = SCIPgetNSolsFound(scip); 414 oldnbestsolsfound = SCIPgetNBestSolsFound(scip); 415 oldnconflictsfound = SCIPgetNConflictConssFound(scip); 416 417 /* link the working solution to the dive set */ 418 SCIPdivesetSetWorkSolution(diveset, worksol); 419 420 if( onlylpbranchcands ) 421 { 422 SCIP_CALL( SCIPallocBufferArray(scip, &lpcandsscores, nlpcands) ); 423 SCIP_CALL( SCIPallocBufferArray(scip, &lpcandroundup, nlpcands) ); 424 425 lpcandsscoressize = nlpcands; 426 } 427 else 428 { 429 lpcandroundup = NULL; 430 lpcandsscores = NULL; 431 lpcandsscoressize = 0; 432 } 433 434 enfosuccess = TRUE; 435 leafsol = FALSE; 436 437 /* LP loop; every time a new LP was solved, conditions are checked 438 * dive as long we are in the given objective, depth and iteration limits and fractional variables exist, but 439 * - if possible, we dive at least with the depth 10 440 * - if the number of fractional variables decreased at least with 1 variable per 2 dive depths, we continue diving 441 */ 442 while( !lperror && !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL && enfosuccess 443 && ( nodelimit == -1 || totalnprobingnodes < nodelimit ) 444 && (SCIPgetProbingDepth(scip) < 10 445 || nlpcands <= startndivecands - SCIPgetProbingDepth(scip) / 2 446 || (SCIPgetProbingDepth(scip) < maxdivedepth && SCIPdivesetGetNLPIterations(diveset, divecontext) < maxnlpiterations && SCIPgetLPObjval(scip) < searchbound)) 447 && !SCIPisStopped(scip) ) 448 { 449 SCIP_Real lastlpobjval; 450 int c; 451 SCIP_Bool allroundable; 452 SCIP_Bool infeasible; 453 int prevcandsinsertidx; 454 455 /* remember the last LP depth */ 456 assert(lastlpdepth < SCIPgetProbingDepth(scip)); 457 lastlpdepth = SCIPgetProbingDepth(scip); 458 domreds = 0; 459 460 SCIPdebugMsg(scip, "%s heuristic continues diving at depth %d, %d candidates left\n", 461 SCIPdivesetGetName(diveset), lastlpdepth, nlpcands); 462 463 /* loop over candidates and determine if they are roundable */ 464 allroundable = TRUE; 465 c = 0; 466 while( allroundable && c < nlpcands ) 467 { 468 if( SCIPvarMayRoundDown(lpcands[c]) || SCIPvarMayRoundUp(lpcands[c]) ) 469 allroundable = TRUE; 470 else 471 allroundable = FALSE; 472 ++c; 473 } 474 475 /* if all candidates are roundable, try to round the solution */ 476 if( allroundable ) 477 { 478 success = FALSE; 479 480 /* working solution must be linked to LP solution */ 481 SCIP_CALL( SCIPlinkLPSol(scip, worksol) ); 482 /* create solution from diving LP and try to round it */ 483 SCIP_CALL( SCIProundSol(scip, worksol, &success) ); 484 485 /* adjust SOS1 constraints */ 486 if( success && sos1conshdlr != NULL ) 487 { 488 SCIP_Bool changed; 489 SCIP_CALL( SCIPmakeSOS1sFeasible(scip, sos1conshdlr, worksol, &changed, &success) ); 490 } 491 492 /* succesfully rounded solutions are tried for primal feasibility */ 493 if( success ) 494 { 495 SCIP_Bool changed = FALSE; 496 SCIPdebugMsg(scip, "%s found roundable primal solution: obj=%g\n", SCIPdivesetGetName(diveset), SCIPgetSolOrigObj(scip, worksol)); 497 498 /* adjust indicator constraints */ 499 if( indconshdlr != NULL ) 500 { 501 SCIP_CALL( SCIPmakeIndicatorsFeasible(scip, indconshdlr, worksol, &changed) ); 502 } 503 504 success = FALSE; 505 506 /* try to add solution to SCIP */ 507 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, changed, &success) ); 508 509 /* check, if solution was feasible and good enough */ 510 if( success ) 511 { 512 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n"); 513 *result = SCIP_FOUNDSOL; 514 leafsol = (nlpcands == 0); 515 516 /* the rounded solution found above led to a cutoff of the node LP solution */ 517 if( SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OBJLIMIT ) 518 { 519 cutoff = TRUE; 520 break; 521 } 522 } 523 } 524 } 525 526 /* working solution must be linked to LP solution */ 527 assert(SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL); 528 lastlpobjval = SCIPgetLPObjval(scip); 529 SCIP_CALL( SCIPlinkLPSol(scip, worksol) ); 530 531 /* we must explicitly store the solution values by unlinking the solution, otherwise, 532 * the working solution may contain wrong entries, if, e.g., a backtrack occurred after an 533 * intermediate LP resolve or the LP was resolved during conflict analysis. 534 */ 535 SCIP_CALL( SCIPunlinkSol(scip, worksol) ); 536 537 /* ensure array sizes for the diving on the fractional variables */ 538 if( onlylpbranchcands && nlpcands > lpcandsscoressize ) 539 { 540 assert(nlpcands > 0); 541 assert(lpcandsscores != NULL); 542 assert(lpcandroundup != NULL); 543 544 SCIP_CALL( SCIPreallocBufferArray(scip, &lpcandsscores, nlpcands) ); 545 SCIP_CALL( SCIPreallocBufferArray(scip, &lpcandroundup, nlpcands) ); 546 547 lpcandsscoressize = nlpcands; 548 } 549 550 enfosuccess = FALSE; 551 /* select the next diving action by selecting appropriate dive bound changes for the preferred and alternative child */ 552 SCIP_CALL( selectNextDiving(scip, diveset, worksol, onlylpbranchcands, SCIPgetProbingDepth(scip) == lastlpdepth, 553 lpcands, lpcandssol, lpcandsfrac, lpcandsscores, lpcandroundup, &nviollpcands, nlpcands, 554 &enfosuccess, &infeasible) ); 555 556 /* if we did not succeed finding an enforcement, the solution is potentially feasible and we break immediately */ 557 if( ! enfosuccess ) 558 break; 559 560 prevcandsinsertidx = -1; 561 562 /* start propagating candidate variables 563 * - until the desired targetdepth is reached, 564 * - or there is no further candidate variable left because of intermediate bound changes, 565 * - or a cutoff is detected 566 */ 567 do 568 { 569 SCIP_VAR* bdchgvar; 570 SCIP_Real bdchgvalue; 571 SCIP_Longint localdomreds; 572 SCIP_BRANCHDIR bdchgdir; 573 int nbdchanges; 574 575 /* ensure that a new candidate was successfully determined (usually at the end of the previous loop iteration) */ 576 assert(enfosuccess); 577 bdchgvar = NULL; 578 bdchgvalue = SCIP_INVALID; 579 bdchgdir = SCIP_BRANCHDIR_AUTO; 580 581 backtracked = FALSE; 582 do 583 { 584 int d; 585 SCIP_VAR** bdchgvars; 586 SCIP_BRANCHDIR* bdchgdirs; 587 SCIP_Real* bdchgvals; 588 589 nbdchanges = 0; 590 /* get the bound change information stored in the dive set */ 591 SCIPgetDiveBoundChangeData(scip, &bdchgvars, &bdchgdirs, &bdchgvals, &nbdchanges, !backtracked); 592 593 assert(nbdchanges > 0); 594 assert(bdchgvars != NULL); 595 assert(bdchgdirs != NULL); 596 assert(bdchgvals != NULL); 597 598 /* return if we reached the depth limit */ 599 if( SCIP_MAXTREEDEPTH <= SCIPgetDepth(scip) ) 600 { 601 SCIPdebugMsg(scip, "depth limit reached, we stop diving immediately.\n"); 602 goto TERMINATE; 603 } 604 605 /* return if we reached the node limit */ 606 if( nodelimit != -1 && totalnprobingnodes >= nodelimit ) 607 { 608 SCIPdebugMsg(scip, "node limit reached, we stop diving immediately.\n"); 609 goto TERMINATE; 610 } 611 612 /* dive deeper into the tree */ 613 SCIP_CALL( SCIPnewProbingNode(scip) ); 614 ++totalnprobingnodes; 615 616 /* apply all suggested domain changes of the variables */ 617 for( d = 0; d < nbdchanges; ++d ) 618 { 619 SCIP_Real lblocal; 620 SCIP_Real ublocal; 621 SCIP_Bool infeasbdchange; 622 623 /* get the next bound change data: variable, direction, and value */ 624 bdchgvar = bdchgvars[d]; 625 bdchgvalue = bdchgvals[d]; 626 bdchgdir = bdchgdirs[d]; 627 628 assert(bdchgvar != NULL); 629 lblocal = SCIPvarGetLbLocal(bdchgvar); 630 ublocal = SCIPvarGetUbLocal(bdchgvar); 631 632 SCIPdebugMsg(scip, " dive %d/%d, LP iter %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ": var <%s>, oldbounds=[%g,%g], newbounds=[%g,%g]\n", 633 SCIPgetProbingDepth(scip), maxdivedepth, SCIPdivesetGetNLPIterations(diveset, divecontext), maxnlpiterations, 634 SCIPvarGetName(bdchgvar), 635 lblocal, ublocal, 636 bdchgdir == SCIP_BRANCHDIR_DOWNWARDS ? lblocal : bdchgvalue, 637 bdchgdir == SCIP_BRANCHDIR_UPWARDS ? ublocal : bdchgvalue); 638 639 infeasbdchange = FALSE; 640 /* tighten the lower and/or upper bound depending on the bound change type */ 641 switch( bdchgdir ) 642 { 643 case SCIP_BRANCHDIR_UPWARDS: 644 /* test if bound change is possible, otherwise propagation might have deduced the same 645 * bound already or numerical troubles might have occurred */ 646 if( SCIPisFeasLE(scip, bdchgvalue, lblocal) || SCIPisFeasGT(scip, bdchgvalue, ublocal) ) 647 infeasbdchange = TRUE; 648 else 649 { 650 /* round variable up */ 651 SCIP_CALL( SCIPchgVarLbProbing(scip, bdchgvar, bdchgvalue) ); 652 } 653 break; 654 case SCIP_BRANCHDIR_DOWNWARDS: 655 /* test if bound change is possible, otherwise propagation might have deduced the same 656 * bound already or numerical troubles might have occurred */ 657 if( SCIPisFeasGE(scip, bdchgvalue, ublocal) || SCIPisFeasLT(scip, bdchgvalue, lblocal) ) 658 infeasbdchange = TRUE; 659 else 660 { 661 /* round variable down */ 662 SCIP_CALL( SCIPchgVarUbProbing(scip, bdchgvar, bdchgvalue) ); 663 } 664 break; 665 case SCIP_BRANCHDIR_FIXED: 666 /* test if bound change is possible, otherwise propagation might have deduced the same 667 * bound already or numerical troubles might have occurred */ 668 if( SCIPisFeasLT(scip, bdchgvalue, lblocal) || SCIPisFeasGT(scip, bdchgvalue, ublocal) || 669 (SCIPisFeasEQ(scip, lblocal, ublocal) && nbdchanges < 2) ) 670 infeasbdchange = TRUE; 671 else 672 { 673 /* fix variable to the bound change value */ 674 if( SCIPisFeasLT(scip, lblocal, bdchgvalue) ) 675 { 676 SCIP_CALL( SCIPchgVarLbProbing(scip, bdchgvar, bdchgvalue) ); 677 } 678 if( SCIPisFeasGT(scip, ublocal, bdchgvalue) ) 679 { 680 SCIP_CALL( SCIPchgVarUbProbing(scip, bdchgvar, bdchgvalue) ); 681 } 682 } 683 break; 684 case SCIP_BRANCHDIR_AUTO: 685 default: 686 SCIPerrorMessage("Error: Unsupported bound change direction <%d> specified for diving, aborting\n",bdchgdirs[d]); 687 SCIPABORT(); 688 return SCIP_INVALIDDATA; /*lint !e527*/ 689 } 690 /* if the variable domain has been shrunk in the meantime, numerical troubles may have 691 * occured or variable was fixed by propagation while backtracking => Abort diving! 692 */ 693 if( infeasbdchange ) 694 { 695 SCIPdebugMsg(scip, "\nSelected variable <%s> domain already [%g,%g] as least as tight as desired bound change, diving aborted \n", 696 SCIPvarGetName(bdchgvar), SCIPvarGetLbLocal(bdchgvar), SCIPvarGetUbLocal(bdchgvar)); 697 cutoff = TRUE; 698 break; 699 } 700 } 701 /* break loop immediately if we detected a cutoff */ 702 if( cutoff ) 703 break; 704 705 /* apply domain propagation */ 706 localdomreds = 0; 707 SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, &localdomreds) ); 708 709 /* add the number of bound changes we applied by ourselves after propagation, otherwise the counter would have been reset */ 710 localdomreds += nbdchanges; 711 712 SCIPdebugMsg(scip, "domain reductions: %" SCIP_LONGINT_FORMAT " (total: %" SCIP_LONGINT_FORMAT ")\n", 713 localdomreds, domreds + localdomreds); 714 715 /* resolve the diving LP if the diving resolve frequency is reached or a sufficient number of intermediate bound changes 716 * was reached 717 */ 718 if( ! cutoff 719 && ((lpsolvefreq > 0 && ((SCIPgetProbingDepth(scip) - lastlpdepth) % lpsolvefreq) == 0) 720 || (domreds + localdomreds > SCIPdivesetGetLPResolveDomChgQuot(diveset) * SCIPgetNVars(scip)) 721 || (onlylpbranchcands && nviollpcands > (int)(SCIPdivesetGetLPResolveDomChgQuot(diveset) * nlpcands))) ) 722 { 723 SCIP_CALL( solveLP(scip, diveset, maxnlpiterations, divecontext, &lperror, &cutoff) ); 724 725 /* lp errors lead to early termination */ 726 if( lperror ) 727 { 728 cutoff = TRUE; 729 break; 730 } 731 } 732 733 /* perform backtracking if a cutoff was detected */ 734 if( cutoff && !backtracked && SCIPdivesetUseBacktrack(diveset) ) 735 { 736 SCIPdebugMsg(scip, " *** cutoff detected at level %d - backtracking\n", SCIPgetProbingDepth(scip)); 737 SCIP_CALL( SCIPbacktrackProbing(scip, SCIPgetProbingDepth(scip) - 1) ); 738 ++totalnbacktracks; 739 backtracked = TRUE; 740 backtrack = TRUE; 741 cutoff = FALSE; 742 } 743 else 744 backtrack = FALSE; 745 } 746 while( backtrack ); 747 748 /* we add the domain reductions from the last evaluated node */ 749 domreds += localdomreds; /*lint !e771 lint thinks localdomreds has not been initialized */ 750 751 /* store candidate for pseudo cost update and choose next candidate only if no cutoff was detected */ 752 if( ! cutoff ) 753 { 754 if( nbdchanges == 1 && (bdchgdir == SCIP_BRANCHDIR_UPWARDS || bdchgdir == SCIP_BRANCHDIR_DOWNWARDS) ) 755 { 756 ++prevcandsinsertidx; 757 assert(prevcandsinsertidx <= SCIPgetProbingDepth(scip) - lastlpdepth - 1); 758 assert(SCIPgetProbingDepth(scip) > 0); 759 assert(bdchgvar != NULL); 760 assert(bdchgvalue != SCIP_INVALID); /*lint !e777 doesn't like comparing floats for equality */ 761 762 /* extend array in case of a dynamic, domain change based LP resolve strategy */ 763 if( prevcandsinsertidx >= previouscandssize ) 764 { 765 previouscandssize *= 2; 766 SCIP_CALL( SCIPreallocBufferArray(scip, &previouscands, previouscandssize) ); 767 SCIP_CALL( SCIPreallocBufferArray(scip, &previousvals, previouscandssize) ); 768 } 769 assert(previouscandssize > prevcandsinsertidx); 770 771 /* store candidate for pseudo cost update */ 772 previouscands[prevcandsinsertidx] = bdchgvar; 773 previousvals[prevcandsinsertidx] = bdchgvalue; 774 } 775 776 /* choose next candidate variable and resolve the LP if none is found. */ 777 if( SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_NOTSOLVED ) 778 { 779 assert(SCIPgetProbingDepth(scip) > lastlpdepth); 780 enfosuccess = FALSE; 781 782 /* select the next diving action */ 783 SCIP_CALL( selectNextDiving(scip, diveset, worksol, onlylpbranchcands, SCIPgetProbingDepth(scip) == lastlpdepth, 784 lpcands, lpcandssol, lpcandsfrac, lpcandsscores, lpcandroundup, &nviollpcands, nlpcands, 785 &enfosuccess, &infeasible) ); 786 787 /* in case of an unsuccesful candidate search, we solve the node LP */ 788 if( !enfosuccess ) 789 { 790 SCIP_CALL( solveLP(scip, diveset, maxnlpiterations, divecontext, &lperror, &cutoff) ); 791 792 /* check for an LP error and terminate in this case, cutoffs lead to termination anyway */ 793 if( lperror ) 794 cutoff = TRUE; 795 796 /* enfosuccess must be set to TRUE for entering the main LP loop again */ 797 enfosuccess = TRUE; 798 } 799 } 800 } 801 } 802 while( !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_NOTSOLVED ); 803 804 assert(cutoff || lperror || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_NOTSOLVED); 805 806 assert(cutoff || (SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OBJLIMIT && SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_INFEASIBLE && 807 (SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL || SCIPisLT(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip))))); 808 809 /* check new LP candidates and use the LP Objective gain to update pseudo cost information */ 810 if( ! cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ) 811 { 812 int v; 813 SCIP_Real gain; 814 815 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL) ); 816 817 SCIPdebugMsg(scip, " -> lpsolstat=%d, objval=%g/%g, nfrac=%d\n", SCIPgetLPSolstat(scip), SCIPgetLPObjval(scip), searchbound, nlpcands); 818 /* distribute the gain equally over all variables that we rounded since the last LP */ 819 gain = SCIPgetLPObjval(scip) - lastlpobjval; 820 gain = MAX(gain, 0.0); 821 gain /= (1.0 * (SCIPgetProbingDepth(scip) - lastlpdepth)); 822 823 /* loop over previously fixed candidates and share gain improvement */ 824 for( v = 0; v <= prevcandsinsertidx; ++v ) 825 { 826 SCIP_VAR* cand = previouscands[v]; 827 SCIP_Real val = previousvals[v]; 828 SCIP_Real solval = SCIPgetSolVal(scip, worksol, cand); 829 830 /* todo: should the gain be shared with a smaller weight, instead of dividing the gain itself? */ 831 /* it may happen that a variable had an integral solution value beforehand, e.g., for indicator variables */ 832 if( ! SCIPisZero(scip, val - solval) ) 833 { 834 SCIP_CALL( SCIPupdateVarPseudocost(scip, cand, val - solval, gain, 1.0) ); 835 } 836 } 837 } 838 else 839 nlpcands = 0; 840 } 841 842 success = FALSE; 843 /* check if a solution has been found */ 844 if( !enfosuccess && !lperror && !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ) 845 { 846 /* create solution from diving LP */ 847 SCIP_CALL( SCIPlinkLPSol(scip, worksol) ); 848 SCIPdebugMsg(scip, "%s found primal solution: obj=%g\n", SCIPdivesetGetName(diveset), SCIPgetSolOrigObj(scip, worksol)); 849 850 /* try to add solution to SCIP */ 851 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) ); 852 853 /* check, if solution was feasible and good enough */ 854 if( success ) 855 { 856 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n"); 857 *result = SCIP_FOUNDSOL; 858 leafsol = TRUE; 859 } 860 } 861 862 SCIPupdateDivesetStats(scip, diveset, totalnprobingnodes, totalnbacktracks, SCIPgetNSolsFound(scip) - oldnsolsfound, 863 SCIPgetNBestSolsFound(scip) - oldnbestsolsfound, SCIPgetNConflictConssFound(scip) - oldnconflictsfound, leafsol, divecontext); 864 865 SCIPdebugMsg(scip, "(node %" SCIP_LONGINT_FORMAT ") finished %s diveset (%s heur): %d fractionals, dive %d/%d, LP iter %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ", objval=%g/%g, lpsolstat=%d, cutoff=%u\n", 866 SCIPgetNNodes(scip), SCIPdivesetGetName(diveset), SCIPheurGetName(heur), nlpcands, SCIPgetProbingDepth(scip), maxdivedepth, SCIPdivesetGetNLPIterations(diveset, divecontext), maxnlpiterations, 867 SCIPretransformObj(scip, SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ? SCIPgetLPObjval(scip) : SCIPinfinity(scip)), SCIPretransformObj(scip, searchbound), SCIPgetLPSolstat(scip), cutoff); 868 869 TERMINATE: 870 /* end probing mode */ 871 SCIP_CALL( SCIPendProbing(scip) ); 872 873 SCIPdivesetSetWorkSolution(diveset, NULL); 874 875 if( onlylpbranchcands ) 876 { 877 SCIPfreeBufferArray(scip, &lpcandroundup); 878 SCIPfreeBufferArray(scip, &lpcandsscores); 879 } 880 SCIPfreeBufferArray(scip, &previousvals); 881 SCIPfreeBufferArray(scip, &previouscands); 882 883 return SCIP_OKAY; 884 } 885 886 887 /** creates the rows of the subproblem */ 888 static 889 SCIP_RETCODE createRows( 890 SCIP* scip, /**< original SCIP data structure */ 891 SCIP* subscip, /**< SCIP data structure for the subproblem */ 892 SCIP_HASHMAP* varmap /**< a hashmap to store the mapping of source variables to the corresponding 893 * target variables */ 894 ) 895 { 896 SCIP_ROW** rows; /* original scip rows */ 897 SCIP_CONS* cons; /* new constraint */ 898 SCIP_VAR** consvars; /* new constraint's variables */ 899 SCIP_COL** cols; /* original row's columns */ 900 901 SCIP_Real constant; /* constant added to the row */ 902 SCIP_Real lhs; /* left hand side of the row */ 903 SCIP_Real rhs; /* left right side of the row */ 904 SCIP_Real* vals; /* variables' coefficient values of the row */ 905 906 int nrows; 907 int nnonz; 908 int i; 909 int j; 910 911 /* get the rows and their number */ 912 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) ); 913 914 /* copy all rows to linear constraints */ 915 for( i = 0; i < nrows; i++ ) 916 { 917 /* ignore rows that are only locally valid */ 918 if( SCIProwIsLocal(rows[i]) ) 919 continue; 920 921 /* get the row's data */ 922 constant = SCIProwGetConstant(rows[i]); 923 lhs = SCIProwGetLhs(rows[i]) - constant; 924 rhs = SCIProwGetRhs(rows[i]) - constant; 925 vals = SCIProwGetVals(rows[i]); 926 nnonz = SCIProwGetNNonz(rows[i]); 927 cols = SCIProwGetCols(rows[i]); 928 929 assert(lhs <= rhs); 930 931 /* allocate memory array to be filled with the corresponding subproblem variables */ 932 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nnonz) ); 933 for( j = 0; j < nnonz; j++ ) 934 consvars[j] = (SCIP_VAR*) SCIPhashmapGetImage(varmap, (SCIPcolGetVar(cols[j]))); 935 936 /* create a new linear constraint and add it to the subproblem */ 937 SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, SCIProwGetName(rows[i]), nnonz, consvars, vals, lhs, rhs, 938 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) ); 939 SCIP_CALL( SCIPaddCons(subscip, cons) ); 940 SCIP_CALL( SCIPreleaseCons(subscip, &cons) ); 941 942 /* free temporary memory */ 943 SCIPfreeBufferArray(scip, &consvars); 944 } 945 946 return SCIP_OKAY; 947 } 948 949 950 /** get a sub-SCIP copy of the transformed problem */ 951 SCIP_RETCODE SCIPcopyLargeNeighborhoodSearch( 952 SCIP* sourcescip, /**< source SCIP data structure */ 953 SCIP* subscip, /**< sub-SCIP used by the heuristic */ 954 SCIP_HASHMAP* varmap, /**< a hashmap to store the mapping of source variables to the corresponding 955 * target variables */ 956 const char* suffix, /**< suffix for the problem name */ 957 SCIP_VAR** fixedvars, /**< source variables whose copies should be fixed in the target SCIP environment, or NULL */ 958 SCIP_Real* fixedvals, /**< array of fixing values for target SCIP variables, or NULL */ 959 int nfixedvars, /**< number of source variables whose copies should be fixed in the target SCIP environment, or NULL */ 960 SCIP_Bool uselprows, /**< should the linear relaxation of the problem defined by LP rows be copied? */ 961 SCIP_Bool copycuts, /**< should cuts be copied (only if uselprows == FALSE) */ 962 SCIP_Bool* success, /**< was the copying successful? */ 963 SCIP_Bool* valid /**< pointer to store whether the copying was valid, or NULL */ 964 ) 965 { 966 assert(sourcescip != NULL); 967 assert(suffix != NULL); 968 assert(subscip != NULL); 969 assert(varmap != NULL); 970 assert(success != NULL); 971 972 if( uselprows ) 973 { 974 char probname[SCIP_MAXSTRLEN]; 975 976 /* copy all plugins */ 977 SCIP_CALL( SCIPincludeDefaultPlugins(subscip) ); 978 979 /* set name to the original problem name and possibly add a suffix */ 980 (void) SCIPsnprintf(probname, SCIP_MAXSTRLEN, "%s_%s", SCIPgetProbName(sourcescip), suffix); 981 982 /* create the subproblem */ 983 SCIP_CALL( SCIPcreateProb(subscip, probname, NULL, NULL, NULL, NULL, NULL, NULL, NULL) ); 984 985 /* copy all variables */ 986 SCIP_CALL( SCIPcopyVars(sourcescip, subscip, varmap, NULL, fixedvars, fixedvals, nfixedvars, TRUE) ); 987 988 /* copy parameter settings */ 989 SCIP_CALL( SCIPcopyParamSettings(sourcescip, subscip) ); 990 991 /* create linear constraints from LP rows of the source problem */ 992 SCIP_CALL( createRows(sourcescip, subscip, varmap) ); 993 } 994 else 995 { 996 SCIP_CALL( SCIPcopyConsCompression(sourcescip, subscip, varmap, NULL, suffix, fixedvars, fixedvals, nfixedvars, 997 TRUE, FALSE, FALSE, TRUE, valid) ); 998 999 if( copycuts ) 1000 { 1001 /* copies all active cuts from cutpool of sourcescip to linear constraints in targetscip */ 1002 SCIP_CALL( SCIPcopyCuts(sourcescip, subscip, varmap, NULL, TRUE, NULL) ); 1003 } 1004 } 1005 1006 *success = TRUE; 1007 1008 return SCIP_OKAY; 1009 } 1010 1011 /** adds a trust region neighborhood constraint to the @p targetscip 1012 * 1013 * a trust region constraint measures the deviation from the current incumbent solution \f$x^*\f$ by an auxiliary 1014 * continuous variable \f$v \geq 0\f$: 1015 * \f[ 1016 * \sum\limits_{j\in B} |x_j^* - x_j| = v 1017 * \f] 1018 * Only binary variables are taken into account. The deviation is penalized in the objective function using 1019 * a positive \p violpenalty. 1020 * 1021 * @note: the trust region constraint creates an auxiliary variable to penalize the deviation from 1022 * the current incumbent solution. This variable can afterwards be accessed using SCIPfindVar() by its name 1023 * 'trustregion_violationvar' 1024 */ 1025 SCIP_RETCODE SCIPaddTrustregionNeighborhoodConstraint( 1026 SCIP* sourcescip, /**< the data structure for the main SCIP instance */ 1027 SCIP* targetscip, /**< SCIP data structure of the subproblem */ 1028 SCIP_VAR** subvars, /**< variables of the subproblem, NULL entries are ignored */ 1029 SCIP_Real violpenalty /**< the penalty for violating the trust region */ 1030 ) 1031 { 1032 SCIP_VAR* violvar; 1033 SCIP_CONS* trustregioncons; 1034 SCIP_VAR** consvars; 1035 SCIP_VAR** vars; 1036 SCIP_SOL* bestsol; 1037 1038 int nvars; 1039 int nbinvars; 1040 int nconsvars; 1041 int i; 1042 SCIP_Real rhs; 1043 SCIP_Real* consvals; 1044 char name[SCIP_MAXSTRLEN]; 1045 1046 /* get the data of the variables and the best solution */ 1047 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, &nvars, &nbinvars, NULL, NULL, NULL) ); 1048 bestsol = SCIPgetBestSol(sourcescip); 1049 assert(bestsol != NULL); 1050 /* otherwise, this neighborhood would not be active in the first place */ 1051 assert(nbinvars > 0); 1052 1053 /* memory allocation */ 1054 SCIP_CALL( SCIPallocBufferArray(sourcescip, &consvars, nbinvars + 1) ); 1055 SCIP_CALL( SCIPallocBufferArray(sourcescip, &consvals, nbinvars + 1) ); 1056 nconsvars = 0; 1057 1058 /* set initial left and right hand sides of trust region constraint */ 1059 rhs = 0.0; 1060 1061 /* create the distance (to incumbent) function of the binary variables */ 1062 for( i = 0; i < nbinvars; i++ ) 1063 { 1064 SCIP_Real solval; 1065 1066 if( subvars[i] == NULL ) 1067 continue; 1068 1069 solval = SCIPgetSolVal(sourcescip, bestsol, vars[i]); 1070 assert( SCIPisFeasIntegral(sourcescip,solval) ); 1071 1072 /* is variable i part of the binary support of bestsol? */ 1073 if( SCIPisFeasEQ(sourcescip, solval, 1.0) ) 1074 { 1075 consvals[nconsvars] = -1.0; 1076 rhs -= 1.0; 1077 } 1078 else 1079 consvals[nconsvars] = 1.0; 1080 consvars[nconsvars] = subvars[i]; 1081 assert(SCIPvarGetType(consvars[nconsvars]) == SCIP_VARTYPE_BINARY); 1082 ++nconsvars; 1083 } 1084 1085 /* adding the violation variable */ 1086 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_trustregionviolvar", SCIPgetProbName(sourcescip)); 1087 SCIP_CALL( SCIPcreateVarBasic(targetscip, &violvar, name, 0.0, SCIPinfinity(targetscip), violpenalty, SCIP_VARTYPE_CONTINUOUS) ); 1088 SCIP_CALL( SCIPaddVar(targetscip, violvar) ); 1089 consvars[nconsvars] = violvar; 1090 consvals[nconsvars] = -1.0; 1091 ++nconsvars; 1092 1093 /* creates trustregion constraint and adds it to subscip */ 1094 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_trustregioncons", SCIPgetProbName(sourcescip)); 1095 1096 SCIP_CALL( SCIPcreateConsLinear(targetscip, &trustregioncons, name, nconsvars, consvars, consvals, 1097 rhs, rhs, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) ); 1098 SCIP_CALL( SCIPaddCons(targetscip, trustregioncons) ); 1099 SCIP_CALL( SCIPreleaseCons(targetscip, &trustregioncons) ); 1100 1101 /* releasing the violation variable */ 1102 SCIP_CALL( SCIPreleaseVar(targetscip, &violvar) ); 1103 1104 /* free local memory */ 1105 SCIPfreeBufferArray(sourcescip, &consvals); 1106 SCIPfreeBufferArray(sourcescip, &consvars); 1107 1108 return SCIP_OKAY; 1109 } 1110