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_dps.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief dynamic partition search 28 * @author Katrin Halbig 29 * 30 * The dynamic partition search (DPS) is a construction heuristic which additionally needs a 31 * user decomposition with linking constraints only. 32 * 33 * This heuristic splits the problem into several sub-SCIPs according to the given decomposition. Thereby the linking constraints 34 * with their right-hand and left-hand sides are also split. DPS searches for a partition of the sides on the blocks 35 * so that a feasible solution is obtained. 36 * For each block the parts of the original linking constraints are extended by slack variables. Moreover, the objective function 37 * is replaced by the sum of these additional variables weighted by penalty parameters lambda. If all blocks have an optimal solution 38 * of zero, the algorithm terminates with a feasible solution for the main problem. Otherwise, the partition and the penalty parameters 39 * are updated, and the sub-SCIPs are solved again. 40 * 41 * A detailed description can be found in 42 * K. Halbig, A. Göß and D. Weninger (2023). Exploiting user-supplied Decompositions inside Heuristics. https://optimization-online.org/?p=23386 43 */ 44 45 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 46 47 #include "scip/heur_dps.h" 48 #include "scip/pub_dcmp.h" 49 #include "scip/pub_heur.h" 50 #include "scip/pub_misc.h" 51 #include "scip/scipdefplugins.h" 52 #include "scip/scip_cons.h" 53 #include "scip/scip_dcmp.h" 54 #include "scip/scip_general.h" 55 #include "scip/scip_heur.h" 56 #include "scip/scip_mem.h" 57 #include "scip/scip_message.h" 58 #include "scip/scip_param.h" 59 #include "scip/scip_prob.h" 60 61 62 #define HEUR_NAME "dps" 63 #define HEUR_DESC "primal heuristic for decomposable MIPs" 64 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_LNS 65 #define HEUR_PRIORITY 75000 66 #define HEUR_FREQ -1 67 #define HEUR_FREQOFS 0 68 #define HEUR_MAXDEPTH -1 69 #define HEUR_TIMING SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_AFTERNODE 70 #define HEUR_USESSUBSCIP TRUE /**< does the heuristic use a secondary SCIP instance? */ 71 72 #define DEFAULT_MAXIT 50 /**< maximum number of iterations */ 73 #define DEFAULT_PENALTY 100.0 /**< multiplier for absolute increase of penalty parameters */ 74 75 /* event handler properties */ 76 #define EVENTHDLR_NAME "Dps" 77 #define EVENTHDLR_DESC "event handler for " HEUR_NAME " heuristic" 78 79 /* 80 * Data structures 81 */ 82 83 /** primal heuristic data */ 84 struct SCIP_HeurData 85 { 86 SCIP_CONS** linkingconss; /**< linking constraints */ 87 int nlinking; /**< number of linking constraints */ 88 int nblocks; /**< number of blocks */ 89 int maxit; /**< maximal number of iterations */ 90 int timing; /**< should the heuristic run before or after the processing of the node? 91 (0: before, 1: after, 2: both)*/ 92 SCIP_Real maxlinkscore; /**< maximal linking score of used decomposition (equivalent to percentage of linking constraints) */ 93 SCIP_Real penalty; /**< multiplier for absolute increase of penalty parameters */ 94 SCIP_Bool reoptimize; /**< should the problem get reoptimized with the original objective function? */ 95 SCIP_Bool reuse; /**< should solutions get reused in subproblems? */ 96 SCIP_Bool reoptlimits; /**< should strict limits for reoptimization be set? */ 97 }; 98 99 /** data related to one block */ 100 struct Blockproblem 101 { 102 SCIP* blockscip; /**< SCIP data structure */ 103 SCIP_VAR** slackvars; /**< slack variables */ 104 SCIP_CONS** linkingconss; /**< linking constraints */ 105 int* linkingindices; /**< indices of linking constraints in original problem */ 106 int nlinking; /**< number of linking constraints */ 107 int nblockvars; /**< number of block variables */ 108 int nslackvars; /**< number of slack variables */ 109 SCIP_Real* origobj; /**< original objective coefficients */ 110 }; 111 typedef struct Blockproblem BLOCKPROBLEM; 112 113 /** data related to one linking constraint */ 114 struct Linking 115 { 116 SCIP_CONS* linkingcons; /**< corresponding linking constraint of original problem */ 117 SCIP_CONS** blockconss; /**< linking constraints of the blocks */ 118 SCIP_VAR** slacks; /**< slackvars of the blocks */ 119 SCIP_Real* minactivity; /**< minimal activity of constraint for each block */ 120 SCIP_Real* maxactivity; /**< maximal activity of constraint for each block */ 121 SCIP_Real* currentrhs; /**< current partition of rhs */ 122 SCIP_Real* currentlhs; /**< current partition of lhs */ 123 int* blocknumbers; /**< number of the blocks */ 124 int nblocks; /**< number of blocks in which this linking constraint participates; dimension of arrays */ 125 int nslacks; /**< number of slack variables */ 126 int nslacksperblock; /**< 2, if ranged constraint; 1, if only rhs or lhs */ 127 int lastviolations; /**< number of iterations in which the constraint was violated in succession */ 128 SCIP_Bool hasrhs; /**< has linking constraint finite right-hand side? */ 129 SCIP_Bool haslhs; /**< has linking constraint finite left-hand side? */ 130 }; 131 typedef struct Linking LINKING; 132 133 /* 134 * Local methods 135 */ 136 137 /** assigns linking variables to last block 138 * 139 * The labels are copied to newdecomp and the linking variables are assigned to the last block (i.e., highest block label). 140 * Constraint labels and statistics are recomputed. 141 */ 142 static 143 SCIP_RETCODE assignLinking( 144 SCIP* scip, /**< SCIP data structure */ 145 SCIP_DECOMP* newdecomp, /**< decomposition with assigned linking variables */ 146 SCIP_VAR** vars, /**< sorted array of variables */ 147 SCIP_CONS** conss, /**< sorted array of constraints */ 148 int* varlabels, /**< sorted array of variable labels */ 149 int* conslabels, /**< sorted array of constraint labels */ 150 int nvars, /**< number of variables */ 151 int nconss, /**< number of constraints */ 152 int nlinkvars /**< number of linking variables */ 153 ) 154 { 155 int newlabel; 156 int v; 157 158 assert(scip != NULL); 159 assert(newdecomp != NULL); 160 assert(vars != NULL); 161 assert(conss != NULL); 162 assert(varlabels != NULL); 163 assert(conslabels != NULL); 164 165 /* copy the labels */ 166 SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, vars, varlabels, nvars) ); 167 SCIP_CALL( SCIPdecompSetConsLabels(newdecomp, conss, conslabels, nconss) ); 168 169 /* assign linking variables */ 170 newlabel = varlabels[nvars - 1]; /* take always label of last block */ 171 assert(newlabel >= 0); 172 for( v = 0; v < nlinkvars; v++ ) 173 { 174 SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, &vars[v], &newlabel, 1) ); 175 } 176 SCIPdebugMsg(scip, "assigned %d linking variables\n", nlinkvars); 177 178 /* recompute constraint labels and statistics */ 179 SCIP_CALL( SCIPcomputeDecompConsLabels(scip, newdecomp, conss, nconss) ); 180 SCIP_CALL( SCIPcomputeDecompStats(scip, newdecomp, TRUE) ); 181 nlinkvars = SCIPdecompGetNBorderVars(newdecomp); 182 183 /* get new labels and sort */ 184 SCIPdecompGetConsLabels(newdecomp, conss, conslabels, nconss); 185 SCIPdecompGetVarsLabels(newdecomp, vars, varlabels, nvars); 186 SCIPsortIntPtr(conslabels, (void**)conss, nconss); 187 SCIPsortIntPtr(varlabels, (void**)vars, nvars); 188 189 /* After assigning the linking variables, blocks can have zero constraints. 190 * So the remaining variables are labeled as linking in SCIPcomputeDecompStats(). 191 * We assign this variables to the same label as above. 192 */ 193 if( nlinkvars >= 1 ) 194 { 195 assert(varlabels[0] == SCIP_DECOMP_LINKVAR); 196 SCIPdebugMsg(scip, "assign again %d linking variables\n", nlinkvars); 197 198 for( v = 0; v < nlinkvars; v++ ) 199 { 200 SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, &vars[v], &newlabel, 1) ); 201 } 202 SCIP_CALL( SCIPcomputeDecompConsLabels(scip, newdecomp, conss, nconss) ); 203 SCIP_CALL( SCIPcomputeDecompStats(scip, newdecomp, TRUE) ); 204 205 SCIPdecompGetConsLabels(newdecomp, conss, conslabels, nconss); 206 SCIPdecompGetVarsLabels(newdecomp, vars, varlabels, nvars); 207 SCIPsortIntPtr(conslabels, (void**)conss, nconss); 208 SCIPsortIntPtr(varlabels, (void**)vars, nvars); 209 } 210 assert(varlabels[0] != SCIP_DECOMP_LINKVAR); 211 212 return SCIP_OKAY; 213 } 214 215 /** creates a sub-SCIP and sets parameters */ 216 static 217 SCIP_RETCODE createSubscip( 218 SCIP* scip, /**< main SCIP data structure */ 219 SCIP** subscip /**< pointer to store created sub-SCIP */ 220 ) 221 { 222 SCIP_Real infvalue; 223 224 assert(scip != NULL); 225 assert(subscip != NULL); 226 227 /* create a new SCIP instance */ 228 SCIP_CALL( SCIPcreate(subscip) ); 229 SCIP_CALL( SCIPincludeDefaultPlugins(*subscip) ); 230 231 /* copy value for infinity */ 232 SCIP_CALL( SCIPgetRealParam(scip, "numerics/infinity", &infvalue) ); 233 SCIP_CALL( SCIPsetRealParam(*subscip, "numerics/infinity", infvalue) ); 234 235 SCIP_CALL( SCIPcopyLimits(scip, *subscip) ); 236 237 /* avoid recursive calls */ 238 SCIP_CALL( SCIPsetSubscipsOff(*subscip, TRUE) ); 239 240 /* disable cutting plane separation */ 241 SCIP_CALL( SCIPsetSeparating(*subscip, SCIP_PARAMSETTING_OFF, TRUE) ); 242 243 /* disable expensive presolving */ 244 SCIP_CALL( SCIPsetPresolving(*subscip, SCIP_PARAMSETTING_FAST, TRUE) ); 245 246 /* disable expensive techniques */ 247 SCIP_CALL( SCIPsetIntParam(*subscip, "misc/usesymmetry", 0) ); 248 249 /* do not abort subproblem on CTRL-C */ 250 SCIP_CALL( SCIPsetBoolParam(*subscip, "misc/catchctrlc", FALSE) ); 251 252 #ifdef SCIP_DEBUG 253 /* for debugging, enable full output */ 254 SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 5) ); 255 SCIP_CALL( SCIPsetIntParam(*subscip, "display/freq", 100000000) ); 256 #else 257 /* disable statistic timing inside sub SCIP and output to console */ 258 SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 0) ); 259 SCIP_CALL( SCIPsetBoolParam(*subscip, "timing/statistictiming", FALSE) ); 260 #endif 261 262 /* speed up sub-SCIP by not checking dual LP feasibility */ 263 SCIP_CALL( SCIPsetBoolParam(*subscip, "lp/checkdualfeas", FALSE) ); 264 265 return SCIP_OKAY; 266 } 267 268 /** copies the given variables and constraints to the given sub-SCIP */ 269 static 270 SCIP_RETCODE copyToSubscip( 271 SCIP* scip, /**< source SCIP */ 272 SCIP* subscip, /**< target SCIP */ 273 const char* name, /**< name for copied problem */ 274 SCIP_VAR** vars, /**< array of variables to copy */ 275 SCIP_CONS** conss, /**< array of constraints to copy */ 276 SCIP_HASHMAP* varsmap, /**< hashmap for copied variables */ 277 SCIP_HASHMAP* conssmap, /**< hashmap for copied constraints */ 278 int nvars, /**< number of variables to copy */ 279 int nconss, /**< number of constraints to copy */ 280 SCIP_Bool* success /**< was copying successful? */ 281 ) 282 { 283 SCIP_CONS* newcons; 284 SCIP_VAR* newvar; 285 int i; 286 287 assert(scip != NULL); 288 assert(subscip != NULL); 289 assert(vars != NULL); 290 assert(conss != NULL); 291 assert(varsmap != NULL); 292 assert(conssmap != NULL); 293 assert(success != NULL); 294 295 SCIPdebugMsg(scip, "copyToSubscip\n"); 296 297 /* create problem in sub-SCIP */ 298 SCIP_CALL( SCIPcreateProb(subscip, name, NULL, NULL, NULL, NULL, NULL, NULL, NULL) ); 299 300 /* copy variables */ 301 for( i = 0; i < nvars; ++i ) 302 { 303 SCIP_CALL( SCIPgetVarCopy(scip, subscip, vars[i], &newvar, varsmap, conssmap, FALSE, success) ); 304 305 /* abort if variable was not successfully copied */ 306 if( !(*success) ) 307 { 308 SCIPwarningMessage(scip, "Abort heuristic dps since not all variables were successfully copied.\n"); 309 SCIPABORT(); 310 return SCIP_OKAY; 311 } 312 } 313 assert(nvars == SCIPgetNOrigVars(subscip)); 314 315 /* copy constraints */ 316 for( i = 0; i < nconss; ++i ) 317 { 318 assert(conss[i] != NULL); 319 assert(!SCIPconsIsModifiable(conss[i])); 320 assert(SCIPconsIsActive(conss[i])); 321 assert(!SCIPconsIsDeleted(conss[i])); 322 323 /* copy the constraint */ 324 SCIP_CALL( SCIPgetConsCopy(scip, subscip, conss[i], &newcons, SCIPconsGetHdlr(conss[i]), varsmap, conssmap, NULL, 325 SCIPconsIsInitial(conss[i]), SCIPconsIsSeparated(conss[i]), SCIPconsIsEnforced(conss[i]), 326 SCIPconsIsChecked(conss[i]), SCIPconsIsPropagated(conss[i]), FALSE, FALSE, 327 SCIPconsIsDynamic(conss[i]), SCIPconsIsRemovable(conss[i]), FALSE, FALSE, success) ); 328 329 /* abort if constraint was not successfully copied */ 330 if( !(*success) ) 331 return SCIP_OKAY; 332 333 SCIP_CALL( SCIPaddCons(subscip, newcons) ); 334 SCIP_CALL( SCIPreleaseCons(subscip, &newcons) ); 335 } 336 337 /* block constraint contains variables which are not part of this block 338 * 339 * todo: maybe they are part of the block, but it is not recognized, because they are, for example, negated or aggregated. 340 */ 341 if( nvars != SCIPgetNOrigVars(subscip) ) 342 *success = FALSE; 343 344 return SCIP_OKAY; 345 } 346 347 /** creates the subscip for a given block */ 348 static 349 SCIP_RETCODE createBlockproblem( 350 SCIP* scip, /**< SCIP data structure */ 351 BLOCKPROBLEM* blockproblem, /**< blockproblem that should be created */ 352 LINKING** linkings, /**< linkings that will be (partially) initialized */ 353 SCIP_CONS** conss, /**< sorted array of constraints of this block */ 354 SCIP_VAR** vars, /**< sorted array of variables of this block */ 355 int nconss, /**< number of constraints of this block */ 356 int nvars, /**< number of variables of this block */ 357 SCIP_CONS** linkingconss, /**< linking constraints in the original problem */ 358 int nlinking, /**< number of linking constraints in the original problem */ 359 int blocknumber, /**< number of block that should be created */ 360 SCIP_Bool* success /**< pointer to store whether creation was successful */ 361 ) 362 { 363 char name[SCIP_MAXSTRLEN]; 364 SCIP_HASHMAP* varsmap; 365 SCIP_HASHMAP* conssmap; 366 SCIP_VAR** consvars; /* all vars in original linking cons */ 367 SCIP_Real* consvals; 368 int nconsvars; 369 SCIP_VAR** blockvars; /* vars of current linking cons of current block */ 370 SCIP_Real* blockvals; 371 int nblockvars; 372 SCIP_VAR** subvars; /* all vars of subscip */ 373 int maxnconsvars; /* current size of arrays */ 374 int c; 375 int v; 376 377 assert(scip != NULL); 378 assert(blockproblem != NULL); 379 assert(conss != NULL); 380 assert(vars != NULL); 381 assert(blockproblem->blockscip != NULL); 382 383 maxnconsvars = 20; /* start size; increase size if necessary */ 384 385 SCIPdebugMsg(scip, "Create blockproblem %d\n", blocknumber); 386 387 /* create the variable/constraint mapping hash map */ 388 SCIP_CALL( SCIPhashmapCreate(&varsmap, SCIPblkmem(scip), nvars) ); 389 SCIP_CALL( SCIPhashmapCreate(&conssmap, SCIPblkmem(scip), nconss) ); 390 391 /* get name of the original problem and add "comp_nr" */ 392 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_comp_%d", SCIPgetProbName(scip), blocknumber); 393 394 SCIP_CALL( copyToSubscip(scip, blockproblem->blockscip, name, vars, conss, varsmap, conssmap, 395 nvars, nconss, success) ); 396 if( !(*success) ) 397 { 398 SCIPdebugMsg(scip, "Copy to subscip failed\n"); 399 SCIPhashmapFree(&conssmap); 400 SCIPhashmapFree(&varsmap); 401 402 return SCIP_OKAY; 403 } 404 405 /* save number of variables that have a corresponding variable in original problem*/ 406 blockproblem->nblockvars = SCIPgetNVars(blockproblem->blockscip); 407 408 /* save original objective and set objective to zero */ 409 subvars = SCIPgetVars(blockproblem->blockscip); 410 for( v = 0; v < nvars; v++ ) 411 { 412 blockproblem->origobj[v] = SCIPvarGetObj(subvars[v]); 413 SCIP_CALL( SCIPchgVarObj(blockproblem->blockscip, subvars[v], 0.0) ); 414 } 415 416 /* allocate memory */ 417 SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &blockvars, nvars + 2) ); /* two entries for the slack variables */ 418 SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &blockvals, nvars + 2) ); 419 SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &consvars, maxnconsvars) ); 420 SCIP_CALL( SCIPallocBufferArray(blockproblem->blockscip, &consvals, maxnconsvars) ); 421 422 /* find and add parts of linking constraints */ 423 SCIPdebugMsg(scip, "add parts of linking constraints\n"); 424 for( c = 0; c < nlinking; c++ ) 425 { 426 const char* conshdlrname; 427 char consname[SCIP_MAXSTRLEN]; 428 SCIP_CONS* newcons; 429 SCIP_Real rhs; 430 SCIP_Real lhs; 431 SCIP_Real minact; 432 SCIP_Real maxact; 433 SCIP_Bool mininfinite; 434 SCIP_Bool maxinfinite; 435 436 assert(linkingconss[c] != NULL); 437 438 newcons = NULL; 439 440 #ifdef SCIP_MORE_DEBUG 441 SCIPdebugMsg(scip, "consider constraint %s\n", SCIPconsGetName(linkingconss[c])); 442 SCIPdebugPrintCons(scip, linkingconss[c], NULL); 443 #endif 444 445 nblockvars = 0; 446 447 /* every constraint with linear representation is allowed */ 448 conshdlrname = SCIPconshdlrGetName(SCIPconsGetHdlr(linkingconss[c])); 449 if( !( (strcmp(conshdlrname, "linear") == 0) || (strcmp(conshdlrname, "setppc") == 0) 450 || (strcmp(conshdlrname, "logicor") == 0) || (strcmp(conshdlrname, "knapsack") == 0) 451 || (strcmp(conshdlrname, "varbound") == 0) ) ) 452 { 453 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, "Heuristic %s cannot handle linking constraints of type %s\n", HEUR_NAME, conshdlrname); 454 /* TODO which other types can we handle/transform in a linear constraint? */ 455 456 *success = FALSE; 457 break; /* releases memory and breaks heuristic */ 458 } 459 460 SCIP_CALL( SCIPgetConsNVars(scip, linkingconss[c], &nconsvars, success) ); 461 462 /* reallocate memory if we have more variables than maxnconsvars */ 463 if( nconsvars > maxnconsvars ) 464 { 465 int newsize; 466 467 /* calculate new size */ 468 newsize = SCIPcalcMemGrowSize(scip, MAX(2 * maxnconsvars, nconsvars)); /* at least double size */ 469 470 SCIP_CALL( SCIPreallocBufferArray(blockproblem->blockscip, &consvars, newsize) ); 471 SCIP_CALL( SCIPreallocBufferArray(blockproblem->blockscip, &consvals, newsize) ); 472 maxnconsvars = newsize; 473 } 474 475 SCIP_CALL( SCIPgetConsVars(scip, linkingconss[c], consvars, nconsvars, success) ); 476 SCIP_CALL( SCIPgetConsVals(scip, linkingconss[c], consvals, nconsvars, success) ); 477 478 if( !(*success) ) 479 { 480 SCIPdebugMsg(scip, "Create blockproblem failed\n"); 481 break; /* releases memory and breaks heuristic */ 482 } 483 484 /* check if constraint contains variables of this block */ 485 for( v = 0; v < nconsvars; v++ ) 486 { 487 if( SCIPhashmapExists(varsmap, (void*)consvars[v]) ) 488 { 489 blockvars[nblockvars] = (SCIP_VAR*) SCIPhashmapGetImage(varsmap, (void*)consvars[v]); 490 blockvals[nblockvars] = consvals[v]; 491 ++nblockvars; 492 } 493 /* handle negated variables*/ 494 else if( SCIPvarGetStatus(consvars[v]) == SCIP_VARSTATUS_NEGATED) 495 { 496 if( SCIPhashmapExists(varsmap, (void*)SCIPvarGetNegationVar(consvars[v])) ) /* negation exists in this block */ 497 { 498 /* save negated variable */ 499 SCIP_VAR* origblockvar = (SCIP_VAR*) SCIPhashmapGetImage(varsmap, (void*)SCIPvarGetNegationVar(consvars[v])); 500 SCIP_VAR* negblockvar = NULL; 501 SCIP_CALL( SCIPgetNegatedVar(blockproblem->blockscip, origblockvar, &negblockvar) ); 502 blockvars[nblockvars] = negblockvar; 503 blockvals[nblockvars] = consvals[v]; 504 ++nblockvars; 505 } 506 } 507 } 508 509 /* continue with next linking constraint if it has no part in current block */ 510 if( nblockvars == 0 ) 511 continue; 512 513 /* get rhs and/or lhs */ 514 rhs = SCIPconsGetRhs(scip, linkingconss[c], success); 515 if( !(*success) ) 516 { 517 SCIPdebugMsg(scip, "Create blockproblem failed\n"); 518 return SCIP_OKAY; 519 } 520 lhs = SCIPconsGetLhs(scip, linkingconss[c], success); 521 if( !(*success) ) 522 { 523 SCIPdebugMsg(scip, "Create blockproblem failed\n"); 524 return SCIP_OKAY; 525 } 526 assert(!SCIPisInfinity(scip, rhs) || !SCIPisInfinity(scip, -lhs)); /* at least one side bounded */ 527 assert(SCIPisLE(scip, lhs, rhs)); 528 529 if( !SCIPisInfinity(scip, rhs) ) 530 linkings[c]->hasrhs = TRUE; 531 if( !SCIPisInfinity(scip, -lhs) ) 532 linkings[c]->haslhs = TRUE; 533 if( !SCIPisInfinity(scip, rhs) && !SCIPisInfinity(scip, -lhs)) 534 linkings[c]->nslacksperblock = 2; 535 else 536 linkings[c]->nslacksperblock = 1; 537 538 /* add slack variable for rhs */ 539 if( linkings[c]->hasrhs ) 540 { 541 /* slack variable z_r >= 0 */ 542 char varname[SCIP_MAXSTRLEN]; 543 (void)SCIPsnprintf(varname, SCIP_MAXSTRLEN, "z_r_%s", SCIPconsGetName(linkingconss[c])); 544 SCIP_CALL( SCIPcreateVarBasic(blockproblem->blockscip, &blockvars[nblockvars], varname, 545 0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) ); 546 blockvals[nblockvars] = -1.0; 547 SCIP_CALL( SCIPaddVar(blockproblem->blockscip, blockvars[nblockvars]) ); 548 #ifdef SCIP_MORE_DEBUG 549 SCIPdebugMsg(scip, "Add variable %s\n", SCIPvarGetName(blockvars[nblockvars])); 550 #endif 551 linkings[c]->slacks[linkings[c]->nslacks] = blockvars[nblockvars]; 552 blockproblem->slackvars[blockproblem->nslackvars] = blockvars[nblockvars]; 553 ++blockproblem->nslackvars; 554 ++linkings[c]->nslacks; 555 ++nblockvars; 556 } 557 558 /* add slack variable for lhs */ 559 if( linkings[c]->haslhs ) 560 { 561 /* slack variable z_l >= 0 */ 562 char varname[SCIP_MAXSTRLEN]; 563 (void)SCIPsnprintf(varname, SCIP_MAXSTRLEN, "z_l_%s", SCIPconsGetName(linkingconss[c])); 564 SCIP_CALL( SCIPcreateVarBasic(blockproblem->blockscip, &blockvars[nblockvars], varname, 565 0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) ); 566 blockvals[nblockvars] = 1.0; 567 SCIP_CALL( SCIPaddVar(blockproblem->blockscip, blockvars[nblockvars]) ); 568 #ifdef SCIP_MORE_DEBUG 569 SCIPdebugMsg(scip, "Add variable %s\n", SCIPvarGetName(blockvars[nblockvars])); 570 #endif 571 linkings[c]->slacks[linkings[c]->nslacks] = blockvars[nblockvars]; 572 blockproblem->slackvars[blockproblem->nslackvars] = blockvars[nblockvars]; 573 ++blockproblem->nslackvars; 574 ++linkings[c]->nslacks; 575 ++nblockvars; 576 } 577 578 /* add linking constraint with slack variable */ 579 (void)SCIPsnprintf(consname, SCIP_MAXSTRLEN, "%s", SCIPconsGetName(linkingconss[c])); 580 SCIP_CALL( SCIPcreateConsBasicLinear(blockproblem->blockscip, &newcons, consname, nblockvars, blockvars, blockvals, lhs, rhs) ); 581 SCIP_CALL( SCIPaddCons(blockproblem->blockscip, newcons) ); 582 #ifdef SCIP_MORE_DEBUG 583 SCIPdebugMsg(blockproblem->blockscip, "add constraint %s\n", SCIPconsGetName(newcons)); 584 SCIPdebugPrintCons(blockproblem->blockscip, newcons, NULL); 585 #endif 586 587 blockproblem->linkingconss[blockproblem->nlinking] = newcons; 588 linkings[c]->blockconss[linkings[c]->nblocks] = newcons; 589 linkings[c]->blocknumbers[linkings[c]->nblocks] = blocknumber; 590 blockproblem->linkingindices[blockproblem->nlinking] = c; 591 592 /* calculate minimal und maximal activity (exclude slack variables) */ 593 minact = 0; 594 maxact = 0; 595 mininfinite = FALSE; 596 maxinfinite = FALSE; 597 for( v = 0; v < nblockvars - linkings[c]->nslacksperblock && (!mininfinite || !maxinfinite); v++ ) 598 { 599 SCIP_Real lb; 600 SCIP_Real ub; 601 lb = SCIPvarGetLbGlobal(blockvars[v]); 602 ub = SCIPvarGetUbGlobal(blockvars[v]); 603 604 if( blockvals[v] >= 0.0 ) 605 { 606 mininfinite = (mininfinite || SCIPisInfinity(scip, -lb)); 607 maxinfinite = (maxinfinite || SCIPisInfinity(scip, ub)); 608 if( !mininfinite ) 609 minact += blockvals[v] * lb; 610 if( !maxinfinite ) 611 maxact += blockvals[v] * ub; 612 } 613 else 614 { 615 mininfinite = (mininfinite || SCIPisInfinity(scip, ub)); 616 maxinfinite = (maxinfinite || SCIPisInfinity(scip, -lb)); 617 if( !mininfinite ) 618 minact += blockvals[v] * ub; 619 if( !maxinfinite ) 620 maxact += blockvals[v] * lb; 621 } 622 } 623 624 if( mininfinite ) 625 linkings[c]->minactivity[linkings[c]->nblocks] = -SCIPinfinity(scip); 626 else 627 linkings[c]->minactivity[linkings[c]->nblocks] = minact; 628 if( maxinfinite ) 629 linkings[c]->maxactivity[linkings[c]->nblocks] = SCIPinfinity(scip); 630 else 631 linkings[c]->maxactivity[linkings[c]->nblocks] = maxact; 632 assert(SCIPisLE(scip, linkings[c]->minactivity[linkings[c]->nblocks], linkings[c]->maxactivity[linkings[c]->nblocks])); 633 634 linkings[c]->nblocks++; 635 blockproblem->nlinking++; 636 637 for( v = 1; v <= linkings[c]->nslacksperblock; v++ ) 638 { 639 SCIP_CALL( SCIPreleaseVar(blockproblem->blockscip, &blockvars[nblockvars - v]) ); 640 } 641 642 SCIP_CALL( SCIPreleaseCons(blockproblem->blockscip, &newcons) ); 643 } 644 assert(blockproblem->nlinking <= nlinking); 645 646 /* free memory */ 647 SCIPfreeBufferArray(blockproblem->blockscip, &consvals); 648 SCIPfreeBufferArray(blockproblem->blockscip, &consvars); 649 SCIPfreeBufferArray(blockproblem->blockscip, &blockvals); 650 SCIPfreeBufferArray(blockproblem->blockscip, &blockvars); 651 652 SCIPhashmapFree(&conssmap); 653 SCIPhashmapFree(&varsmap); 654 655 return SCIP_OKAY; 656 } 657 658 /** creates data structures and splits problem into blocks */ 659 static 660 SCIP_RETCODE createAndSplitProblem( 661 SCIP* scip, /**< SCIP data structure */ 662 SCIP_HEURDATA* heurdata, /**< heuristic data */ 663 SCIP_DECOMP* decomp, /**< decomposition data structure */ 664 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 665 LINKING** linkings, /**< array of linking data structures */ 666 SCIP_VAR** vars, /**< sorted array of variables */ 667 SCIP_CONS** conss, /**< sorted array of constraints */ 668 SCIP_Bool* success /**< pointer to store whether splitting was successful */ 669 ) 670 { 671 int* nconssblock; 672 int* nvarsblock; 673 int conssoffset; 674 int varsoffset; 675 int i; /* blocknumber */ 676 677 assert(scip != NULL); 678 assert(heurdata != NULL); 679 assert(vars != NULL); 680 assert(conss != NULL); 681 682 SCIP_CALL( SCIPallocBufferArray(scip, &nvarsblock, heurdata->nblocks + 1) ); 683 SCIP_CALL( SCIPallocBufferArray(scip, &nconssblock, heurdata->nblocks + 1) ); 684 SCIP_CALL( SCIPdecompGetVarsSize(decomp, nvarsblock, heurdata->nblocks + 1) ); 685 SCIP_CALL( SCIPdecompGetConssSize(decomp, nconssblock, heurdata->nblocks + 1) ); 686 assert(0 == nvarsblock[0]); 687 688 varsoffset = 0; 689 conssoffset = 0; 690 691 for( i = 0; i < heurdata->nblocks; i++) 692 { 693 conssoffset += nconssblock[i]; 694 varsoffset += nvarsblock[i]; 695 696 SCIP_CALL( createBlockproblem(scip, blockproblem[i], linkings, &conss[conssoffset], &vars[varsoffset], nconssblock[i+1], nvarsblock[i+1], 697 heurdata->linkingconss, heurdata->nlinking, i, success) ); 698 if( !(*success) ) 699 break; 700 } 701 702 SCIPfreeBufferArray(scip, &nconssblock); 703 SCIPfreeBufferArray(scip, &nvarsblock); 704 705 return SCIP_OKAY; 706 } 707 708 /** rounds partition for one linking constraint to integer value if variables and coefficients are integer 709 * 710 * changes only currentrhs/currentlhs 711 */ 712 static 713 SCIP_RETCODE roundPartition( 714 SCIP* scip, /**< SCIP data structure */ 715 LINKING* linking, /**< one linking data structure */ 716 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 717 SCIP_Bool roundbyrhs /**< round by right hand side? */ 718 ) 719 { 720 SCIP_Real* fracPart; 721 int* sorting; 722 int* isinteger; 723 SCIP_Real sumbefor; /* includes value at idx */ 724 SCIP_Real sumafter; 725 SCIP_Real diff; 726 int nnonintblocks; /* number of non integer blocks */ 727 int idx; 728 int b; 729 int i; 730 int k; 731 732 assert(scip != NULL); 733 assert(linking != NULL); 734 assert(blockproblem != NULL); 735 736 nnonintblocks = 0; 737 idx = 0; 738 739 SCIP_CALL( SCIPallocBufferArray(scip, &fracPart, linking->nblocks) ); 740 SCIP_CALL( SCIPallocBufferArray(scip, &sorting, linking->nblocks) ); 741 SCIP_CALL( SCIPallocBufferArray(scip, &isinteger, linking->nblocks) ); 742 743 /* get integer blocks and fractional parts */ 744 for( b = 0; b < linking->nblocks; b++ ) 745 { 746 SCIP* subscip; 747 SCIP_CONS* blockcons; 748 SCIP_VAR** blockvars; 749 SCIP_Real* blockvals; 750 int nblockvars; 751 int length; /* number of block variables without slack variables */ 752 SCIP_Bool success; 753 754 subscip = blockproblem[linking->blocknumbers[b]]->blockscip; 755 blockcons = linking->blockconss[b]; 756 sorting[b] = b; /* store current sorting to sort back */ 757 758 SCIP_CALL( SCIPgetConsNVars(subscip, blockcons, &nblockvars, &success) ); 759 assert(success); 760 SCIP_CALL( SCIPallocBufferArray(scip, &blockvars, nblockvars) ); 761 SCIP_CALL( SCIPallocBufferArray(scip, &blockvals, nblockvars) ); 762 763 SCIP_CALL( SCIPgetConsVars(subscip, blockcons, blockvars, nblockvars, &success) ); 764 assert(success); 765 SCIP_CALL( SCIPgetConsVals(subscip, blockcons, blockvals, nblockvars, &success) ); 766 assert(success); 767 768 /* get number of block variables in this constraint without slack variables */ 769 length = nblockvars - linking->nslacksperblock; 770 771 /* get blocks with integer value */ 772 isinteger[b] = 1; 773 for( i = 0; i < length; i++ ) 774 { 775 if( SCIPvarGetType(blockvars[i]) == SCIP_VARTYPE_CONTINUOUS || !SCIPisIntegral(scip, blockvals[i]) ) 776 { 777 isinteger[b] = 0; 778 nnonintblocks++; 779 break; 780 } 781 } 782 783 /* get fractional part of blockconstraints */ 784 if( roundbyrhs ) 785 fracPart[b] = linking->currentrhs[b] - floor(linking->currentrhs[b]); /* do not use SCIPfrac()! */ 786 else 787 fracPart[b] = linking->currentlhs[b] - floor(linking->currentlhs[b]); 788 789 SCIPfreeBufferArray(scip, &blockvals); 790 SCIPfreeBufferArray(scip, &blockvars); 791 } 792 793 /* sort non-integer blocks to the front */ 794 SCIPsortIntIntReal(isinteger, sorting, fracPart, linking->nblocks); 795 796 /* sort by fractional part */ 797 SCIPsortRealInt(fracPart, sorting, nnonintblocks); 798 SCIPsortRealInt(&fracPart[nnonintblocks], &sorting[nnonintblocks], linking->nblocks - nnonintblocks); 799 800 /* detect blocks for rounding down and rounding up: 801 * integer blocks with small fractional parts are rounded down 802 * integer blocks with big fractional parts are rounded up 803 */ 804 805 sumbefor = 0.0; 806 sumafter = 0.0; 807 808 for( i = 0; i < linking->nblocks - nnonintblocks; i++ ) 809 sumafter += 1 - fracPart[nnonintblocks + i]; 810 811 for( i = 0; i < linking->nblocks - nnonintblocks; i++ ) 812 { 813 sumbefor += fracPart[nnonintblocks + i]; 814 sumafter -= 1 - fracPart[nnonintblocks + i]; 815 816 if( sumbefor >= sumafter ) 817 { 818 for( k = 0; k <= i; k++ ) 819 fracPart[nnonintblocks + k] = -fracPart[nnonintblocks + k]; 820 821 for( k = i + 1; k < linking->nblocks - nnonintblocks; k++ ) 822 fracPart[nnonintblocks + k] = 1 - fracPart[nnonintblocks + k]; 823 824 idx = i; 825 break; 826 } 827 } 828 diff = sumbefor - sumafter; 829 assert(SCIPisGE(scip, diff, 0.0)); 830 831 /* add difference to last non integer block */ 832 for( i = nnonintblocks - 1; i >= 0; i-- ) 833 { 834 if( SCIPisGT(scip, diff, 0.0) ) 835 { 836 fracPart[i] = diff; 837 diff = 0; 838 } 839 else 840 fracPart[i] = 0; 841 } 842 843 /* add difference to last rounded down block if no non integer block exists */ 844 if( SCIPisGT(scip, diff, 0.0)) 845 { 846 assert(nnonintblocks == 0); 847 fracPart[idx] += diff; 848 } 849 850 /* sort back */ 851 SCIPsortIntReal(sorting, fracPart, linking->nblocks); 852 853 /* round partition 854 * if we have a ranged constraint, both sides get rounded in the same way 855 */ 856 for( b = 0; b < linking->nblocks; b++ ) 857 { 858 if( linking->hasrhs ) 859 linking->currentrhs[b] += fracPart[b]; 860 861 if( linking->haslhs ) 862 linking->currentlhs[b] += fracPart[b]; 863 } 864 865 SCIPfreeBufferArray(scip, &isinteger); 866 SCIPfreeBufferArray(scip, &sorting); 867 SCIPfreeBufferArray(scip, &fracPart); 868 869 return SCIP_OKAY; 870 } 871 872 /** calculates initial partition and sets rhs/lhs in blockproblems */ 873 static 874 SCIP_RETCODE initCurrent( 875 SCIP* scip, /**< SCIP data structure of main scip */ 876 LINKING** linkings, /**< array of linking data structures */ 877 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 878 SCIP_HEURTIMING heurtiming, /**< current point in the node solving process */ 879 int nlinking, /**< number of linking constraints */ 880 SCIP_Bool* success /**< pointer to store whether initialization was successful */ 881 ) 882 { 883 LINKING* linking; 884 SCIP_Real rhs; 885 SCIP_Real lhs; 886 SCIP_Real residualrhs; 887 SCIP_Real residuallhs; 888 SCIP_Real goalvalue; 889 int c; 890 int b; 891 892 assert(scip != NULL); 893 assert(linkings != NULL); 894 assert(blockproblem != NULL); 895 assert(nlinking > 0); 896 897 SCIPdebugMsg(scip, "initialize partition\n"); 898 899 /* for each linking constraint the rhs/lhs is split between the blocks */ 900 for( c = 0; c < nlinking; c++ ) 901 { 902 linking = linkings[c]; 903 rhs = SCIPconsGetRhs(scip, linking->linkingcons, success); 904 assert(*success); 905 lhs = SCIPconsGetLhs(scip, linking->linkingcons, success); 906 assert(*success); 907 residualrhs = rhs; 908 residuallhs = lhs; 909 910 /* LP value for each block plus part of remaining amount if sum is not equal to rhs/lhs */ 911 if( (heurtiming & SCIP_HEURTIMING_AFTERNODE) && (linking->hasrhs || linking->haslhs) ) 912 { 913 SCIP_Real sumrhs = 0.0; 914 SCIP_Real sumlhs = 0.0; 915 for( b = 0; b < linking->nblocks; b++ ) 916 { 917 SCIP_VAR** consvars; 918 SCIP_Real* consvals; 919 SCIP_CONS* cons; 920 int nconsvars; 921 SCIP_Real lpvalue; 922 int i; 923 924 /* get variables of linking constraint of one block */ 925 cons = linking->blockconss[b]; 926 SCIP_CALL( SCIPgetConsNVars(blockproblem[linking->blocknumbers[b]]->blockscip, cons, &nconsvars, success) ); 927 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nconsvars) ); 928 SCIP_CALL( SCIPallocBufferArray(scip, &consvals, nconsvars) ); 929 SCIP_CALL( SCIPgetConsVars(blockproblem[linking->blocknumbers[b]]->blockscip, cons, consvars, nconsvars, success) ); 930 SCIP_CALL( SCIPgetConsVals(blockproblem[linking->blocknumbers[b]]->blockscip, cons, consvals, nconsvars, success) ); 931 assert(success); 932 933 /* calculate value of partition variable in lp solution */ 934 lpvalue = 0.0; 935 for( i = 0; i < nconsvars - linking->nslacksperblock; i++ ) 936 { 937 SCIP_VAR* origvar; 938 SCIP_Real varlpvalue; 939 940 origvar = SCIPfindVar(scip, SCIPvarGetName(consvars[i])); 941 if( origvar == NULL ) /* e.g. variable is negated */ 942 { 943 *success = FALSE; 944 return SCIP_OKAY; 945 } 946 varlpvalue = SCIPvarGetLPSol(origvar); 947 lpvalue += varlpvalue * consvals[i]; 948 } 949 sumrhs += lpvalue; 950 sumlhs += lpvalue; 951 952 /* set currentrhs/lhs at lp value */ 953 if( linking->hasrhs ) 954 linking->currentrhs[b] = lpvalue; 955 if( linking->haslhs ) 956 linking->currentlhs[b] = lpvalue; 957 958 SCIPfreeBufferArray(scip, &consvars); 959 SCIPfreeBufferArray(scip, &consvals); 960 } 961 assert(SCIPisLE(scip, sumrhs, rhs)); 962 assert(SCIPisGE(scip, sumlhs, lhs)); 963 964 /* distribute remaining amount */ 965 if( !SCIPisEQ(scip, rhs, sumrhs) && linking->hasrhs ) 966 { 967 SCIP_Real diff; 968 SCIP_Real part; 969 SCIP_Real residual; 970 diff = rhs - sumrhs; 971 residual = 0.0; 972 973 for( b = 0; b < linking->nblocks; b++ ) 974 { 975 goalvalue = linking->currentrhs[b] + ( diff / linking->nblocks ) + residual; 976 part = MIN(goalvalue, linking->maxactivity[b]); 977 residual = goalvalue - part; 978 linking->currentrhs[b] = part; 979 } 980 if( !SCIPisZero(scip, residual) ) 981 linking->currentrhs[0] += residual; 982 } 983 if( !SCIPisEQ(scip, lhs, sumlhs) && linking->haslhs ) 984 { 985 SCIP_Real diff; 986 SCIP_Real part; 987 SCIP_Real residual; 988 diff = sumlhs - lhs; /* always positive*/ 989 residual = 0.0; 990 991 for( b = 0; b < linking->nblocks; b++ ) 992 { 993 goalvalue = linking->currentlhs[b] - ( diff / linking->nblocks ) + residual; 994 part = MAX(goalvalue, linking->minactivity[b]); 995 residual = goalvalue - part; 996 linking->currentlhs[b] = part; 997 } 998 if( !SCIPisZero(scip, residual) ) 999 linking->currentlhs[0] += residual; 1000 } 1001 } 1002 /* equal parts for each block with respect to minimal/maximal activity */ 1003 else if( linking->hasrhs || linking->haslhs ) 1004 { 1005 if( linking->hasrhs ) 1006 { 1007 for( b = 0; b < linking->nblocks; b++ ) 1008 { 1009 goalvalue = residualrhs / (linking->nblocks - b); 1010 linking->currentrhs[b] = MIN(MAX(goalvalue, linking->minactivity[b]), linking->maxactivity[b]); 1011 residualrhs -= linking->currentrhs[b]; 1012 } 1013 /* add residual partition to first block */ 1014 linking->currentrhs[0] += residualrhs; 1015 } 1016 if( linking->haslhs ) 1017 { 1018 for( b = 0; b < linking->nblocks; b++ ) 1019 { 1020 goalvalue = residuallhs / (linking->nblocks - b); 1021 linking->currentlhs[b] = MIN(MAX(goalvalue, linking->minactivity[b]), linking->maxactivity[b]); 1022 residuallhs -= linking->currentlhs[b]; 1023 } 1024 /* add residual partition to first block */ 1025 linking->currentlhs[0] += residuallhs; 1026 } 1027 } 1028 else 1029 { 1030 assert(linking->nblocks == 0 && !SCIPconsIsChecked(linking->linkingcons)); 1031 } 1032 1033 SCIP_CALL( roundPartition(scip, linking, blockproblem, linking->hasrhs) ); 1034 1035 /* set sides in blockproblem at initial partition */ 1036 for( b = 0; b < linking->nblocks; b++ ) 1037 { 1038 if( linking->hasrhs ) 1039 { 1040 SCIP_CALL( SCIPchgRhsLinear(blockproblem[linking->blocknumbers[b]]->blockscip, 1041 linking->blockconss[b], linking->currentrhs[b]) ); 1042 #ifdef SCIP_MORE_DEBUG 1043 SCIPdebugMsg(scip, "change rhs of %s in block %d to %f\n", 1044 SCIPconsGetName(linking->linkingcons), linking->blocknumbers[b], linking->currentrhs[b]); 1045 #endif 1046 } 1047 if( linking->haslhs ) 1048 { 1049 SCIP_CALL( SCIPchgLhsLinear(blockproblem[linking->blocknumbers[b]]->blockscip, 1050 linking->blockconss[b], linking->currentlhs[b]) ); 1051 #ifdef SCIP_MORE_DEBUG 1052 SCIPdebugMsg(scip, "change lhs of %s in block %d to %f\n", 1053 SCIPconsGetName(linking->linkingcons), linking->blocknumbers[b], linking->currentlhs[b]); 1054 #endif 1055 } 1056 } 1057 } 1058 1059 return SCIP_OKAY; 1060 } 1061 1062 /** calculates shift */ 1063 static 1064 SCIP_RETCODE calculateShift( 1065 SCIP* scip, /**< SCIP data structure of main scip */ 1066 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 1067 LINKING* linking, /**< one linking data structure */ 1068 SCIP_Real** shift, /**< pointer to store direction vector */ 1069 int* nviolatedblocksrhs, /**< pointer to store number of blocks which violate rhs */ 1070 int* nviolatedblockslhs, /**< pointer to store number of blocks which violate lhs */ 1071 SCIP_Bool* update /**< should we update the partition? */ 1072 ) 1073 { 1074 int* nonviolatedblocksrhs; 1075 int* nonviolatedblockslhs; 1076 SCIP_Real sumviols = 0.0; 1077 int v; 1078 1079 if( linking->hasrhs ) 1080 { 1081 SCIP_CALL( SCIPallocBufferArray(scip, &nonviolatedblocksrhs, linking->nblocks) ); 1082 } 1083 if( linking->haslhs ) 1084 { 1085 SCIP_CALL( SCIPallocBufferArray(scip, &nonviolatedblockslhs, linking->nblocks) ); 1086 } 1087 1088 /* get violated blocks and calculate shift */ 1089 for( v = 0; v < linking->nblocks; v++ ) 1090 { 1091 SCIP* subscip; 1092 SCIP_SOL* subsol; 1093 SCIP_Real slackval; 1094 1095 subscip = blockproblem[linking->blocknumbers[v]]->blockscip; 1096 subsol = SCIPgetBestSol(subscip); 1097 1098 /* if we have ranged constraints, the slack variables of the rhs are in front; 1099 * get slack variable of block; save violated blocks 1100 */ 1101 if( linking->hasrhs ) 1102 { 1103 slackval = SCIPgetSolVal(subscip, subsol, linking->slacks[v * linking->nslacksperblock]); 1104 1105 /* block is violated */ 1106 if( SCIPisPositive(scip, slackval) ) 1107 { 1108 (*nviolatedblocksrhs)++; 1109 1110 (*shift)[v] += slackval; 1111 sumviols += slackval; 1112 } 1113 else 1114 { 1115 nonviolatedblocksrhs[v - *nviolatedblocksrhs] = v; /*lint !e644*/ 1116 } 1117 } 1118 if( linking->haslhs ) 1119 { 1120 slackval = SCIPgetSolVal(subscip, subsol, linking->slacks[(v * linking->nslacksperblock) + linking->nslacksperblock - 1]); 1121 1122 /* block is violated */ 1123 if( SCIPisPositive(scip, slackval) ) 1124 { 1125 (*nviolatedblockslhs)++; 1126 1127 (*shift)[v] -= slackval; 1128 sumviols -= slackval; 1129 } 1130 else 1131 { 1132 nonviolatedblockslhs[v - *nviolatedblockslhs] = v; /*lint !e644*/ 1133 } 1134 } 1135 } 1136 1137 /* linking constraint is in no block violated or always violated -> do not update partition */ 1138 if( *nviolatedblocksrhs + *nviolatedblockslhs == 0 || 1139 linking->nblocks == *nviolatedblocksrhs || linking->nblocks == *nviolatedblockslhs ) 1140 { 1141 *update = FALSE; 1142 /* free memory */ 1143 if( linking->haslhs ) 1144 SCIPfreeBufferArray(scip, &nonviolatedblockslhs); 1145 if( linking->hasrhs ) 1146 SCIPfreeBufferArray(scip, &nonviolatedblocksrhs); 1147 return SCIP_OKAY; 1148 } 1149 1150 /* set values of non violated blocks */ 1151 if( SCIPisPositive(scip, sumviols) ) 1152 { 1153 /* rhs of original linking constraint is violated */ 1154 SCIP_Real residual = sumviols; 1155 SCIP_Real part; 1156 SCIP_Real shift_tmp; 1157 1158 assert(linking->hasrhs); 1159 assert(*nviolatedblocksrhs != 0); 1160 1161 /* substract from each non violated block the same amount with respect to minimal/maximal activity, 1162 * so that the shift is zero in sum 1163 */ 1164 for( v = 0; v < linking->nblocks - *nviolatedblocksrhs; v++ ) 1165 { 1166 part = linking->currentrhs[nonviolatedblocksrhs[v]] - residual/(linking->nblocks - *nviolatedblocksrhs - v); 1167 part = MIN(MAX(part, linking->minactivity[nonviolatedblocksrhs[v]]), linking->maxactivity[nonviolatedblocksrhs[v]]); 1168 shift_tmp = part - linking->currentrhs[nonviolatedblocksrhs[v]]; 1169 residual += shift_tmp; 1170 (*shift)[nonviolatedblocksrhs[v]] += shift_tmp; 1171 } 1172 if( !SCIPisZero(scip, residual) ) 1173 { 1174 /* assign residual to ... */ 1175 if( linking->nblocks - *nviolatedblocksrhs == 1 ) 1176 (*shift)[nonviolatedblocksrhs[0] == 0 ? 1 : 0] -= residual; /* first violated block */ 1177 else 1178 (*shift)[nonviolatedblocksrhs[0]] -= residual; /* first nonviolated block */ 1179 } 1180 } 1181 if( SCIPisNegative(scip, sumviols) ) 1182 { 1183 /* lhs of original linking constraint is violated */ 1184 SCIP_Real residual = sumviols; 1185 SCIP_Real part; 1186 SCIP_Real shift_tmp; 1187 1188 assert(linking->haslhs); 1189 assert(*nviolatedblockslhs != 0); 1190 1191 /* add to each non violated block the same amount with respect to minimal/maximal activity, 1192 * so that the shift is zero in sum 1193 */ 1194 for( v = 0; v < linking->nblocks - *nviolatedblockslhs; v++ ) 1195 { 1196 part = linking->currentlhs[nonviolatedblockslhs[v]] - residual/(linking->nblocks - *nviolatedblockslhs - v); 1197 part = MIN(MAX(part, linking->minactivity[nonviolatedblockslhs[v]]), linking->maxactivity[nonviolatedblockslhs[v]]); 1198 shift_tmp = part - linking->currentlhs[nonviolatedblockslhs[v]]; 1199 residual += shift_tmp; 1200 (*shift)[nonviolatedblockslhs[v]] += shift_tmp; 1201 } 1202 if( !SCIPisZero(scip, residual) ) 1203 { 1204 /* assign residual to ... */ 1205 if( linking->nblocks - *nviolatedblockslhs == 1 ) 1206 (*shift)[nonviolatedblockslhs[0] == 0 ? 1 : 0] -= residual; /* first violated block */ 1207 else 1208 (*shift)[nonviolatedblockslhs[0]] -= residual; /* first nonviolated block */ 1209 } 1210 } 1211 1212 *update = TRUE; 1213 1214 /* free memory */ 1215 if( linking->haslhs ) 1216 SCIPfreeBufferArray(scip, &nonviolatedblockslhs); 1217 if( linking->hasrhs ) 1218 SCIPfreeBufferArray(scip, &nonviolatedblocksrhs); 1219 1220 return SCIP_OKAY; 1221 } 1222 1223 /** update partition */ 1224 static 1225 SCIP_RETCODE updatePartition( 1226 SCIP* scip, /**< SCIP data structure of main scip */ 1227 LINKING** linkings, /**< linking data structure */ 1228 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 1229 int** nviolatedblocksrhs, /**< pointer to store number of blocks which violate rhs */ 1230 int** nviolatedblockslhs, /**< pointer to store number of blocks which violate lhs */ 1231 int nlinking, /**< number of linking constraints */ 1232 int nblocks, /**< number of blocks */ 1233 int iteration, /**< number of iteration */ 1234 SCIP_Bool* oneupdate /**< is at least partition of one constraint updated? */ 1235 ) 1236 { 1237 SCIP_Real* shift; /* direction vector */ 1238 int v; 1239 int c; 1240 SCIP_Bool update; /* is partition of current constraint updated? */ 1241 1242 assert(scip != NULL); 1243 assert(linkings != NULL); 1244 assert(blockproblem != NULL); 1245 assert(iteration >= 0); 1246 assert(!*oneupdate); 1247 1248 SCIP_CALL( SCIPallocBufferArray(scip, &shift, nblocks) ); /* allocate memory only once */ 1249 1250 for( c = 0; c < nlinking; c++ ) 1251 { 1252 LINKING* linking; /* one linking data structure */ 1253 1254 linking = linkings[c]; 1255 (*nviolatedblocksrhs)[c] = 0; 1256 (*nviolatedblockslhs)[c] = 0; 1257 update = TRUE; 1258 BMSclearMemoryArray(shift, linking->nblocks); 1259 assert(nblocks >= linking->nblocks); 1260 1261 SCIP_CALL( calculateShift(scip, blockproblem, linking, &shift, &(*nviolatedblocksrhs)[c], &(*nviolatedblockslhs)[c], &update) ); 1262 1263 if( !update ) 1264 continue; 1265 1266 *oneupdate = TRUE; 1267 1268 #ifdef SCIP_DEBUG 1269 SCIP_Real sum = 0.0; 1270 /* check sum of shift; must be zero */ 1271 for( v = 0; v < linking->nblocks; v++ ) 1272 sum += shift[v]; 1273 assert(SCIPisFeasZero(scip, sum)); 1274 #endif 1275 1276 /* add shift to both sides */ 1277 for( v = 0; v < linking->nblocks; v++) 1278 { 1279 if( linking->hasrhs ) 1280 linking->currentrhs[v] += shift[v]; 1281 1282 if( linking->haslhs ) 1283 linking->currentlhs[v] += shift[v]; 1284 } 1285 1286 SCIP_CALL( roundPartition(scip, linking, blockproblem, ((linking->hasrhs && ((*nviolatedblocksrhs)[c] != 0)) || !linking->haslhs)) ); 1287 1288 /* set sides in blockproblems to new partition */ 1289 for( v = 0; v < linking->nblocks; v++ ) 1290 { 1291 SCIP* subscip; 1292 subscip = blockproblem[linking->blocknumbers[v]]->blockscip; 1293 1294 if( linking->hasrhs ) 1295 { 1296 SCIP_CALL( SCIPchgRhsLinear(subscip, linking->blockconss[v], linking->currentrhs[v]) ); 1297 } 1298 if( linking->haslhs ) 1299 { 1300 SCIP_CALL( SCIPchgLhsLinear(subscip, linking->blockconss[v], linking->currentlhs[v]) ); 1301 } 1302 } 1303 } 1304 1305 /* free memory */ 1306 SCIPfreeBufferArray(scip, &shift); 1307 1308 return SCIP_OKAY; 1309 } 1310 1311 /** update penalty parameters lambda 1312 * 1313 * if a linking constraint is violated two times in succession, the corresponding penalty parameter is increased in each block 1314 */ 1315 static 1316 SCIP_RETCODE updateLambda( 1317 SCIP* scip, /**< SCIP data structure */ 1318 SCIP_HEURDATA* heurdata, /**< heuristic data */ 1319 LINKING** linkings, /**< array of linking data structures */ 1320 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 1321 int* nviolatedblocksrhs, /**< number of blocks which violate rhs */ 1322 int* nviolatedblockslhs, /**< number of blocks which violate lhs */ 1323 int nlinking /**< number of linking constraints */ 1324 ) 1325 { 1326 SCIP_VAR* slackvar; 1327 SCIP_Real old_obj; 1328 SCIP_Real new_obj; 1329 int c; 1330 int b; 1331 1332 assert(scip != NULL); 1333 assert(linkings != NULL); 1334 assert(blockproblem != NULL); 1335 1336 for( c = 0; c < nlinking; c++ ) 1337 { 1338 assert(nviolatedblocksrhs[c] >= 0); 1339 assert(nviolatedblockslhs[c] >= 0); 1340 1341 /* skip constraint if it is not violated */ 1342 if( nviolatedblocksrhs[c] + nviolatedblockslhs[c] == 0 ) 1343 { 1344 linkings[c]->lastviolations = 0; /* reset flag */ 1345 continue; 1346 } 1347 1348 /* add number of violated blocks multiplied with parameter "penalty" to lambda (initial value is 1) */ 1349 for( b = 0; b < linkings[c]->nblocks; b++ ) 1350 { 1351 SCIP* subscip; 1352 subscip = blockproblem[linkings[c]->blocknumbers[b]]->blockscip; 1353 assert(subscip != NULL); 1354 1355 if( linkings[c]->hasrhs && (nviolatedblocksrhs[c] >= 1) && (linkings[c]->lastviolations >= 1) ) 1356 { 1357 slackvar = linkings[c]->slacks[b * linkings[c]->nslacksperblock]; 1358 old_obj = SCIPvarGetObj(slackvar); 1359 new_obj = old_obj + heurdata->penalty * nviolatedblocksrhs[c]; 1360 1361 SCIP_CALL( SCIPchgVarObj(subscip, slackvar, new_obj) ); 1362 } 1363 if( linkings[c]->haslhs && (nviolatedblockslhs[c] >= 1) && (linkings[c]->lastviolations >= 1) ) 1364 { 1365 slackvar = linkings[c]->slacks[b * linkings[c]->nslacksperblock + linkings[c]->nslacksperblock - 1]; 1366 old_obj = SCIPvarGetObj(slackvar); 1367 new_obj = old_obj + heurdata->penalty * nviolatedblockslhs[c]; 1368 1369 SCIP_CALL( SCIPchgVarObj(subscip, slackvar, new_obj) ); 1370 } 1371 } 1372 1373 /* update number of violations in the last iterations */ 1374 linkings[c]->lastviolations += 1; 1375 } 1376 1377 return SCIP_OKAY; 1378 } 1379 1380 /** computes feasible solution from last stored solution for each block and adds it to the solution storage */ 1381 static 1382 SCIP_RETCODE reuseSolution( 1383 LINKING** linkings, /**< array of linking data structures */ 1384 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 1385 int nblocks /**< number of blocks */ 1386 ) 1387 { 1388 SCIP_SOL** sols; 1389 SCIP_SOL* sol; /* solution of block that will be repaired */ 1390 SCIP_SOL* newsol; 1391 SCIP_VAR** blockvars; 1392 SCIP_Real* blockvals; 1393 int nsols; 1394 int nvars; 1395 int b; 1396 int c; 1397 int i; 1398 SCIP_Bool success; 1399 1400 assert(linkings != NULL); 1401 assert(blockproblem != NULL); 1402 1403 for( b = 0; b < nblocks; b++ ) 1404 { 1405 SCIP* subscip; 1406 1407 subscip = blockproblem[b]->blockscip; 1408 nsols = SCIPgetNSols(subscip); 1409 1410 /* no solution in solution candidate storage found */ 1411 if( nsols == 0 ) 1412 return SCIP_OKAY; 1413 1414 /* take last solution */ 1415 sols = SCIPgetSols(subscip); 1416 sol = sols[nsols - 1]; 1417 1418 /* copy the solution */ 1419 nvars = SCIPgetNVars(subscip); 1420 blockvars = SCIPgetVars(subscip); 1421 SCIP_CALL( SCIPallocBufferArray(subscip, &blockvals, nvars) ); 1422 SCIP_CALL( SCIPgetSolVals(subscip, sol, nvars, blockvars, blockvals) ); 1423 SCIP_CALL( SCIPcreateOrigSol(subscip, &newsol, NULL) ); 1424 SCIP_CALL( SCIPsetSolVals(subscip, newsol, nvars, blockvars, blockvals) ); 1425 1426 /* correct each coupling constraint: 1427 * lhs <= orig_var - z_r + z_l <= rhs 1428 * adapt slack variables so that constraint is feasible 1429 */ 1430 for( c = 0; c < blockproblem[b]->nlinking; c++ ) 1431 { 1432 LINKING* linking; /* linking data structure of this constraint */ 1433 SCIP_VAR* rvar; /* slack variable z_r */ 1434 SCIP_VAR* lvar; /* slack variable z_l */ 1435 SCIP_Real rval; /* value of slack variable z_r */ 1436 SCIP_Real lval; /* value of slack variable z_l */ 1437 SCIP_Real activitycons; /* activity of constraint*/ 1438 SCIP_Real activity; /* activity of constraint without slack variables */ 1439 SCIP_Real rhs; /* current right hand side */ 1440 SCIP_Real lhs; /* current left hand side */ 1441 1442 linking = linkings[blockproblem[b]->linkingindices[c]]; 1443 rhs = SCIPgetRhsLinear(subscip, blockproblem[b]->linkingconss[c]); 1444 lhs = SCIPgetLhsLinear(subscip, blockproblem[b]->linkingconss[c]); 1445 assert(SCIPisGE(subscip, rhs, lhs)); 1446 1447 activitycons = SCIPgetActivityLinear(subscip, blockproblem[b]->linkingconss[c], sol); 1448 1449 /* get slack variables and subtract their value from the activity; 1450 * calculate and set values of slack variables 1451 */ 1452 for( i = 0; i < linking->nblocks; i++ ) 1453 { 1454 if( linking->blocknumbers[i] == b ) 1455 { 1456 if( linking->hasrhs && linking->haslhs ) 1457 { 1458 rvar = linking->slacks[2 * i]; 1459 lvar = linking->slacks[2 * i + 1]; 1460 rval = SCIPgetSolVal(subscip, sol, rvar); 1461 lval = SCIPgetSolVal(subscip, sol, lvar); 1462 activity = activitycons + rval - lval; 1463 SCIP_CALL( SCIPsetSolVal(subscip, newsol, rvar, MAX(0.0, activity - rhs)) ); 1464 SCIP_CALL( SCIPsetSolVal(subscip, newsol, lvar, MAX(0.0, lhs - activity)) ); 1465 } 1466 else if( linking->hasrhs ) 1467 { 1468 rvar = linking->slacks[i]; 1469 rval = SCIPgetSolVal(subscip, sol, rvar); 1470 activity = activitycons + rval; 1471 SCIP_CALL( SCIPsetSolVal(subscip, newsol, rvar, MAX(0.0, activity - rhs)) ); 1472 } 1473 else /* linking->haslhs */ 1474 { 1475 assert(linking->haslhs); 1476 lvar = linking->slacks[i]; 1477 lval = SCIPgetSolVal(subscip, sol, lvar); 1478 activity = activitycons - lval; 1479 SCIP_CALL( SCIPsetSolVal(subscip, newsol, lvar, MAX(0.0, lhs - activity)) ); 1480 } 1481 break; 1482 } 1483 } 1484 } 1485 1486 SCIPdebugMsg(subscip, "Try adding solution with objective value %.2f\n", SCIPgetSolOrigObj(subscip, newsol)); 1487 SCIP_CALL( SCIPaddSolFree(subscip, &newsol, &success) ); 1488 1489 if( !success ) 1490 SCIPdebugMsg(subscip, "Correcting solution failed\n"); /* maybe not better than old solutions */ 1491 else 1492 SCIPdebugMsg(subscip, "Correcting solution successful\n"); 1493 1494 SCIPfreeBufferArray(subscip, &blockvals); 1495 } 1496 1497 return SCIP_OKAY; 1498 } 1499 1500 /** reoptimizes the heuristic solution with original objective function */ 1501 static 1502 SCIP_RETCODE reoptimize( 1503 SCIP* scip, /**< SCIP data structure */ 1504 SCIP_HEUR* heur, /**< pointer to heuristic */ 1505 SCIP_SOL* sol, /**< heuristic solution */ 1506 BLOCKPROBLEM** blockproblem, /**< array of blockproblem data structures */ 1507 int nblocks, /**< number of blockproblems */ 1508 SCIP_Bool limits, /**< should strict limits be set? */ 1509 SCIP_SOL** newsol, /**< pointer to store improved solution */ 1510 SCIP_Bool* success /**< pointer to store whether reoptimization was successful */ 1511 ) 1512 { 1513 SCIP_Real time; 1514 SCIP_Real timesubscip; 1515 SCIP_Bool check; 1516 int b; 1517 int v; 1518 1519 assert(scip != NULL); 1520 assert(heur != NULL); 1521 assert(sol != NULL); 1522 assert(blockproblem != NULL); 1523 1524 *success = FALSE; 1525 check = FALSE; 1526 1527 /* get used time */ 1528 timesubscip = SCIPgetTotalTime(blockproblem[0]->blockscip); 1529 1530 /* for each blockproblem: 1531 * - change back to original objective function 1532 * - fix slack variables to zero 1533 * - set limits and solve problem 1534 */ 1535 for( b = 0; b < nblocks; b++ ) 1536 { 1537 SCIP* subscip; 1538 SCIP_VAR** blockvars; 1539 SCIP_Real infvalue; 1540 int nvars; 1541 int nsols; 1542 1543 subscip = blockproblem[b]->blockscip; 1544 blockvars = SCIPgetOrigVars(subscip); 1545 nvars = SCIPgetNOrigVars(subscip); 1546 nsols = 0; 1547 1548 /* in order to change objective function */ 1549 SCIP_CALL( SCIPfreeTransform(subscip) ); 1550 1551 /* change back to original objective function */ 1552 for( v = 0; v < blockproblem[b]->nblockvars; v++ ) 1553 { 1554 SCIP_CALL( SCIPchgVarObj(subscip, blockvars[v], blockproblem[b]->origobj[v]) ); 1555 } 1556 1557 /* fix slack variables to zero */ 1558 for( v = blockproblem[b]->nblockvars; v < nvars; v++ ) 1559 { 1560 SCIP_CALL( SCIPchgVarUb(subscip, blockvars[v], 0.0) ); 1561 SCIP_CALL( SCIPchgVarLb(subscip, blockvars[v], 0.0) ); 1562 } 1563 1564 /* do not abort subproblem on CTRL-C */ 1565 SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", FALSE) ); 1566 1567 /* forbid recursive call of heuristics and separators solving sub-SCIPs */ 1568 SCIP_CALL( SCIPsetSubscipsOff(subscip, TRUE) ); 1569 1570 #ifdef SCIP_DEBUG 1571 /* for debugging, enable full output */ 1572 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) ); 1573 SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 100000000) ); 1574 #else 1575 /* disable statistic timing inside sub SCIP and output to console */ 1576 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) ); 1577 SCIP_CALL( SCIPsetBoolParam(subscip, "timing/statistictiming", FALSE) ); 1578 #endif 1579 1580 /* disable cutting plane separation */ 1581 SCIP_CALL( SCIPsetSeparating(subscip, SCIP_PARAMSETTING_OFF, TRUE) ); 1582 1583 /* disable expensive presolving */ 1584 SCIP_CALL( SCIPsetPresolving(subscip, SCIP_PARAMSETTING_FAST, TRUE) ); 1585 1586 /* disable expensive techniques */ 1587 SCIP_CALL( SCIPsetIntParam(subscip, "misc/usesymmetry", 0) ); 1588 1589 /* speed up sub-SCIP by not checking dual LP feasibility */ 1590 SCIP_CALL( SCIPsetBoolParam(subscip, "lp/checkdualfeas", FALSE) ); 1591 1592 /* copy value for infinity */ 1593 SCIP_CALL( SCIPgetRealParam(scip, "numerics/infinity", &infvalue) ); 1594 SCIP_CALL( SCIPsetRealParam(subscip, "numerics/infinity", infvalue) ); 1595 1596 /* set limits; do not use more time in each subproblem than the heuristic has already used for first solution */ 1597 SCIP_CALL( SCIPcopyLimits(scip, subscip) ); 1598 if( limits ) 1599 { 1600 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", 1LL) ); 1601 SCIP_CALL( SCIPgetRealParam(subscip, "limits/time", &time) ); 1602 if( timesubscip < time - 1.0 ) 1603 { 1604 SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timesubscip + 1.0) ); 1605 } 1606 SCIP_CALL( SCIPtransformProb(subscip) ); 1607 nsols = SCIPgetNSols(subscip); 1608 /* one improving solution */ 1609 SCIP_CALL( SCIPsetIntParam(subscip, "limits/bestsol", nsols + 1) ); 1610 } 1611 else 1612 { 1613 /* only 50% of remaining time */ 1614 SCIP_CALL( SCIPgetRealParam(subscip, "limits/time", &time) ); 1615 SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", time * 0.5) ); 1616 /* set big node limits */ 1617 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", 1000LL) ); 1618 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", 100LL) ); 1619 } 1620 1621 /* reoptimize problem */ 1622 SCIP_CALL_ABORT( SCIPsolve(subscip) ); 1623 1624 if( SCIPgetNSols(subscip) == 0 ) 1625 { 1626 /* we found no solution */ 1627 return SCIP_OKAY; 1628 } 1629 else if( SCIPgetStatus(subscip) == SCIP_STATUS_BESTSOLLIMIT || 1630 SCIPgetStatus(subscip) == SCIP_STATUS_OPTIMAL || 1631 SCIPgetNSols(subscip) > nsols ) 1632 { 1633 check = TRUE; 1634 } 1635 } 1636 1637 if( !check ) 1638 { 1639 /* we have no better solution */ 1640 return SCIP_OKAY; 1641 } 1642 1643 /* create sol of main scip */ 1644 SCIP_CALL( SCIPcreateSol(scip, newsol, heur) ); 1645 1646 /* copy solution to main scip */ 1647 for( b = 0; b < nblocks; b++ ) 1648 { 1649 SCIP_SOL* blocksol; 1650 SCIP_VAR** blockvars; 1651 SCIP_Real* blocksolvals; 1652 int nblockvars; 1653 1654 /* get solution of block variables (without slack variables) */ 1655 blocksol = SCIPgetBestSol(blockproblem[b]->blockscip); 1656 blockvars = SCIPgetOrigVars(blockproblem[b]->blockscip); 1657 nblockvars = blockproblem[b]->nblockvars; 1658 1659 SCIP_CALL( SCIPallocBufferArray(scip, &blocksolvals, nblockvars) ); 1660 SCIP_CALL( SCIPgetSolVals(blockproblem[b]->blockscip, blocksol, nblockvars, blockvars, blocksolvals) ); 1661 1662 for( v = 0; v < nblockvars; v++ ) 1663 { 1664 SCIP_VAR* origvar; 1665 1666 origvar = SCIPfindVar(scip, SCIPvarGetName(blockvars[v])); 1667 SCIP_CALL( SCIPsetSolVal(scip, *newsol, origvar, blocksolvals[v]) ); 1668 } 1669 1670 SCIPfreeBufferArray(scip, &blocksolvals); 1671 } 1672 1673 *success = TRUE; 1674 1675 return SCIP_OKAY; 1676 } 1677 1678 1679 /* ---------------- Callback methods of event handler ---------------- */ 1680 1681 /* exec the event handler 1682 * 1683 * interrupt solution process of sub-SCIP if dual bound is greater than zero and a solution is available 1684 */ 1685 static 1686 SCIP_DECL_EVENTEXEC(eventExecDps) 1687 { 1688 assert(eventhdlr != NULL); 1689 assert(eventdata != NULL); 1690 assert(strcmp(SCIPeventhdlrGetName(eventhdlr), EVENTHDLR_NAME) == 0); 1691 assert(event != NULL); 1692 assert(SCIPeventGetType(event) & SCIP_EVENTTYPE_LPSOLVED); 1693 1694 SCIPdebugMsg(scip, "dual bound: %0.2f\n", SCIPgetDualbound(scip)); 1695 1696 if( SCIPisFeasGT(scip, SCIPgetDualbound(scip), 0.0) && SCIPgetNSols(scip) >= 1 ) 1697 { 1698 SCIPdebugMsg(scip, "DPS: interrupt subscip\n"); 1699 SCIP_CALL( SCIPinterruptSolve(scip) ); 1700 } 1701 1702 return SCIP_OKAY; 1703 } 1704 1705 1706 /* 1707 * Callback methods of primal heuristic 1708 */ 1709 1710 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 1711 static 1712 SCIP_DECL_HEURCOPY(heurCopyDps) 1713 { /*lint --e{715}*/ 1714 assert(scip != NULL); 1715 assert(heur != NULL); 1716 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 1717 1718 /* call inclusion method of primal heuristic */ 1719 SCIP_CALL( SCIPincludeHeurDps(scip) ); 1720 1721 return SCIP_OKAY; 1722 } 1723 1724 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 1725 static 1726 SCIP_DECL_HEURFREE(heurFreeDps) 1727 { /*lint --e{715}*/ 1728 SCIP_HEURDATA* heurdata; 1729 1730 assert(heur != NULL); 1731 assert(scip != NULL); 1732 1733 /* free heuristic data */ 1734 heurdata = SCIPheurGetData(heur); 1735 assert(heurdata != NULL); 1736 1737 SCIPfreeBlockMemory(scip, &heurdata); 1738 SCIPheurSetData(heur, NULL); 1739 1740 return SCIP_OKAY; 1741 } 1742 1743 /** execution method of primal heuristic */ 1744 static 1745 SCIP_DECL_HEUREXEC(heurExecDps) 1746 { /*lint --e{715}*/ 1747 SCIP_DECOMP** alldecomps; 1748 SCIP_DECOMP* decomp; 1749 SCIP_DECOMP* assigneddecomp; 1750 SCIP_VAR** vars; 1751 SCIP_CONS** conss; 1752 SCIP_VAR** sortedvars; 1753 SCIP_CONS** sortedconss; 1754 SCIP_HEURDATA* heurdata; 1755 BLOCKPROBLEM** blockproblem; 1756 LINKING** linkings; 1757 int* sortedvarlabels; 1758 int* sortedconslabels; 1759 SCIP_EVENTHDLR* eventhdlr; /* event handler */ 1760 SCIP_Real memory; /* in MB */ 1761 SCIP_Real timelimit; 1762 SCIP_Real allslacksval; 1763 SCIP_Real blocksolval; 1764 SCIP_STATUS status; 1765 SCIP_Bool avoidmemout; 1766 SCIP_Bool disablemeasures; 1767 int maxgraphedge; 1768 int ndecomps; 1769 int nvars; 1770 int nconss; 1771 int nblocks; 1772 SCIP_Bool success; 1773 int b; 1774 int c; 1775 int k; 1776 1777 assert( heur != NULL ); 1778 assert( scip != NULL ); 1779 assert( result != NULL ); 1780 1781 heurdata = SCIPheurGetData(heur); 1782 assert(heurdata != NULL); 1783 1784 assigneddecomp = NULL; 1785 blockproblem = NULL; 1786 linkings = NULL; 1787 eventhdlr = NULL; 1788 1789 *result = SCIP_DIDNOTRUN; 1790 1791 if( !((heurtiming & SCIP_HEURTIMING_BEFORENODE) && (heurdata->timing == 0 || heurdata->timing == 2)) 1792 && !((heurtiming & SCIP_HEURTIMING_AFTERNODE) && (heurdata->timing == 1 || heurdata->timing == 2)) ) 1793 { 1794 return SCIP_OKAY; 1795 } 1796 1797 /* call heuristic only once */ 1798 if( SCIPheurGetNCalls(heur) > 0 ) 1799 return SCIP_OKAY; 1800 1801 /* -------------------------------------------------------------------- */ 1802 SCIPdebugMsg(scip, "initialize dps heuristic\n"); 1803 1804 /* take the first transformed decomposition */ 1805 SCIPgetDecomps(scip, &alldecomps, &ndecomps, FALSE); 1806 if( ndecomps == 0 ) 1807 return SCIP_OKAY; 1808 1809 decomp = alldecomps[0]; 1810 assert(decomp != NULL); 1811 SCIPdebugMsg(scip, "First transformed decomposition is selected\n"); 1812 1813 nblocks = SCIPdecompGetNBlocks(decomp); 1814 nconss = SCIPgetNConss(scip); 1815 nvars = SCIPgetNVars(scip); 1816 1817 /* if problem has no constraints, no variables or less than two blocks, return */ 1818 if( nconss == 0 || nvars == 0 || nblocks <= 1 ) 1819 { 1820 SCIPdebugMsg(scip, "problem has no constraints, no variables or less than two blocks\n"); 1821 return SCIP_OKAY; 1822 } 1823 1824 /* estimate required memory for all blocks and terminate if not enough memory is available */ 1825 SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &memory) ); 1826 SCIP_CALL( SCIPgetBoolParam(scip, "misc/avoidmemout", &avoidmemout) ); 1827 if( avoidmemout && (((SCIPgetMemUsed(scip) + SCIPgetMemExternEstim(scip))/1048576.0) * (nblocks/4.0 + 2) >= memory) ) 1828 { 1829 SCIPdebugMsg(scip, "The estimated memory usage for %d blocks is too large.\n", nblocks); 1830 return SCIP_OKAY; 1831 } 1832 1833 /* we do not need the block decomposition graph and expensive measures of the decomposition statistics */ 1834 SCIP_CALL( SCIPgetIntParam(scip, "decomposition/maxgraphedge", &maxgraphedge) ); 1835 if( !SCIPisParamFixed(scip, "decomposition/maxgraphedge") ) 1836 { 1837 SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", 0) ); 1838 } 1839 SCIP_CALL( SCIPgetBoolParam(scip, "decomposition/disablemeasures", &disablemeasures) ); 1840 if( !SCIPisParamFixed(scip, "decomposition/disablemeasures") ) 1841 { 1842 SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", TRUE) ); 1843 } 1844 1845 vars = SCIPgetVars(scip); 1846 conss = SCIPgetConss(scip); 1847 SCIP_CALL( SCIPduplicateBufferArray(scip, &sortedvars, vars, nvars) ); 1848 SCIP_CALL( SCIPduplicateBufferArray(scip, &sortedconss, conss, nconss) ); 1849 SCIP_CALL( SCIPallocBufferArray(scip, &sortedvarlabels, nvars) ); 1850 SCIP_CALL( SCIPallocBufferArray(scip, &sortedconslabels, nconss) ); 1851 1852 /* get labels and sort in increasing order */ 1853 SCIPdecompGetVarsLabels(decomp, sortedvars, sortedvarlabels, nvars); 1854 SCIPdecompGetConsLabels(decomp, sortedconss, sortedconslabels, nconss); 1855 SCIPsortIntPtr(sortedvarlabels, (void**)sortedvars, nvars); 1856 SCIPsortIntPtr(sortedconslabels, (void**)sortedconss, nconss); 1857 1858 if( sortedvarlabels[0] == SCIP_DECOMP_LINKVAR ) 1859 { 1860 /* create new decomposition; don't change the decompositions in the decompstore */ 1861 SCIP_CALL( SCIPdecompCreate(&assigneddecomp, SCIPblkmem(scip), nblocks, SCIPdecompIsOriginal(decomp), SCIPdecompUseBendersLabels(decomp)) ); 1862 1863 SCIP_CALL( assignLinking(scip, assigneddecomp, sortedvars, sortedconss, sortedvarlabels, sortedconslabels, nvars, nconss, SCIPdecompGetNBorderVars(decomp)) ); 1864 assert(SCIPdecompGetNBlocks(decomp) >= SCIPdecompGetNBlocks(assigneddecomp)); 1865 decomp = assigneddecomp; 1866 1867 /* number of blocks can get smaller */ 1868 nblocks = SCIPdecompGetNBlocks(decomp); 1869 } 1870 else 1871 { 1872 /* The decomposition statistics were computed during transformation of the decomposition store. 1873 * Since propagators can have changed the number of constraints/variables, 1874 * the statistics are no longer up-to-date and have to be recomputed. 1875 */ 1876 SCIP_CALL( SCIPcomputeDecompStats(scip, decomp, TRUE) ); 1877 nblocks = SCIPdecompGetNBlocks(decomp); 1878 } 1879 1880 /* reset parameters */ 1881 SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", maxgraphedge) ); 1882 SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", disablemeasures) ); 1883 1884 #ifdef SCIP_DEBUG 1885 char buffer[SCIP_MAXSTRLEN]; 1886 SCIPdebugMsg(scip, "DPS used decomposition:\n%s\n", SCIPdecompPrintStats(decomp, buffer)); 1887 #endif 1888 1889 /* check properties of used decomposition */ 1890 if( heurdata->maxlinkscore != 1.0 ) 1891 { 1892 SCIP_Real linkscore; 1893 int nlinkconss; 1894 1895 nlinkconss = SCIPdecompGetNBorderConss(decomp); 1896 1897 linkscore = (SCIP_Real)nlinkconss / (SCIP_Real)nconss; 1898 1899 if( linkscore > heurdata->maxlinkscore ) 1900 { 1901 SCIPdebugMsg(scip, "decomposition has not the required properties\n"); 1902 goto TERMINATE; 1903 } 1904 } 1905 1906 if( sortedvarlabels[0] == SCIP_DECOMP_LINKVAR || 1907 sortedconslabels[0] != SCIP_DECOMP_LINKCONS || 1908 nblocks <= 1 ) 1909 { 1910 SCIPdebugMsg(scip, "Problem has linking variables or no linking constraints or less than two blocks\n"); 1911 goto TERMINATE; 1912 } 1913 1914 /* initialize heurdata */ 1915 heurdata->linkingconss = sortedconss; 1916 heurdata->nlinking = SCIPdecompGetNBorderConss(decomp); 1917 heurdata->nblocks = nblocks; 1918 1919 /* allocate memory for blockproblems and initialize partially */ 1920 SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem, nblocks) ); 1921 for( b = 0; b < nblocks; b++ ) 1922 { 1923 SCIP_CALL( SCIPallocBlockMemory(scip, &(blockproblem[b])) ); /*lint !e866*/ 1924 SCIP_CALL( createSubscip(scip, &blockproblem[b]->blockscip) ); 1925 1926 SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->linkingconss, heurdata->nlinking) ); 1927 SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->linkingindices, heurdata->nlinking) ); 1928 SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->slackvars, heurdata->nlinking * 2) ); /* maximum two slacks per linking constraint */ 1929 SCIP_CALL( SCIPallocBufferArray(scip, &blockproblem[b]->origobj, nvars) ); 1930 blockproblem[b]->nblockvars = 0; 1931 blockproblem[b]->nlinking = 0; 1932 blockproblem[b]->nslackvars = 0; 1933 } 1934 1935 /* allocate memory for linkings and initialize partially */ 1936 SCIP_CALL( SCIPallocBufferArray(scip, &linkings, heurdata->nlinking) ); 1937 for( c = 0; c < heurdata->nlinking; c++ ) 1938 { 1939 SCIP_CALL( SCIPallocBlockMemory(scip, &(linkings[c])) ); /*lint !e866*/ 1940 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->blockconss, heurdata->nblocks) ); 1941 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->slacks, heurdata->nblocks*2) ); /* maximum two slacks per block */ 1942 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->blocknumbers, heurdata->nblocks) ); 1943 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->minactivity, heurdata->nblocks) ); 1944 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->maxactivity, heurdata->nblocks) ); 1945 1946 linkings[c]->linkingcons = heurdata->linkingconss[c]; 1947 linkings[c]->currentrhs = NULL; 1948 linkings[c]->currentlhs = NULL; 1949 linkings[c]->nblocks = 0; 1950 linkings[c]->nslacks = 0; 1951 linkings[c]->nslacksperblock = 0; 1952 linkings[c]->lastviolations = 0; 1953 linkings[c]->hasrhs = FALSE; 1954 linkings[c]->haslhs = FALSE; 1955 } 1956 1957 SCIP_CALL( createAndSplitProblem(scip, heurdata, decomp, blockproblem, linkings, sortedvars, sortedconss, &success) ); 1958 if( !success ) 1959 { 1960 SCIPdebugMsg(scip, "Create and split problem failed\n"); 1961 goto TERMINATE; 1962 } 1963 1964 /* allocate memory for current partition*/ 1965 for( c = 0; c < heurdata->nlinking; c++ ) 1966 { 1967 if( linkings[c]->hasrhs ) 1968 { 1969 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->currentrhs, linkings[c]->nblocks ) ); 1970 } 1971 1972 if( linkings[c]->haslhs ) 1973 { 1974 SCIP_CALL( SCIPallocBufferArray(scip, &(linkings[c])->currentlhs, linkings[c]->nblocks ) ); 1975 } 1976 } 1977 1978 /* initialize partition */ 1979 SCIP_CALL( initCurrent(scip, linkings, blockproblem, heurtiming, heurdata->nlinking, &success) ); 1980 if( !success ) 1981 { 1982 SCIPdebugMsg(scip, "Initialization of partition failed\n"); 1983 goto TERMINATE; 1984 } 1985 1986 /* ------------------------------------------------------------------------ */ 1987 SCIPdebugMsg(scip, "Start heuristik DPS\n"); 1988 *result = SCIP_DIDNOTFIND; 1989 1990 for( k = 0; k < heurdata->maxit; k++ ) 1991 { 1992 /* do not exceed the timelimit */ 1993 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) ); 1994 if( (timelimit - SCIPgetSolvingTime(scip)) <= 0 ) 1995 goto TERMINATE; 1996 1997 /* solve the subproblems */ 1998 allslacksval = 0.0; 1999 for( b = 0; b < heurdata->nblocks; b++ ) 2000 { 2001 SCIP* subscip; 2002 subscip = blockproblem[b]->blockscip; 2003 2004 /* update time and memory limit of subproblem */ 2005 SCIP_CALL( SCIPcopyLimits(scip, subscip) ); 2006 2007 /* create event handler for LP events */ 2008 if( k == 0 ) 2009 { 2010 SCIP_CALL( SCIPincludeEventhdlrBasic(subscip, &eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC, eventExecDps, NULL) ); 2011 if( eventhdlr == NULL ) 2012 { 2013 SCIPerrorMessage("event handler for " HEUR_NAME " heuristic not found.\n"); 2014 return SCIP_PLUGINNOTFOUND; 2015 } 2016 } 2017 2018 /* catch LP events of sub-SCIP */ 2019 SCIP_CALL( SCIPtransformProb(subscip) ); 2020 SCIP_CALL( SCIPcatchEvent(subscip, SCIP_EVENTTYPE_LPSOLVED, eventhdlr, (SCIP_EVENTDATA*) heurdata, NULL) ); 2021 2022 SCIPdebugMsg(scip, "Solve blockproblem %d\n", b); 2023 SCIP_CALL_ABORT( SCIPsolve(subscip) ); 2024 2025 /* drop LP events of sub-SCIP */ 2026 SCIP_CALL( SCIPdropEvent(subscip, SCIP_EVENTTYPE_LPSOLVED, eventhdlr, (SCIP_EVENTDATA*) heurdata, -1) ); 2027 2028 /* get status and objective value if available */ 2029 status = SCIPgetStatus(subscip); 2030 if( status == SCIP_STATUS_INFEASIBLE ) 2031 { 2032 SCIPdebugMsg(scip, "Subproblem is infeasible\n"); 2033 goto TERMINATE; 2034 } 2035 else if( status == SCIP_STATUS_UNBOUNDED ) 2036 { 2037 SCIPdebugMsg(scip, "Subproblem is unbounded\n"); 2038 goto TERMINATE; 2039 } 2040 else if( SCIPgetNSols(subscip) >= 1 ) 2041 { 2042 blocksolval = SCIPgetPrimalbound(subscip); 2043 2044 if( status == SCIP_STATUS_TIMELIMIT && !SCIPisZero(scip, blocksolval) ) 2045 { 2046 SCIPdebugMsg(scip, "Subproblem reached timelimit without optimal solution\n"); 2047 goto TERMINATE; 2048 } 2049 SCIPdebugMsg(scip, "Solution value: %f\n", blocksolval); 2050 allslacksval += blocksolval; 2051 } 2052 else 2053 { 2054 SCIPdebugMsg(scip, "No subproblem solution available\n"); 2055 goto TERMINATE; 2056 } 2057 } 2058 2059 /* all slack variables are zero -> we found a feasible solution */ 2060 if( SCIPisZero(scip, allslacksval) ) 2061 { 2062 SCIP_SOL* newsol; 2063 2064 SCIPdebugMsg(scip, "Feasible solution found after %i iterations\n", k); 2065 2066 /* create new solution */ 2067 SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) ); 2068 for( b = 0; b < heurdata->nblocks; b++ ) 2069 { 2070 SCIP_SOL* blocksol; 2071 SCIP_VAR** blockvars; 2072 SCIP_Real* blocksolvals; 2073 int nblockvars; 2074 2075 /* get solution of block variables (without slack variables) */ 2076 blocksol = SCIPgetBestSol(blockproblem[b]->blockscip); 2077 blockvars = SCIPgetOrigVars(blockproblem[b]->blockscip); 2078 nblockvars = blockproblem[b]->nblockvars; 2079 2080 SCIP_CALL( SCIPallocBufferArray(scip, &blocksolvals, nblockvars) ); 2081 SCIP_CALL( SCIPgetSolVals(blockproblem[b]->blockscip, blocksol, nblockvars, blockvars, blocksolvals) ); 2082 2083 for( c = 0; c < nblockvars; c++ ) 2084 { 2085 SCIP_VAR* origvar; 2086 2087 origvar = SCIPfindVar(scip, SCIPvarGetName(blockvars[c])); 2088 SCIP_CALL( SCIPsetSolVal(scip, newsol, origvar, blocksolvals[c]) ); 2089 } 2090 2091 SCIPfreeBufferArray(scip, &blocksolvals); 2092 } 2093 2094 /* if reoptimization is activated, fix partition and reoptimize with original objective function */ 2095 if( heurdata->reoptimize ) 2096 { 2097 SCIP_SOL* improvedsol = NULL; 2098 SCIP_CALL( reoptimize(scip, heur, newsol, blockproblem, heurdata->nblocks, heurdata->reoptlimits, &improvedsol, &success) ); 2099 assert(improvedsol != NULL || success == FALSE); 2100 2101 if( success ) 2102 { 2103 SCIP_CALL( SCIPtrySolFree(scip, &improvedsol, TRUE, FALSE, TRUE, TRUE, TRUE, &success) ); 2104 if( success ) 2105 { 2106 SCIPdebugMsg(scip, "Reoptimizing solution successful\n"); 2107 *result = SCIP_FOUNDSOL; 2108 } 2109 } 2110 } 2111 2112 /* if reoptimization is turned off or reoptimization found no solution, try initial solution */ 2113 if( *result != SCIP_FOUNDSOL ) 2114 { 2115 SCIPdebugMsg(scip, "Solution has value: %.2f\n", SCIPgetSolOrigObj(scip, newsol)); 2116 SCIP_CALL( SCIPtrySolFree(scip, &newsol, TRUE, FALSE, TRUE, TRUE, TRUE, &success) ); 2117 if( success ) 2118 { 2119 SCIPdebugMsg(scip, "Solution copy successful\n"); 2120 *result = SCIP_FOUNDSOL; 2121 } 2122 } 2123 else 2124 { 2125 SCIP_CALL( SCIPfreeSol(scip, &newsol) ); 2126 } 2127 2128 goto TERMINATE; 2129 } 2130 /* current partition is not feasible: 2131 * - update partition 2132 * - update penalty parameters lambda 2133 * - reuse last solution 2134 */ 2135 else 2136 { 2137 int* nviolatedblocksrhs; /* number of blocks which violate rhs for each linking constraint */ 2138 int* nviolatedblockslhs; /* number of blocks which violate lhs for each linking constraint */ 2139 SCIP_Bool update = FALSE; 2140 2141 SCIP_CALL( SCIPallocBufferArray(scip, &nviolatedblocksrhs, heurdata->nlinking) ); 2142 SCIP_CALL( SCIPallocBufferArray(scip, &nviolatedblockslhs, heurdata->nlinking) ); 2143 BMSclearMemoryArray(nviolatedblocksrhs, heurdata->nlinking); 2144 BMSclearMemoryArray(nviolatedblockslhs, heurdata->nlinking); 2145 2146 SCIPdebugMsg(scip, "update partition\n"); 2147 SCIP_CALL( updatePartition(scip, linkings, blockproblem, &nviolatedblocksrhs, &nviolatedblockslhs, 2148 heurdata->nlinking, nblocks, k, &update) ); 2149 2150 /* terminate if partition is not updated */ 2151 if( !update ) 2152 { 2153 SCIPfreeBufferArray(scip, &nviolatedblockslhs); 2154 SCIPfreeBufferArray(scip, &nviolatedblocksrhs); 2155 goto TERMINATE; 2156 } 2157 2158 /* in order to change objective function */ 2159 for( b = 0; b < heurdata->nblocks; b++ ) 2160 { 2161 SCIP_CALL( SCIPfreeTransform(blockproblem[b]->blockscip) ); 2162 } 2163 2164 SCIPdebugMsg(scip, "update lambda\n"); 2165 SCIP_CALL( updateLambda(scip, heurdata, linkings, blockproblem, nviolatedblocksrhs, nviolatedblockslhs, heurdata->nlinking) ); 2166 2167 if( heurdata->reuse ) 2168 { 2169 /* reuse old solution in each block if available */ 2170 SCIP_CALL( reuseSolution(linkings, blockproblem, nblocks) ); 2171 } 2172 2173 SCIPfreeBufferArray(scip, &nviolatedblockslhs); 2174 SCIPfreeBufferArray(scip, &nviolatedblocksrhs); 2175 } 2176 } 2177 SCIPdebugMsg(scip, "maximum number of iterations reached\n"); 2178 2179 /* ------------------------------------------------------------------------ */ 2180 /* free memory */ 2181 TERMINATE: 2182 if( linkings != NULL ) 2183 { 2184 for( c = heurdata->nlinking - 1; c >= 0; c-- ) 2185 { 2186 if( linkings[c]->currentlhs != NULL ) 2187 SCIPfreeBufferArray(scip, &(linkings[c])->currentlhs); 2188 2189 if( linkings[c]->currentrhs != NULL ) 2190 SCIPfreeBufferArray(scip, &(linkings[c])->currentrhs); 2191 } 2192 2193 for( c = heurdata->nlinking - 1; c >= 0; c-- ) 2194 { 2195 linkings[c]->linkingcons = NULL; 2196 SCIPfreeBufferArray(scip, &(linkings[c])->maxactivity); 2197 SCIPfreeBufferArray(scip, &(linkings[c])->minactivity); 2198 SCIPfreeBufferArray(scip, &(linkings[c])->blocknumbers); 2199 SCIPfreeBufferArray(scip, &(linkings[c])->slacks); 2200 SCIPfreeBufferArray(scip, &(linkings[c])->blockconss); 2201 SCIPfreeBlockMemory(scip, &(linkings[c])); /*lint !e866*/ 2202 } 2203 SCIPfreeBufferArray(scip, &linkings); 2204 } 2205 2206 if( blockproblem != NULL ) 2207 { 2208 for( b = nblocks - 1; b >= 0; b-- ) 2209 { 2210 SCIPfreeBufferArray(scip, &(blockproblem[b])->origobj); 2211 SCIPfreeBufferArray(scip, &(blockproblem[b])->slackvars); 2212 SCIPfreeBufferArray(scip, &(blockproblem[b])->linkingindices); 2213 SCIPfreeBufferArray(scip, &(blockproblem[b])->linkingconss); 2214 SCIP_CALL( SCIPfree(&blockproblem[b]->blockscip) ); 2215 SCIPfreeBlockMemory(scip, &(blockproblem[b])); /*lint !e866*/ 2216 } 2217 SCIPfreeBufferArray(scip, &blockproblem); 2218 } 2219 2220 if( assigneddecomp != NULL ) 2221 SCIPdecompFree(&assigneddecomp, SCIPblkmem(scip)); 2222 2223 if( sortedconslabels != NULL ) 2224 SCIPfreeBufferArray(scip, &sortedconslabels); 2225 2226 if( sortedvarlabels != NULL ) 2227 SCIPfreeBufferArray(scip, &sortedvarlabels); 2228 2229 if( sortedconss != NULL ) 2230 SCIPfreeBufferArray(scip, &sortedconss); 2231 2232 if( sortedvars != NULL ) 2233 SCIPfreeBufferArray(scip, &sortedvars); 2234 2235 SCIPdebugMsg(scip, "Leave DPS heuristic\n"); 2236 2237 return SCIP_OKAY; 2238 } 2239 2240 2241 /* 2242 * primal heuristic specific interface methods 2243 */ 2244 2245 /** creates the dps primal heuristic and includes it in SCIP */ 2246 SCIP_RETCODE SCIPincludeHeurDps( 2247 SCIP* scip /**< SCIP data structure */ 2248 ) 2249 { 2250 SCIP_HEURDATA* heurdata; 2251 SCIP_HEUR* heur; 2252 2253 /* create dps primal heuristic data */ 2254 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 2255 2256 heur = NULL; 2257 2258 /* include primal heuristic */ 2259 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 2260 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 2261 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecDps, heurdata) ); 2262 2263 assert(heur != NULL); 2264 2265 /* set non fundamental callbacks via setter functions */ 2266 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyDps) ); 2267 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeDps) ); 2268 2269 /* add dps primal heuristic parameters */ 2270 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiterations", 2271 "maximal number of iterations", &heurdata->maxit, FALSE, DEFAULT_MAXIT, 1, INT_MAX, NULL, NULL) ); 2272 2273 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxlinkscore", 2274 "maximal linking score of used decomposition (equivalent to percentage of linking constraints)", 2275 &heurdata->maxlinkscore, FALSE, 1.0, 0.0, 1.0, NULL, NULL) ); 2276 2277 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/penalty", 2278 "multiplier for absolute increase of penalty parameters (0: no increase)", 2279 &heurdata->penalty, FALSE, DEFAULT_PENALTY, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 2280 2281 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reoptimize", 2282 "should the problem get reoptimized with the original objective function?", &heurdata->reoptimize, FALSE, FALSE, NULL, NULL) ); 2283 2284 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reuse", 2285 "should solutions get reused in subproblems?", &heurdata->reuse, FALSE, FALSE, NULL, NULL) ); 2286 2287 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reoptlimits", 2288 "should strict limits for reoptimization be set?", &heurdata->reoptlimits, FALSE, TRUE, NULL, NULL) ); 2289 2290 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/timing", 2291 "should the heuristic run before or after the processing of the node? (0: before, 1: after, 2: both)", 2292 &heurdata->timing, FALSE, 0, 0, 2, NULL, NULL) ); 2293 2294 return SCIP_OKAY; 2295 } 2296