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_padm.c 26 * @brief PADM primal heuristic based on ideas published in the paper 27 * "A Decomposition Heuristic for Mixed-Integer Supply Chain Problems" 28 * by Martin Schmidt, Lars Schewe, and Dieter Weninger 29 * @author Dieter Weninger 30 * @author Katrin Halbig 31 * 32 * The penalty alternating direction method (PADM) heuristic is a construction heuristic which additionally needs a 33 * user decomposition with linking variables only. 34 * 35 * PADM splits the problem into several sub-SCIPs according to the decomposition, whereby the linking variables get 36 * copied and the difference is penalized. Then the sub-SCIPs are solved on an alternating basis until they arrive at 37 * the same values of the linking variables (ADM-loop). If they don't reconcile after a couple of iterations, 38 * the penalty parameters are increased (penalty-loop) and the sub-SCIPs are solved again on an alternating basis. 39 */ 40 41 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 42 43 #include <assert.h> 44 #include <string.h> 45 46 #include "blockmemshell/memory.h" 47 #include "scip/cons_linear.h" 48 #include "scip/debug.h" 49 #include "scip/heur_padm.h" 50 #include "scip/heuristics.h" 51 #include "scip/pub_cons.h" 52 #include "scip/pub_tree.h" 53 #include "scip/pub_heur.h" 54 #include "scip/pub_message.h" 55 #include "scip/pub_misc.h" 56 #include "scip/pub_misc_select.h" 57 #include "scip/pub_sol.h" 58 #include "scip/pub_var.h" 59 #include "scip/scipdefplugins.h" 60 #include "scip/scip_branch.h" 61 #include "scip/scip_cons.h" 62 #include "scip/scip_copy.h" 63 #include "scip/scip_dcmp.h" 64 #include "scip/scip_general.h" 65 #include "scip/scip_heur.h" 66 #include "scip/scip_lp.h" 67 #include "scip/scip_mem.h" 68 #include "scip/scip_message.h" 69 #include "scip/scip_nodesel.h" 70 #include "scip/scip_numerics.h" 71 #include "scip/scip_param.h" 72 #include "scip/scip_prob.h" 73 #include "scip/scip_randnumgen.h" 74 #include "scip/scip_sol.h" 75 #include "scip/scip_solve.h" 76 #include "scip/scip_solvingstats.h" 77 #include "scip/scip_table.h" 78 #include "scip/scip_timing.h" 79 #include "scip/scip_tree.h" 80 #include "scip/scip_var.h" 81 82 #define HEUR_NAME "padm" 83 #define HEUR_DESC "penalty alternating direction method primal heuristic" 84 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_LNS 85 #define HEUR_PRIORITY 70000 86 #define HEUR_FREQ 0 87 #define HEUR_FREQOFS 0 88 #define HEUR_MAXDEPTH -1 89 #define HEUR_TIMING SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_AFTERNODE 90 #define HEUR_USESSUBSCIP TRUE /**< does the heuristic use a secondary SCIP instance? */ 91 92 #define COUPLINGSIZE 3 93 #define DEFAULT_MINNODES 50LL 94 #define DEFAULT_MAXNODES 5000LL 95 #define DEFAULT_NODEFAC 0.8 96 #define DEFAULT_ADMIT 4 97 #define DEFAULT_PENALTYIT 100 98 #define DEFAULT_GAP 2.0 99 100 /* 101 * Data structures 102 */ 103 104 /** data related to one problem (see below) */ 105 typedef struct Problem PROBLEM; 106 107 /** data related to one block */ 108 typedef struct Block 109 { 110 PROBLEM* problem; /**< the problem this block belongs to */ 111 SCIP* subscip; /**< sub-SCIP representing this block */ 112 int number; /**< component number */ 113 SCIP_VAR** subvars; /**< variables belonging to this block (without slack variables) */ 114 int nsubvars; /**< number of variables belonging to this block (without slack variables) */ 115 SCIP_VAR** slackspos; /**< positive slack variables */ 116 SCIP_VAR** slacksneg; /**< negative slack variables */ 117 SCIP_CONS** couplingcons; /**< coupling contraints */ 118 int ncoupling; /**< number of coupling contraints (equal to positive/negative slack variables) */ 119 SCIP_Real size; /**< share of total problem */ 120 } BLOCK; 121 122 /** data related to one problem */ 123 struct Problem 124 { 125 SCIP* scip; /**< the SCIP instance this problem belongs to */ 126 char* name; /**< name of the problem */ 127 BLOCK* blocks; /**< blocks into which the problem will be divided */ 128 int nblocks; /**< number of blocks */ 129 }; 130 131 /** set data structure */ 132 typedef struct set 133 { 134 int size; /**< size of the set */ 135 int* indexes; /**< set of indexes */ 136 } SET; 137 138 /** data of one linking variable related to one block */ 139 typedef struct blockinfo 140 { 141 int block; /**< index of this block */ 142 int otherblock; /**< index of the other connected block */ 143 int linkvaridx; /**< linking variable index */ 144 SCIP_Real linkvarval; /**< value of linking variable */ 145 SCIP_VAR* linkvar; /**< linking variable */ 146 SCIP_Real slackposobjcoeff; /**< penalty coefficient of positive slack variable */ 147 SCIP_VAR* slackposvar; /**< positive slack variable */ 148 SCIP_Real slacknegobjcoeff; /**< penalty coefficient of negative slack variable */ 149 SCIP_VAR* slacknegvar; /**< negative slack variable */ 150 SCIP_CONS* couplingCons; /**< coupling contraint (equation) */ 151 } BLOCKINFO; 152 153 /** returns TRUE iff both keys are equal */ 154 static 155 SCIP_DECL_HASHKEYEQ(indexesEqual) 156 { /*lint --e{715}*/ 157 BLOCKINFO* binfo1; 158 BLOCKINFO* binfo2; 159 160 binfo1 = (BLOCKINFO*) key1; 161 binfo2 = (BLOCKINFO*) key2; 162 163 if( binfo1->block != binfo2->block || binfo1->otherblock != binfo2->otherblock || 164 binfo1->linkvaridx != binfo2->linkvaridx ) 165 return FALSE; 166 167 return TRUE; 168 } 169 170 /** returns the hash value of the key */ 171 static 172 SCIP_DECL_HASHKEYVAL(indexesHashval) 173 { /*lint --e{715}*/ 174 BLOCKINFO* binfo; 175 binfo = (BLOCKINFO*) key; 176 177 return SCIPhashFour(SCIPrealHashCode((double)binfo->block), SCIPrealHashCode((double)binfo->otherblock), 178 SCIPrealHashCode((double)binfo->linkvaridx), SCIPrealHashCode((double)binfo->linkvaridx)); 179 } 180 181 /** primal heuristic data */ 182 struct SCIP_HeurData 183 { 184 SCIP_Longint maxnodes; /**< maximum number of nodes to regard in all subproblems */ 185 SCIP_Longint minnodes; /**< minimum number of nodes to regard in one subproblem */ 186 int admiterations; /**< maximal number of ADM iterations in each penalty loop */ 187 int penaltyiterations; /**< maximal number of penalty iterations */ 188 int timing; /**< should the heuristic run before or after the processing of the node? 189 (0: before, 1: after, 2: both) */ 190 SCIP_Real nodefac; /**< factor to control nodelimits of subproblems */ 191 SCIP_Real gap; /**< mipgap at start */ 192 SCIP_Bool reoptimize; /**< should the problem get reoptimized with the original objective function? */ 193 SCIP_Bool scaling; /**< enable sigmoid rescaling of penalty parameters */ 194 SCIP_Bool assignlinking; /**< should linking constraints be assigned? */ 195 SCIP_Bool original; /**< should the original problem be used? */ 196 }; 197 198 /* 199 * Local methods 200 */ 201 202 /** initializes one block */ 203 static 204 SCIP_RETCODE initBlock( 205 PROBLEM* problem /**< problem structure */ 206 ) 207 { 208 BLOCK* block; 209 210 assert(problem != NULL); 211 assert(problem->scip != NULL); 212 213 block = &problem->blocks[problem->nblocks]; 214 215 block->problem = problem; 216 block->subscip = NULL; 217 block->subvars = NULL; 218 block->nsubvars = 0; 219 block->number = problem->nblocks; 220 block->slackspos = NULL; 221 block->slacksneg = NULL; 222 block->couplingcons = NULL; 223 block->ncoupling = 0; 224 block->size = 0; 225 226 ++problem->nblocks; 227 228 return SCIP_OKAY; 229 } 230 231 /** frees component structure */ 232 static 233 SCIP_RETCODE freeBlock( 234 BLOCK* block /**< block structure */ 235 ) 236 { 237 assert(block != NULL); 238 239 block->ncoupling = 0; 240 241 if( block->subvars != NULL ) 242 { 243 SCIPfreeBufferArray(block->problem->scip, &(block->subvars)); 244 } 245 246 if( block->subscip != NULL ) 247 { 248 SCIP_CALL( SCIPfree(&block->subscip) ); 249 } 250 251 return SCIP_OKAY; 252 } 253 254 /** initializes subproblem structure */ 255 static 256 SCIP_RETCODE initProblem( 257 SCIP* scip, /**< SCIP data structure */ 258 PROBLEM** problem, /**< pointer to problem structure */ 259 int nblocks /**< number of blocks */ 260 ) 261 { 262 char name[SCIP_MAXSTRLEN]; 263 264 assert(scip != NULL); 265 assert(problem != NULL); 266 267 SCIP_CALL( SCIPallocBlockMemory(scip, problem) ); 268 assert(*problem != NULL); 269 270 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s", SCIPgetProbName(scip)); 271 272 SCIP_CALL( SCIPduplicateMemoryArray(scip, &(*problem)->name, name, strlen(name) + 1) ); 273 274 SCIPdebugMessage("initialized problem <%s>\n", (*problem)->name); 275 276 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*problem)->blocks, nblocks) ); 277 278 (*problem)->scip = scip; 279 (*problem)->nblocks = 0; 280 281 return SCIP_OKAY; 282 } 283 284 /** frees subproblem structure */ 285 static 286 SCIP_RETCODE freeProblem( 287 PROBLEM** problem, /**< pointer to problem to free */ 288 int nblocks /**< number of blocks in decomposition */ 289 ) 290 { 291 SCIP* scip; 292 int c; 293 294 assert(problem != NULL); 295 assert(*problem != NULL); 296 297 scip = (*problem)->scip; 298 assert(scip != NULL); 299 300 /* free all blocks */ 301 for( c = nblocks - 1; c >= 0; --c ) 302 { 303 SCIP_CALL( freeBlock(&(*problem)->blocks[c]) ); 304 } 305 if( (*problem)->blocks != NULL ) 306 { 307 SCIPfreeBlockMemoryArray(scip, &(*problem)->blocks, nblocks); 308 } 309 310 /* free problem name */ 311 SCIPfreeMemoryArray(scip, &(*problem)->name); 312 313 /* free PROBLEM struct and set the pointer to NULL */ 314 SCIPfreeBlockMemory(scip, problem); 315 *problem = NULL; 316 317 return SCIP_OKAY; 318 } 319 320 /** creates a sub-SCIP for the given variables and constraints */ 321 static 322 SCIP_RETCODE createSubscip( 323 SCIP* scip, /**< main SCIP data structure */ 324 SCIP** subscip /**< pointer to store created sub-SCIP */ 325 ) 326 { 327 SCIP_Real infvalue; 328 329 /* create a new SCIP instance */ 330 SCIP_CALL( SCIPcreate(subscip) ); 331 332 SCIP_CALL( SCIPincludeDefaultPlugins(*subscip) ); 333 334 /* copy value for infinity */ 335 SCIP_CALL( SCIPgetRealParam(scip, "numerics/infinity", &infvalue) ); 336 SCIP_CALL( SCIPsetRealParam(*subscip, "numerics/infinity", infvalue) ); 337 338 SCIP_CALL( SCIPcopyLimits(scip, *subscip) ); 339 340 /* avoid recursive calls */ 341 SCIP_CALL( SCIPsetSubscipsOff(*subscip, TRUE) ); 342 343 /* disable cutting plane separation */ 344 SCIP_CALL( SCIPsetSeparating(*subscip, SCIP_PARAMSETTING_OFF, TRUE) ); 345 346 /* disable expensive presolving */ 347 SCIP_CALL( SCIPsetPresolving(*subscip, SCIP_PARAMSETTING_FAST, TRUE) ); 348 349 /* disable expensive techniques */ 350 SCIP_CALL( SCIPsetIntParam(*subscip, "misc/usesymmetry", 0) ); 351 352 /* do not abort subproblem on CTRL-C */ 353 SCIP_CALL( SCIPsetBoolParam(*subscip, "misc/catchctrlc", FALSE) ); 354 355 #ifdef SCIP_DEBUG 356 /* for debugging, enable full output */ 357 SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 5) ); 358 SCIP_CALL( SCIPsetIntParam(*subscip, "display/freq", 100000000) ); 359 #else 360 /* disable statistic timing inside sub SCIP and output to console */ 361 SCIP_CALL( SCIPsetIntParam(*subscip, "display/verblevel", 0) ); 362 SCIP_CALL( SCIPsetBoolParam(*subscip, "timing/statistictiming", FALSE) ); 363 #endif 364 365 return SCIP_OKAY; 366 } 367 368 /** copies the given constraints and the corresponding variables to the given sub-SCIP */ 369 static 370 SCIP_RETCODE copyToSubscip( 371 SCIP* scip, /**< source SCIP */ 372 SCIP* subscip, /**< target SCIP */ 373 const char* name, /**< name for copied problem */ 374 SCIP_CONS** conss, /**< constraints to copy */ 375 SCIP_HASHMAP* varmap, /**< hashmap used for the copy process of variables */ 376 SCIP_HASHMAP* consmap, /**< hashmap used for the copy process of constraints */ 377 int nconss, /**< number of constraints to copy */ 378 SCIP_Bool useorigprob, /**< do we use the original problem? */ 379 SCIP_Bool* success /**< pointer to store whether copying was successful */ 380 ) 381 { 382 SCIP_CONS* newcons; 383 int i; 384 385 assert(scip != NULL); 386 assert(subscip != NULL); 387 assert(conss != NULL); 388 assert(consmap != NULL); 389 assert(success != NULL); 390 391 *success = TRUE; 392 newcons = NULL; 393 394 /* create problem in sub-SCIP */ 395 SCIP_CALL( SCIPcopyProb(scip, subscip, varmap, consmap, FALSE, name) ); 396 397 /* copy constraints */ 398 for( i = 0; i < nconss; ++i ) 399 { 400 assert(conss[i] != NULL); 401 402 /* do not check this if we use the original problem 403 * Since constraints can be deleted etc. during presolving, these assertions would fail. 404 */ 405 if( !useorigprob ) 406 { 407 assert(!SCIPconsIsModifiable(conss[i])); 408 assert(SCIPconsIsActive(conss[i])); 409 assert(!SCIPconsIsDeleted(conss[i])); 410 } 411 412 /* copy the constraint */ 413 SCIP_CALL( SCIPgetConsCopy(scip, subscip, conss[i], &newcons, SCIPconsGetHdlr(conss[i]), varmap, consmap, NULL, 414 SCIPconsIsInitial(conss[i]), SCIPconsIsSeparated(conss[i]), SCIPconsIsEnforced(conss[i]), 415 SCIPconsIsChecked(conss[i]), SCIPconsIsPropagated(conss[i]), FALSE, FALSE, 416 SCIPconsIsDynamic(conss[i]), SCIPconsIsRemovable(conss[i]), FALSE, FALSE, success) ); 417 418 /* abort if constraint was not successfully copied */ 419 if( !(*success) || newcons == NULL) 420 return SCIP_OKAY; 421 422 SCIP_CALL( SCIPaddCons(subscip, newcons) ); 423 SCIP_CALL( SCIPreleaseCons(subscip, &newcons) ); 424 } 425 426 return SCIP_OKAY; 427 } 428 429 /** creates the subscip for a given block */ 430 static 431 SCIP_RETCODE blockCreateSubscip( 432 BLOCK* block, /**< block structure */ 433 SCIP_HASHMAP* varmap, /**< variable hashmap used to improve performance */ 434 SCIP_HASHMAP* consmap, /**< constraint hashmap used to improve performance */ 435 SCIP_CONS** conss, /**< constraints contained in this block */ 436 int nconss, /**< number of constraints contained in this block */ 437 SCIP_Bool useorigprob, /**< do we use the original problem? */ 438 SCIP_Bool* success /**< pointer to store whether the copying process was successful */ 439 ) 440 { 441 char name[SCIP_MAXSTRLEN]; 442 PROBLEM* problem; 443 SCIP* scip; 444 SCIP_VAR** subscipvars; 445 int nsubscipvars; 446 int i; 447 448 assert(block != NULL); 449 assert(varmap != NULL); 450 assert(consmap != NULL); 451 assert(conss != NULL); 452 assert(success != NULL); 453 454 problem = block->problem; 455 assert(problem != NULL); 456 457 scip = problem->scip; 458 assert(scip != NULL); 459 460 (*success) = TRUE; 461 462 SCIP_CALL( createSubscip(scip, &block->subscip) ); 463 464 if( block->subscip != NULL ) 465 { 466 /* get name of the original problem and add "comp_nr" */ 467 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_comp_%d", problem->name, block->number); 468 469 SCIP_CALL( copyToSubscip(scip, block->subscip, name, conss, varmap, consmap, nconss, useorigprob, success) ); 470 471 if( !(*success) ) 472 { 473 SCIP_CALL( SCIPfree(&block->subscip) ); 474 block->subscip = NULL; 475 } 476 else 477 { 478 /* save variables of subscip (without slack variables) */ 479 nsubscipvars = SCIPgetNOrigVars(block->subscip); 480 subscipvars = SCIPgetOrigVars(block->subscip); 481 SCIP_CALL( SCIPallocBufferArray(scip, &(block->subvars), nsubscipvars) ); 482 block->nsubvars = nsubscipvars; 483 for( i = 0; i < nsubscipvars; i++ ) 484 block->subvars[i] = subscipvars[i]; 485 486 /* calculate size of sub-SCIP with focus on the number of integer variables 487 * we use this value to determine the nodelimit 488 */ 489 block->size = (SCIP_Real)(SCIPgetNOrigVars(block->subscip) + SCIPgetNOrigIntVars(block->subscip) + SCIPgetNOrigBinVars(block->subscip)) / 490 (SCIPgetNVars(scip) + SCIPgetNIntVars(scip) + SCIPgetNBinVars(scip)); 491 } 492 } 493 else 494 (*success) = FALSE; 495 496 SCIPdebugMsg(scip, "created subscip of block %d\n", block->number); 497 498 return SCIP_OKAY; 499 } 500 501 /** creates problem structure and split it into blocks */ 502 static 503 SCIP_RETCODE createAndSplitProblem( 504 SCIP* scip, /**< SCIP data structure */ 505 SCIP_CONS** sortedconss, /**< array of (checked) constraints sorted by blocks */ 506 int nconss, /**< number of constraints */ 507 int* consssize, /**< number of constraints per block (and border at index 0) */ 508 int nblocks, /**< number of blocks */ 509 PROBLEM** problem, /**< pointer to store problem structure */ 510 SCIP_Bool useorigprob, /**< do we use the original problem? */ 511 SCIP_Bool* success /**< pointer to store whether the process was successful */ 512 ) 513 { 514 BLOCK* block; 515 SCIP_HASHMAP* varmap; 516 SCIP_HASHMAP* consmap; 517 SCIP_CONS** blockconss; 518 int nhandledconss; 519 int nblockconss; 520 int b; 521 522 (*success) = TRUE; 523 524 /* init subproblem data structure */ 525 SCIP_CALL( initProblem(scip, problem, nblocks) ); 526 assert((*problem)->blocks != NULL); 527 528 /* hashmap mapping from original constraints to constraints in the sub-SCIPs (for performance reasons) */ 529 SCIP_CALL( SCIPhashmapCreate(&consmap, SCIPblkmem(scip), nconss) ); 530 531 for( b = 0; b < nblocks; b++ ) 532 { 533 SCIP_CALL( initBlock(*problem) ); 534 } 535 536 /* loop over all blocks and create subscips */ 537 nhandledconss = 0; 538 for( b = 0; b < nblocks; b++ ) 539 { 540 block = &(*problem)->blocks[b]; 541 542 SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 543 544 /* get block constraints */ 545 blockconss = &(sortedconss[nhandledconss]); 546 nblockconss = consssize[b + 1]; 547 548 /* build subscip for block */ 549 SCIP_CALL( blockCreateSubscip(block, varmap, consmap, blockconss, nblockconss, useorigprob, success) ); 550 551 SCIPhashmapFree(&varmap); 552 nhandledconss += nblockconss; 553 554 if( !(*success) ) 555 break; 556 } 557 558 SCIPhashmapFree(&consmap); 559 560 if( !(*success) ) 561 { 562 /* free subproblem data structure since not all blocks could be copied */ 563 SCIP_CALL( freeProblem(problem, nblocks) ); 564 } 565 566 return SCIP_OKAY; 567 } 568 569 /** copies labels to newdecomp and assigns linking constraints if possible*/ 570 static 571 SCIP_RETCODE assignLinking( 572 SCIP* scip, /**< SCIP data structure */ 573 SCIP_DECOMP* newdecomp, /**< decomposition with (partially) assigned linking constraints */ 574 SCIP_VAR** vars, /**< array of variables */ 575 SCIP_CONS** sortedconss, /**< sorted array of constraints */ 576 int* varlabels, /**< array of variable labels */ 577 int* conslabels, /**< sorted array of constraint labels */ 578 int nvars, /**< number of variables */ 579 int nconss, /**< number of constraints */ 580 int nlinkconss /**< number of linking constraints */ 581 ) 582 { 583 assert(scip != NULL); 584 assert(vars != NULL); 585 assert(sortedconss != NULL); 586 assert(varlabels != NULL); 587 assert(conslabels != NULL); 588 589 /* copy the labels */ 590 SCIP_CALL( SCIPdecompSetVarsLabels(newdecomp, vars, varlabels, nvars) ); 591 SCIP_CALL( SCIPdecompSetConsLabels(newdecomp, sortedconss, conslabels, nconss) ); 592 593 SCIPdebugMsg(scip, "try to assign %d linking constraints\n", nlinkconss); 594 595 /* reassign linking constraints */ 596 SCIP_CALL( SCIPassignDecompLinkConss(scip, newdecomp, &sortedconss[0], nlinkconss, NULL) ); 597 598 SCIP_CALL( SCIPcomputeDecompVarsLabels(scip, newdecomp, sortedconss, nconss) ); 599 600 SCIP_CALL( SCIPcomputeDecompStats(scip, newdecomp, TRUE) ); 601 602 SCIPdecompGetConsLabels(newdecomp, sortedconss, conslabels, nconss); 603 SCIPdecompGetVarsLabels(newdecomp, vars, varlabels, nvars); 604 605 SCIPsortIntPtr(conslabels, (void**)sortedconss, nconss); 606 607 return SCIP_OKAY; 608 } 609 610 /** computes feasible solution from last stored solution of the block*/ 611 static 612 SCIP_RETCODE reuseSolution( 613 SCIP* subscip, /**< SCIP data structure */ 614 BLOCK* block /**< block structure*/ 615 ) 616 { 617 SCIP_SOL** sols; 618 SCIP_SOL* sol; /* solution of block that will be repaired */ 619 SCIP_SOL* newsol; 620 SCIP_VAR** blockvars; 621 SCIP_VAR** consvars; 622 SCIP_Real* blockvals; 623 int nsols; 624 int nvars; 625 int c; 626 SCIP_Bool success; 627 628 assert(subscip != NULL); 629 assert(block != NULL); 630 631 nsols = SCIPgetNSols(subscip); 632 633 /* no solution in solution candidate storage found */ 634 if( nsols == 0 ) 635 return SCIP_OKAY; 636 637 SCIP_CALL( SCIPallocBufferArray(subscip, &consvars, COUPLINGSIZE) ); 638 639 sols = SCIPgetSols(subscip); 640 sol = sols[nsols - 1]; 641 642 /* copy the solution */ 643 nvars = SCIPgetNVars(subscip); 644 blockvars = SCIPgetVars(subscip); 645 SCIP_CALL( SCIPallocBufferArray(subscip, &blockvals, nvars) ); 646 SCIP_CALL( SCIPgetSolVals(subscip, sol, nvars, blockvars, blockvals) ); 647 SCIP_CALL( SCIPcreateOrigSol(subscip, &newsol, NULL) ); 648 SCIP_CALL( SCIPsetSolVals(subscip, newsol, nvars, blockvars, blockvals) ); 649 650 /* correct each coupling constraint; 651 * orig_var + slackpos - slackneg == side 652 * adapt slack variables so that constraint is feasible */ 653 for( c = 0; c < block->ncoupling; c++ ) 654 { 655 SCIP_Real solval; /* old solution values of variables; [0] original variable, [1] slackpos, [2] slackneg */ 656 SCIP_Real side; /* current right hand side */ 657 SCIP_Real diff; 658 659 SCIP_CALL( SCIPgetConsVars(subscip, block->couplingcons[c], consvars, COUPLINGSIZE, &success) ); 660 solval = SCIPgetSolVal(subscip, sol, consvars[0]); 661 662 side = SCIPgetRhsLinear(subscip, block->couplingcons[c]); 663 assert(SCIPisEQ(subscip, SCIPgetRhsLinear(subscip, block->couplingcons[c]), SCIPgetLhsLinear(subscip, block->couplingcons[c]))); 664 665 diff = side - solval; 666 667 /* slackpos is strict positiv */ 668 if( diff > 0 ) 669 { 670 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slackspos[c], diff) ); 671 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slacksneg[c], 0.0) ); 672 } 673 /* slackneg is strict positiv */ 674 else if( diff < 0 ) 675 { 676 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slacksneg[c], -diff) ); 677 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slackspos[c], 0.0) ); 678 } 679 /* no slack variable necessary */ 680 else 681 { 682 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slackspos[c], 0.0) ); 683 SCIP_CALL( SCIPsetSolVal(subscip, newsol, block->slacksneg[c], 0.0) ); 684 } 685 } 686 687 SCIPdebugMsg(subscip, "Try adding solution with objective value %.2f\n", SCIPgetSolOrigObj(subscip, newsol)); 688 SCIP_CALL( SCIPaddSolFree(subscip, &newsol, &success) ); 689 690 if( !success ) 691 SCIPdebugMsg(subscip, "Correcting solution failed\n"); /* maybe not better than old solutions */ 692 else 693 { 694 SCIPdebugMsg(subscip, "Correcting solution successful\n"); 695 } 696 697 SCIPfreeBufferArray(subscip, &blockvals); 698 SCIPfreeBufferArray(subscip, &consvars); 699 700 return SCIP_OKAY; 701 } 702 703 /** reoptimizes the heuristic solution with original objective function 704 * 705 * Since the main algorithm of padm ignores the objective function, this method can be called to obtain better solutions. 706 * It copies the main scip, fixes the linking variables at the values of the already found solution 707 * and solves the new problem with small limits. 708 */ 709 static 710 SCIP_RETCODE reoptimize( 711 SCIP* scip, /**< SCIP data structure */ 712 SCIP_HEUR* heur, /**< pointer to heuristic*/ 713 SCIP_SOL* sol, /**< heuristic solution */ 714 SCIP_VAR** vars, /**< pointer to variables */ 715 int nvars, /**< number of variables */ 716 SCIP_VAR** linkvars, /**< pointer to linking variables */ 717 int nlinkvars, /**< number of linking variables */ 718 SCIP_SOL** newsol, /**< pointer to store improved solution */ 719 SCIP_Bool* success /**< pointer to store whether reoptimization was successful */ 720 ) 721 { 722 SCIP* scipcopy; 723 SCIP_SOL* startsol; 724 SCIP_HASHMAP* varmap; 725 SCIP_VAR** subvars; 726 SCIP_STATUS status; 727 SCIP_Real* linkvals; 728 SCIP_Real time; 729 int v; 730 731 assert(scip != NULL); 732 assert(heur != NULL); 733 assert(sol != NULL); 734 assert(vars != NULL); 735 assert(linkvars != NULL); 736 737 SCIP_CALL( SCIPallocBufferArray(scip, &linkvals, nlinkvars) ); 738 SCIP_CALL( SCIPallocBufferArray(scip, &subvars, nvars) ); 739 740 SCIP_CALL( SCIPgetSolVals(scip, sol, nlinkvars, linkvars, linkvals) ); 741 742 /* initializing the problem copy*/ 743 SCIP_CALL( SCIPcreate(&scipcopy) ); 744 745 /* - create the variable mapping hash map 746 * - create a problem copy of main SCIP 747 */ 748 if( SCIPheurGetData(heur)->original ) 749 { 750 SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(scipcopy), SCIPgetNOrigVars(scip)) ); 751 SCIP_CALL( SCIPcopyOrigConsCompression(scip, scipcopy, varmap, NULL, "reopt_padm", linkvars, linkvals, nlinkvars, 752 FALSE, FALSE, TRUE, success) ); 753 } 754 else 755 { 756 SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(scipcopy), SCIPgetNVars(scip)) ); 757 SCIP_CALL( SCIPcopyConsCompression(scip, scipcopy, varmap, NULL, "reopt_padm", linkvars, linkvals, nlinkvars, 758 TRUE, FALSE, FALSE, TRUE, success) ); 759 } 760 for( v = 0; v < nvars; v++ ) 761 { 762 subvars[v] = (SCIP_VAR*) SCIPhashmapGetImage(varmap, vars[v]); 763 } 764 765 /* do not abort subproblem on CTRL-C */ 766 SCIP_CALL( SCIPsetBoolParam(scipcopy, "misc/catchctrlc", FALSE) ); 767 768 /* forbid recursive call of heuristics and separators solving sub-SCIPs */ 769 SCIP_CALL( SCIPsetSubscipsOff(scipcopy, TRUE) ); 770 771 #ifdef SCIP_DEBUG 772 /* for debugging, enable full output */ 773 SCIP_CALL( SCIPsetIntParam(scipcopy, "display/verblevel", 5) ); 774 SCIP_CALL( SCIPsetIntParam(scipcopy, "display/freq", 100000000) ); 775 #else 776 /* disable statistic timing inside sub SCIP and output to console */ 777 SCIP_CALL( SCIPsetIntParam(scipcopy, "display/verblevel", 0) ); 778 SCIP_CALL( SCIPsetBoolParam(scipcopy, "timing/statistictiming", FALSE) ); 779 #endif 780 781 /* disable cutting plane separation */ 782 SCIP_CALL( SCIPsetSeparating(scipcopy, SCIP_PARAMSETTING_OFF, TRUE) ); 783 784 /* disable expensive presolving but enable components presolver */ 785 SCIP_CALL( SCIPsetPresolving(scipcopy, SCIP_PARAMSETTING_FAST, TRUE) ); 786 SCIP_CALL( SCIPsetIntParam(scipcopy, "constraints/components/maxprerounds", 1) ); 787 SCIP_CALL( SCIPsetLongintParam(scipcopy, "constraints/components/nodelimit", 0LL) ); 788 789 /* disable expensive techniques */ 790 SCIP_CALL( SCIPsetIntParam(scipcopy, "misc/usesymmetry", 0) ); 791 792 /* speed up sub-SCIP by not checking dual LP feasibility */ 793 SCIP_CALL( SCIPsetBoolParam(scipcopy, "lp/checkdualfeas", FALSE) ); 794 795 /* add heuristic solution as start solution */ 796 SCIP_CALL( SCIPtransformProb(scipcopy) ); 797 *success = FALSE; 798 if( SCIPisLT(scip, SCIPgetSolOrigObj(scip, sol), SCIPinfinity(scip)) ) 799 { 800 SCIP_CALL( SCIPcreateSol(scipcopy, &startsol, heur) ); 801 for( v = 0; v < nvars; v++ ) 802 { 803 SCIP_VAR* subvar; 804 subvar = (SCIP_VAR*) SCIPhashmapGetImage(varmap, vars[v]); 805 if( subvar != NULL ) 806 { 807 SCIP_CALL( SCIPsetSolVal(scipcopy, startsol, subvar, SCIPgetSolVal(scip, sol, vars[v])) ); 808 } 809 } 810 811 SCIP_CALL( SCIPtrySolFree(scipcopy, &startsol, FALSE, FALSE, FALSE, FALSE, FALSE, success) ); 812 if( *success ) 813 SCIPdebugMsg(scip, "set start solution\n"); 814 else 815 SCIPdebugMsg(scip, "start solution for reoptimizing is not feasible\n"); 816 } 817 818 /* set limits; do not use more time than the heuristic has already used */ 819 SCIP_CALL( SCIPcopyLimits(scip, scipcopy) ); 820 SCIP_CALL( SCIPsetLongintParam(scipcopy, "limits/nodes", 1LL) ); 821 SCIP_CALL( SCIPgetRealParam(scipcopy, "limits/time", &time) ); 822 if( SCIPheurGetTime(heur) < time - 1.0 ) 823 { 824 SCIP_CALL( SCIPsetRealParam(scipcopy, "limits/time", SCIPheurGetTime(heur) + 1.0) ); 825 } 826 if( *success ) 827 { 828 SCIP_CALL( SCIPsetIntParam(scipcopy, "limits/bestsol", 2) ); /* first solution is start solution */ 829 } 830 else 831 { 832 SCIP_CALL( SCIPsetIntParam(scipcopy, "limits/bestsol", 1) ); 833 } 834 835 /* reoptimize problem */ 836 SCIP_CALL_ABORT( SCIPsolve(scipcopy) ); 837 status = SCIPgetStatus(scipcopy); 838 839 /* copy solution to main scip */ 840 if( status == SCIP_STATUS_BESTSOLLIMIT || status == SCIP_STATUS_OPTIMAL ) 841 { 842 SCIP_SOL* solcopy; 843 844 solcopy = SCIPgetBestSol(scipcopy); 845 SCIP_CALL( SCIPtranslateSubSol(scip, scipcopy, solcopy, heur, subvars, newsol) ); 846 847 SCIPdebugMsg(scip, "Objective value of reoptimized solution %.2f\n", SCIPgetSolOrigObj(scip, *newsol)); 848 *success = TRUE; 849 } 850 else 851 { 852 *success = FALSE; 853 } 854 855 /* free data */ 856 SCIPhashmapFree(&varmap); 857 SCIP_CALL( SCIPfree(&scipcopy) ); 858 SCIPfreeBufferArray(scip, &subvars); 859 SCIPfreeBufferArray(scip, &linkvals); 860 861 return SCIP_OKAY; 862 } 863 864 /** rescales the penalty parameters 865 * 866 * A sigmoid function is a function with an "S"-shaped graph, e.g. S(x) = x/(1+|x|). 867 * In order to avoid numerical instabilities due to large penalty parameters we rescale them 868 * using the sigmoid function 869 * S(x) = (x - shift)/(flatness + |x - shift|) * (range/2) + offset. 870 * The parameters are mapped into the more controllable interval [lowestslack, range + lowestslack]. 871 */ 872 static 873 SCIP_RETCODE scalePenalties( 874 PROBLEM* problem, /**< block structure */ 875 SET* linkvartoblocks, /**< linking variable to blocks set */ 876 SET* blocktolinkvars, /**< block to linking variable set */ 877 SCIP_HASHTABLE* htable, /**< hashtable containing blockinfo*/ 878 SCIP_Real maxpenalty /**< maximum penalty parameter */ 879 ) 880 { 881 SCIP_Real shift; 882 SCIP_Real lowestslack; 883 SCIP_Real range; 884 SCIP_Real offset; 885 SCIP_Real flatness; 886 int b; 887 int i; 888 int k; 889 890 shift = maxpenalty / 2.0; 891 lowestslack = 0.1; 892 range = 10.0; 893 offset = range / 2.0 + lowestslack; 894 flatness = maxpenalty / 10.0; 895 896 for( b = 0; b < problem->nblocks; b++ ) 897 { 898 for( i = 0; i < blocktolinkvars[b].size; i++ ) 899 { 900 int linkvaridx; 901 linkvaridx = blocktolinkvars[b].indexes[i]; 902 903 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 904 { 905 int b2; 906 b2 = linkvartoblocks[linkvaridx].indexes[k]; 907 908 if( b2 != b ) 909 { 910 BLOCKINFO binfo; 911 BLOCKINFO* binfoout; 912 SCIP_Real oldcoeff; 913 914 binfo.block = b; 915 binfo.otherblock = b2; 916 binfo.linkvaridx = linkvaridx; 917 binfoout = (BLOCKINFO*) SCIPhashtableRetrieve(htable, (void*) &binfo); 918 assert(binfoout != NULL); 919 920 /* scale coefficient of positive slack variable */ 921 oldcoeff = binfoout->slackposobjcoeff; 922 binfoout->slackposobjcoeff = ((oldcoeff - shift) / (flatness + REALABS(oldcoeff - shift))) * range / 2.0 + offset; 923 924 /* scale coefficient of negative slack variable */ 925 oldcoeff = binfoout->slacknegobjcoeff; 926 binfoout->slacknegobjcoeff = ((oldcoeff - shift) / (flatness + REALABS(oldcoeff - shift))) * range / 2.0 + offset; 927 } 928 } 929 } 930 } 931 return SCIP_OKAY; 932 } 933 934 /** returns the available time limit that is left */ 935 static 936 SCIP_RETCODE getTimeLeft( 937 SCIP* scip, /**< SCIP data structure */ 938 SCIP_Real* time /**< pointer to store remaining time */ 939 ) 940 { 941 SCIP_Real timelim; 942 SCIP_Real solvingtime; 943 944 assert(scip != NULL); 945 946 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelim) ); 947 solvingtime = SCIPgetSolvingTime(scip); 948 949 if( !SCIPisInfinity(scip, timelim) ) 950 *time = MAX(0.0, (timelim - solvingtime)); 951 else 952 *time = SCIPinfinity(scip); 953 954 return SCIP_OKAY; 955 } 956 957 /* 958 * Callback methods of primal heuristic 959 */ 960 961 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 962 static 963 SCIP_DECL_HEURCOPY(heurCopyPADM) 964 { /*lint --e{715}*/ 965 assert(scip != NULL); 966 assert(heur != NULL); 967 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 968 969 /* call inclusion method of primal heuristic */ 970 SCIP_CALL( SCIPincludeHeurPADM(scip) ); 971 972 return SCIP_OKAY; 973 } 974 975 976 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 977 static 978 SCIP_DECL_HEURFREE(heurFreePADM) 979 { /*lint --e{715}*/ 980 SCIP_HEURDATA* heurdata; 981 982 heurdata = SCIPheurGetData(heur); 983 SCIPheurSetData(heur, NULL); 984 985 assert(heurdata != NULL); 986 987 SCIPfreeBlockMemory(scip, &heurdata); 988 989 return SCIP_OKAY; 990 } 991 992 993 /** execution method of primal heuristic */ 994 static SCIP_DECL_HEUREXEC(heurExecPADM) 995 { /*lint --e{715}*/ 996 char name[SCIP_MAXSTRLEN]; 997 char info[SCIP_MAXSTRLEN]; 998 SCIP_HEURDATA* heurdata; 999 PROBLEM* problem; 1000 SCIP_DECOMP** alldecomps; 1001 SCIP_DECOMP* decomp; 1002 SCIP_DECOMP* assigneddecomp; 1003 SCIP_VAR** vars; 1004 SCIP_VAR** linkvars; 1005 SCIP_VAR** tmpcouplingvars; 1006 SCIP_CONS** conss; 1007 SCIP_CONS** sortedconss; 1008 SET* linkvartoblocks; 1009 SET* blocktolinkvars; 1010 BLOCKINFO* blockinfolist; 1011 SCIP_HASHTABLE* htable; 1012 int* varlabels; 1013 int* conslabels; 1014 int* consssize; 1015 int* alllinkvartoblocks; /* for efficient memory allocation */ 1016 SCIP_Bool* varonlyobj; 1017 SCIP_Real* tmpcouplingcoef; 1018 SCIP_Real gap; 1019 SCIP_Real maxpenalty; 1020 SCIP_Real slackthreshold; 1021 SCIP_Real memory; /* in MB */ 1022 SCIP_Real timeleft; 1023 SCIP_STATUS status; 1024 SCIP_Bool solutionsdiffer; 1025 SCIP_Bool solved; 1026 SCIP_Bool doscaling; 1027 SCIP_Bool istimeleft; 1028 SCIP_Bool success; 1029 SCIP_Bool avoidmemout; 1030 SCIP_Bool disablemeasures; 1031 int maxgraphedge; 1032 int ndecomps; 1033 int nconss; 1034 int nvars; 1035 int nblocks; 1036 int numlinkvars; 1037 int nentries; 1038 int aiter; 1039 int piter; 1040 int increasedslacks; 1041 int blockinfolistfill; 1042 SCIP_Longint nodesleft; 1043 int i; 1044 int b; 1045 int k; 1046 int j; 1047 1048 assert(scip != NULL); 1049 assert(heur != NULL); 1050 assert(result != NULL); 1051 1052 heurdata = SCIPheurGetData(heur); 1053 assert(heurdata != NULL); 1054 1055 *result = SCIP_DIDNOTRUN; 1056 1057 problem = NULL; 1058 assigneddecomp = NULL; 1059 sortedconss = NULL; 1060 varlabels = NULL; 1061 conslabels = NULL; 1062 consssize = NULL; 1063 alllinkvartoblocks = NULL; 1064 linkvars = NULL; 1065 linkvartoblocks = NULL; 1066 blocktolinkvars = NULL; 1067 tmpcouplingvars = NULL; 1068 tmpcouplingcoef = NULL; 1069 varonlyobj = NULL; 1070 blockinfolist = NULL; 1071 htable = NULL; 1072 1073 nodesleft = heurdata->maxnodes; 1074 gap = heurdata->gap; 1075 1076 if( (heurtiming & SCIP_HEURTIMING_BEFORENODE) && heurdata->timing !=1 ) 1077 { 1078 SCIPdebugMsg(scip, "Initialize padm heuristic before node\n"); 1079 } 1080 else if( (heurtiming & SCIP_HEURTIMING_AFTERNODE) && heurdata->timing >=1 ) 1081 { 1082 SCIPdebugMsg(scip, "Initialize padm heuristic after node\n"); 1083 } 1084 else 1085 { 1086 return SCIP_OKAY; 1087 } 1088 1089 #ifdef PADM_WRITE_PROBLEMS 1090 SCIP_CALL( SCIPwriteOrigProblem(scip, "orig_problem", NULL, FALSE) ); 1091 SCIP_CALL( SCIPwriteTransProblem(scip, "trans_problem", NULL, FALSE) ); 1092 #endif 1093 1094 /* take original problem (This is only for testing and not recommended!) */ 1095 if( heurdata->original ) 1096 { 1097 /* multiaggregation of variables has to be switched off */ 1098 if( !SCIPdoNotMultaggr(scip) ) 1099 { 1100 SCIPwarningMessage(scip, "Heuristic %s does not support multiaggregation when the original problem is used.\nPlease turn multiaggregation off to use this feature.\n", HEUR_NAME); 1101 return SCIP_OKAY; 1102 } 1103 1104 SCIPgetDecomps(scip, &alldecomps, &ndecomps, TRUE); 1105 if( ndecomps == 0) 1106 return SCIP_OKAY; 1107 1108 /* it takes the first decomposition */ 1109 decomp = alldecomps[0]; 1110 SCIPdebugMsg(scip, "First original decomposition is selected\n"); 1111 assert(decomp != NULL); 1112 1113 nconss = SCIPgetNOrigConss(scip); 1114 conss = SCIPgetOrigConss(scip); 1115 nvars = SCIPgetNOrigVars(scip); 1116 vars = SCIPgetOrigVars(scip); 1117 } 1118 /* take transformed problem */ 1119 else 1120 { 1121 SCIPgetDecomps(scip, &alldecomps, &ndecomps, FALSE); 1122 if( ndecomps == 0) 1123 return SCIP_OKAY; 1124 1125 /* it takes the first decomposition */ 1126 decomp = alldecomps[0]; 1127 SCIPdebugMsg(scip, "First transformed decomposition is selected\n"); 1128 assert(decomp != NULL); 1129 1130 nconss = SCIPgetNConss(scip); 1131 conss = SCIPgetConss(scip); 1132 nvars = SCIPgetNVars(scip); 1133 vars = SCIPgetVars(scip); 1134 } 1135 1136 nblocks = SCIPdecompGetNBlocks(decomp); 1137 1138 /* if problem has no constraints, no variables or less than two blocks, return */ 1139 if( nconss == 0 || nvars == 0 || nblocks <= 1 ) 1140 { 1141 SCIPdebugMsg(scip, "problem has no constraints, no variables or less than two blocks\n"); 1142 goto TERMINATE; 1143 } 1144 1145 /* estimate required memory for all blocks and terminate if not enough memory is available */ 1146 SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &memory) ); 1147 SCIP_CALL( SCIPgetBoolParam(scip, "misc/avoidmemout", &avoidmemout) ); 1148 if( avoidmemout && (((SCIPgetMemUsed(scip) + SCIPgetMemExternEstim(scip))/1048576.0) * nblocks >= memory) ) 1149 { 1150 SCIPdebugMsg(scip, "The estimated memory usage for %d blocks is too large.\n", nblocks); 1151 goto TERMINATE; 1152 } 1153 1154 /* we do not need the block decomposition graph and expensive measures of the decomposition statistics */ 1155 SCIP_CALL( SCIPgetIntParam(scip, "decomposition/maxgraphedge", &maxgraphedge) ); 1156 if( !SCIPisParamFixed(scip, "decomposition/maxgraphedge") ) 1157 { 1158 SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", 0) ); 1159 } 1160 SCIP_CALL( SCIPgetBoolParam(scip, "decomposition/disablemeasures", &disablemeasures) ); 1161 if( !SCIPisParamFixed(scip, "decomposition/disablemeasures") ) 1162 { 1163 SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", TRUE) ); 1164 } 1165 1166 /* don't change problem by sorting constraints */ 1167 SCIP_CALL( SCIPduplicateBufferArray(scip, &sortedconss, conss, nconss) ); 1168 1169 SCIP_CALL( SCIPallocBufferArray(scip, &varlabels, nvars) ); 1170 SCIP_CALL( SCIPallocBufferArray(scip, &conslabels, nconss) ); 1171 SCIP_CALL( SCIPallocBufferArray(scip, &consssize, nblocks + 1) ); 1172 1173 SCIPdecompGetConsLabels(decomp, conss, conslabels, nconss); 1174 SCIPdecompGetVarsLabels(decomp, vars, varlabels, nvars); 1175 1176 /* sort constraints by blocks */ 1177 SCIPsortIntPtr(conslabels, (void**)sortedconss, nconss); 1178 1179 /* try to assign linking constraints */ 1180 if( heurdata->assignlinking && conslabels[0] == SCIP_DECOMP_LINKCONS ) 1181 { 1182 /* create new decomposition; don't change the decompositions in the decompstore */ 1183 SCIP_CALL( SCIPcreateDecomp(scip, &assigneddecomp, nblocks, heurdata->original, SCIPdecompUseBendersLabels(decomp)) ); 1184 1185 SCIP_CALL( assignLinking(scip, assigneddecomp, vars, sortedconss, varlabels, conslabels, nvars, nconss, SCIPdecompGetNBorderConss(decomp)) ); 1186 assert(SCIPdecompGetNBlocks(decomp) >= SCIPdecompGetNBlocks(assigneddecomp)); 1187 decomp = assigneddecomp; 1188 1189 /* number of blocks can get smaller (since assigning constraints can lead to empty blocks) */ 1190 nblocks = SCIPdecompGetNBlocks(decomp); 1191 } 1192 else 1193 { 1194 /* The decomposition statistics were computed during transformation of the decomposition store. 1195 * Since propagators can have changed the number of constraints/variables, 1196 * the statistics are no longer up-to-date and have to be recomputed. 1197 */ 1198 SCIP_CALL( SCIPcomputeDecompStats(scip, decomp, TRUE) ); 1199 nblocks = SCIPdecompGetNBlocks(decomp); 1200 } 1201 1202 /* reset parameters */ 1203 SCIP_CALL( SCIPsetIntParam(scip, "decomposition/maxgraphedge", maxgraphedge) ); 1204 SCIP_CALL( SCIPsetBoolParam(scip, "decomposition/disablemeasures", disablemeasures) ); 1205 1206 /* @note the terms 'linking' and 'border' (constraints/variables) are used interchangeably */ 1207 1208 if( SCIPdecompGetNBorderConss(decomp) != 0 ) 1209 { 1210 SCIPdebugMsg(scip, "No support for linking contraints\n"); 1211 goto TERMINATE; 1212 } 1213 1214 /* get number of linking variables */ 1215 numlinkvars = SCIPdecompGetNBorderVars(decomp); 1216 SCIPdebugMsg(scip, "%d linking variables\n", numlinkvars); 1217 1218 if( numlinkvars == 0 ) 1219 { 1220 SCIPdebugMsg(scip, "No linking variables\n"); 1221 goto TERMINATE; 1222 } 1223 1224 *result = SCIP_DIDNOTFIND; 1225 1226 /* get for every block the number of constraints (first entry belongs to border) */ 1227 SCIP_CALL( SCIPdecompGetConssSize(decomp, consssize, nblocks + 1) ); 1228 1229 /* create blockproblems */ 1230 SCIP_CALL( createAndSplitProblem(scip, sortedconss, nconss, consssize, nblocks, &problem, heurdata->original, &success) ); 1231 1232 if( !success ) 1233 { 1234 SCIPdebugMsg(scip, "Some subscips could not be created successfully.\n"); 1235 goto TERMINATE; 1236 } 1237 1238 SCIP_CALL( SCIPallocBufferArray(scip, &linkvartoblocks, numlinkvars) ); 1239 SCIP_CALL( SCIPallocBufferArray(scip, &blocktolinkvars, problem->nblocks) ); 1240 1241 /* set pointer to NULL for safe memory release */ 1242 for( i = 0; i < numlinkvars; i++ ) 1243 linkvartoblocks[i].indexes = NULL; 1244 for( i = 0; i < problem->nblocks; i++ ) 1245 blocktolinkvars[i].indexes = NULL; 1246 1247 /* extract linking variables and init linking variable to blocks set */ 1248 SCIP_CALL( SCIPallocBufferArray(scip, &alllinkvartoblocks, problem->nblocks * numlinkvars) ); 1249 SCIP_CALL( SCIPallocBufferArray(scip, &linkvars, numlinkvars) ); 1250 1251 b = 0; 1252 for( i = 0; i < nvars; i++ ) 1253 { 1254 if( varlabels[i] == SCIP_DECOMP_LINKVAR ) 1255 { 1256 linkvars[b] = vars[i]; 1257 linkvartoblocks[b].indexes = &alllinkvartoblocks[b * problem->nblocks]; 1258 linkvartoblocks[b].size = 0; 1259 b++; 1260 } 1261 } 1262 1263 /* fill linking variable to blocks set */ 1264 for( i = 0; i < numlinkvars; i++ ) 1265 { 1266 SCIP_VAR* var; 1267 const char* vname; 1268 1269 vname = SCIPvarGetName(linkvars[i]); 1270 k = 0; 1271 for( b = 0; b < problem->nblocks; b++ ) 1272 { 1273 var = SCIPfindVar((problem->blocks[b]).subscip, vname); 1274 if( var != NULL ) 1275 { 1276 linkvartoblocks[i].indexes[k] = b; 1277 linkvartoblocks[i].size = k + 1; 1278 k++; 1279 } 1280 } 1281 } 1282 1283 /* check whether there is enough time left */ 1284 SCIP_CALL( getTimeLeft(scip, &timeleft) ); 1285 if( timeleft <= 0 ) 1286 { 1287 SCIPdebugMsg(scip, "no time left\n"); 1288 goto TERMINATE; 1289 } 1290 1291 /* init varonlyobj; true if variable is only part of the objective function */ 1292 SCIP_CALL( SCIPallocBufferArray(scip, &varonlyobj, numlinkvars) ); 1293 for( i = 0; i < numlinkvars; ++i) 1294 varonlyobj[i] = TRUE; 1295 1296 /* init and fill block to linking variables set */ 1297 for( b = 0; b < problem->nblocks; b++ ) 1298 { 1299 SCIP_CALL( SCIPallocBufferArray(scip, &(blocktolinkvars[b].indexes), numlinkvars) ); 1300 blocktolinkvars[b].size = 0; 1301 1302 k = 0; 1303 for( i = 0; i < numlinkvars; i++ ) 1304 { 1305 SCIP_VAR* var; 1306 const char* vname; 1307 1308 vname = SCIPvarGetName(linkvars[i]); 1309 var = SCIPfindVar((problem->blocks[b]).subscip, vname); 1310 if( var != NULL ) 1311 { 1312 varonlyobj[i] = FALSE; 1313 blocktolinkvars[b].indexes[k] = i; 1314 blocktolinkvars[b].size = k + 1; 1315 k++; 1316 } 1317 } 1318 } 1319 1320 /* init arrays for slack variables and coupling constraints */ 1321 for( b = 0; b < problem->nblocks; b++ ) 1322 { 1323 SCIP_CALL( SCIPallocBufferArray(scip, &((problem->blocks[b]).slackspos), blocktolinkvars[b].size * (nblocks - 1)) ); 1324 SCIP_CALL( SCIPallocBufferArray(scip, &((problem->blocks[b]).slacksneg), blocktolinkvars[b].size * (nblocks - 1)) ); 1325 SCIP_CALL( SCIPallocBufferArray(scip, &((problem->blocks[b]).couplingcons), blocktolinkvars[b].size * (nblocks - 1)) ); 1326 } 1327 1328 SCIP_CALL( SCIPallocBufferArray(scip, &tmpcouplingvars, COUPLINGSIZE) ); 1329 SCIP_CALL( SCIPallocBufferArray(scip, &tmpcouplingcoef, COUPLINGSIZE) ); 1330 tmpcouplingcoef[0] = 1.0; 1331 tmpcouplingcoef[1] = 1.0; 1332 tmpcouplingcoef[2] = -1.0; 1333 1334 /* count hashtable entries */ 1335 nentries = 0; 1336 for( b = 0; b < problem->nblocks; b++ ) 1337 { 1338 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1339 { 1340 int linkvaridx = blocktolinkvars[b].indexes[i]; 1341 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1342 { 1343 if( linkvartoblocks[linkvaridx].indexes[k] != b ) 1344 nentries++; 1345 } 1346 } 1347 } 1348 1349 SCIP_CALL( SCIPallocBufferArray(scip, &blockinfolist, nentries) ); 1350 SCIP_CALL( SCIPhashtableCreate(&htable, SCIPblkmem(scip), 1, SCIPhashGetKeyStandard, indexesEqual, indexesHashval, (void*) scip) ); 1351 blockinfolistfill = 0; 1352 1353 /* extend submips */ 1354 SCIPdebugMsg(scip, "Extending %d block models\n", problem->nblocks); 1355 for( b = 0; b < problem->nblocks; b++ ) 1356 { 1357 SCIP_VAR** blockvars; 1358 int nblockvars; 1359 1360 blockvars = SCIPgetVars((problem->blocks[b]).subscip); 1361 nblockvars = SCIPgetNVars((problem->blocks[b]).subscip); 1362 1363 /* set objective function of each block to zero */ 1364 for( i = 0; i < nblockvars; i++ ) 1365 { 1366 SCIP_CALL( SCIPchgVarObj((problem->blocks[b]).subscip, blockvars[i], 0.0) ); 1367 } 1368 1369 /* add two slack variables for each linking variable in block */ 1370 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1371 { 1372 int linkvaridx; 1373 linkvaridx = blocktolinkvars[b].indexes[i]; 1374 1375 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1376 { 1377 int b2; 1378 b2 = linkvartoblocks[linkvaridx].indexes[k]; 1379 1380 /* handle different blocks with common linking variable */ 1381 if( b2 != b ) 1382 { 1383 BLOCKINFO* binfo; 1384 binfo = &blockinfolist[blockinfolistfill]; 1385 blockinfolistfill++; 1386 binfo->block = b; 1387 binfo->otherblock = b2; 1388 binfo->linkvaridx = linkvaridx; 1389 binfo->linkvar = SCIPfindVar((problem->blocks[b]).subscip, SCIPvarGetName(linkvars[linkvaridx])); 1390 j = (problem->blocks[b]).ncoupling; 1391 1392 /* create positive slack variable */ 1393 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_slackpos_block_%d", SCIPvarGetName(linkvars[linkvaridx]), b2); 1394 (problem->blocks[b]).slackspos[j] = NULL; 1395 SCIP_CALL( SCIPcreateVarBasic((problem->blocks[b]).subscip, 1396 &((problem->blocks[b]).slackspos[j]), name, 1397 0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) ); 1398 SCIP_CALL( SCIPaddVar((problem->blocks[b]).subscip, (problem->blocks[b]).slackspos[j]) ); 1399 assert((problem->blocks[b]).slackspos[j] != NULL); 1400 binfo->slackposobjcoeff = 1.0; 1401 binfo->slackposvar = (problem->blocks[b]).slackspos[j]; 1402 1403 /* create negative slack variable */ 1404 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_slackneg_block_%d", SCIPvarGetName(linkvars[linkvaridx]), b2); 1405 (problem->blocks[b]).slacksneg[j] = NULL; 1406 SCIP_CALL( SCIPcreateVarBasic((problem->blocks[b]).subscip, 1407 &((problem->blocks[b]).slacksneg[j]), name, 1408 0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_CONTINUOUS) ); 1409 SCIP_CALL( SCIPaddVar((problem->blocks[b]).subscip, (problem->blocks[b]).slacksneg[j]) ); 1410 assert((problem->blocks[b]).slacksneg[j] != NULL); 1411 binfo->slacknegobjcoeff = 1.0; 1412 binfo->slacknegvar = (problem->blocks[b]).slacksneg[j]; 1413 1414 /* fill variables for linking constraint */ 1415 tmpcouplingvars[0] = binfo->linkvar; 1416 tmpcouplingvars[1] = binfo->slackposvar; 1417 tmpcouplingvars[2] = binfo->slacknegvar; 1418 1419 /* create linking constraint */ 1420 (void)SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_coupling_block_%d", 1421 SCIPvarGetName(linkvars[linkvaridx]), b2); 1422 (problem->blocks[b]).couplingcons[j] = NULL; 1423 1424 /* create linking constraint with initial side equal to zero (or lower bound of linking variable) */ 1425 if( heurtiming & SCIP_HEURTIMING_BEFORENODE ) 1426 { 1427 SCIP_Real initval; 1428 SCIP_Real lb; 1429 1430 lb = SCIPvarGetLbOriginal(binfo->linkvar); 1431 initval = MAX(lb, 0.0); 1432 1433 SCIP_CALL( SCIPcreateConsBasicLinear((problem->blocks[b]).subscip, &((problem->blocks[b]).couplingcons[j]), 1434 name, COUPLINGSIZE, tmpcouplingvars, tmpcouplingcoef, initval, initval) ); 1435 1436 /* set initial value of linking variable */ 1437 binfo->linkvarval = initval; 1438 } 1439 1440 /* create linking constraint with initial side equal to LP solution (rounded if variable is integer) */ 1441 if( heurtiming & SCIP_HEURTIMING_AFTERNODE ) 1442 { 1443 SCIP_Real initval; 1444 1445 initval = SCIPvarGetLPSol(linkvars[linkvaridx]); 1446 if( SCIPvarGetType(binfo->linkvar) != SCIP_VARTYPE_CONTINUOUS ) 1447 initval = SCIPround(scip, initval); 1448 1449 SCIP_CALL( SCIPcreateConsBasicLinear((problem->blocks[b]).subscip, &((problem->blocks[b]).couplingcons[j]), 1450 name, COUPLINGSIZE, tmpcouplingvars, tmpcouplingcoef, initval, initval) ); 1451 1452 /* set initial value of linking variable */ 1453 binfo->linkvarval = initval; 1454 } 1455 1456 SCIP_CALL( SCIPaddCons((problem->blocks[b]).subscip, (problem->blocks[b]).couplingcons[j]) ); 1457 assert((problem->blocks[b]).couplingcons[j] != NULL); 1458 binfo->couplingCons = (problem->blocks[b]).couplingcons[j]; 1459 1460 (problem->blocks[b]).ncoupling++; 1461 1462 /* feed hashtable */ 1463 SCIP_CALL( SCIPhashtableSafeInsert(htable, (void*) binfo) ); 1464 } 1465 } 1466 } 1467 } 1468 assert(nentries == SCIPhashtableGetNElements(htable)); 1469 1470 #ifdef PADM_WRITE_PROBLEMS 1471 /* write extended submips */ 1472 for( b = 0; b < problem->nblocks; b++ ) 1473 { 1474 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "extended_block_%d.lp", b); 1475 SCIP_CALL( SCIPwriteOrigProblem((problem->blocks[b]).subscip, name, NULL, FALSE) ); 1476 } 1477 #endif 1478 1479 /* determine threshold for penalty coefficients via maximum norm */ 1480 slackthreshold = SCIP_REAL_MIN; 1481 for( i = 0; i < nvars; i++ ) 1482 { 1483 SCIP_Real obj; 1484 1485 obj = REALABS(SCIPvarGetObj(vars[i])); 1486 if( obj > slackthreshold ) 1487 slackthreshold = obj; 1488 } 1489 1490 /* ------------------------------------------------------------------------------------------------- */ 1491 1492 /* check whether there is enough time left */ 1493 SCIP_CALL( getTimeLeft(scip, &timeleft) ); 1494 if( timeleft <= 0 ) 1495 { 1496 SCIPdebugMsg(scip, "no time left\n"); 1497 goto TERMINATE; 1498 } 1499 1500 SCIPdebugMsg(scip, "Starting iterations\n"); 1501 SCIPdebugMsg(scip, "PIt\tADMIt\tSlacks\tInfo\n"); 1502 1503 piter = 0; 1504 increasedslacks = 0; 1505 (void) SCIPsnprintf(info, SCIP_MAXSTRLEN, "-"); 1506 solved = FALSE; 1507 istimeleft = TRUE; 1508 1509 /* Penalty loop */ 1510 while( !solved && piter < heurdata->penaltyiterations && istimeleft ) 1511 { 1512 piter++; 1513 solutionsdiffer = TRUE; 1514 aiter = 0; 1515 1516 /* Alternating direction method loop */ 1517 while( solutionsdiffer && aiter < heurdata->admiterations && istimeleft ) 1518 { 1519 aiter++; 1520 solutionsdiffer = FALSE; 1521 SCIPdebugMsg(scip, "%d\t%d\t%d\t%s\n", piter, aiter, increasedslacks, info); 1522 1523 /* Loop through the blocks and solve each sub-SCIP, potentially multiple times */ 1524 for( b = 0; b < problem->nblocks; b++ ) 1525 { 1526 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1527 { 1528 int linkvaridx; 1529 linkvaridx = blocktolinkvars[b].indexes[i]; 1530 1531 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1532 { 1533 int b2; 1534 b2 = linkvartoblocks[linkvaridx].indexes[k]; 1535 1536 if( b2 != b ) 1537 { 1538 BLOCKINFO binfo; 1539 BLOCKINFO* binfoout; 1540 BLOCKINFO binfo2; 1541 BLOCKINFO* binfo2out; 1542 1543 SCIP_CONS* couplingcons; 1544 SCIP_Real newrhs; 1545 1546 binfo.block = b; 1547 binfo.otherblock = b2; 1548 binfo.linkvaridx = linkvaridx; 1549 1550 binfoout = (BLOCKINFO*) SCIPhashtableRetrieve(htable, (void *)&binfo); 1551 assert(binfoout != NULL); 1552 couplingcons = binfoout->couplingCons; 1553 1554 /* interchange blocks b and b2 for getting new right hand side */ 1555 binfo2.block = b2; 1556 binfo2.otherblock = b; 1557 binfo2.linkvaridx = linkvaridx; 1558 binfo2out = (BLOCKINFO*) SCIPhashtableRetrieve(htable, (void*) &binfo2); 1559 assert(binfo2out != NULL); 1560 newrhs = binfo2out->linkvarval; 1561 1562 /* change side of coupling constraint equation with linking variable value of the other block */ 1563 SCIP_CALL( SCIPchgLhsLinear((problem->blocks[b]).subscip, couplingcons, newrhs) ); 1564 SCIP_CALL( SCIPchgRhsLinear((problem->blocks[b]).subscip, couplingcons, newrhs) ); 1565 1566 /* change penalty coefficients of slack variables */ 1567 SCIP_CALL( SCIPchgVarObj((problem->blocks[b]).subscip, binfoout->slackposvar, binfoout->slackposobjcoeff) ); 1568 SCIP_CALL( SCIPchgVarObj((problem->blocks[b]).subscip, binfoout->slacknegvar, binfoout->slacknegobjcoeff) ); 1569 } 1570 } 1571 } 1572 1573 /* increase slack penalty coeffs until each subproblem can be solved to optimality */ 1574 do 1575 { 1576 SCIP_Longint nnodes; 1577 int iteration; 1578 1579 #ifdef PADM_WRITE_PROBLEMS 1580 SCIPdebugMsg(scip, "write subscip of block %d in piter=%d and aiter=%d\n", b, piter, aiter); 1581 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "blockproblem_%d_%d_%d.lp", b, piter, aiter); 1582 SCIP_CALL( SCIPwriteOrigProblem((problem->blocks[b]).subscip, name, NULL, FALSE) ); 1583 #endif 1584 1585 SCIP_CALL( SCIPsetRealParam((problem->blocks[b]).subscip, "limits/gap", gap) ); 1586 1587 /* reuse old solution if available */ 1588 SCIP_CALL( reuseSolution((problem->blocks[b]).subscip, &problem->blocks[b]) ); 1589 1590 /* update time and memory limit of subproblem */ 1591 SCIP_CALL( SCIPcopyLimits(scip, (problem->blocks[b]).subscip) ); 1592 1593 /* stop if there are not enough nodes left */ 1594 if( nodesleft < heurdata->minnodes ) 1595 { 1596 SCIPdebugMsg(scip, "Node limit reached.\n"); 1597 goto TERMINATE; 1598 } 1599 1600 /* update node limit of subproblem 1601 * in the first iterations we have a smaller node limit 1602 */ 1603 iteration = ((piter - 1) * heurdata->admiterations) + aiter; 1604 nnodes = (SCIP_Longint)SCIPceil(scip, (problem->blocks[b]).size * nodesleft * ( 1 - pow(heurdata->nodefac, (double)iteration) )); 1605 nnodes = MAX( heurdata->minnodes, nnodes ); 1606 SCIP_CALL( SCIPsetLongintParam((problem->blocks[b]).subscip, "limits/nodes", nnodes) ); 1607 1608 /* solve block 1609 * 1610 * errors in solving the subproblem should not kill the overall solving process; 1611 * hence, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. 1612 */ 1613 SCIP_CALL_ABORT( SCIPsolve((problem->blocks[b]).subscip) ); 1614 status = SCIPgetStatus((problem->blocks[b]).subscip); 1615 1616 /* subtract used nodes from the total nodelimit */ 1617 nodesleft -= (SCIP_Longint)SCIPceil(scip, SCIPgetNNodes((problem->blocks[b]).subscip) * (problem->blocks[b]).size); 1618 1619 /* check solution if one of the four cases occurs 1620 * - solution is optimal 1621 * - solution reached gaplimit 1622 * - node limit reached with at least one feasible solution 1623 * - time limit is reached but best solution needs no slack variables (no dual solution available) 1624 */ 1625 if( status == SCIP_STATUS_OPTIMAL || status == SCIP_STATUS_GAPLIMIT || 1626 (status == SCIP_STATUS_NODELIMIT && SCIPgetNSols((problem->blocks[b]).subscip) > 0) || 1627 (status == SCIP_STATUS_TIMELIMIT && SCIPgetNSols((problem->blocks[b]).subscip) > 0 && 1628 SCIPisEQ(scip, SCIPgetSolOrigObj((problem->blocks[b]).subscip, SCIPgetBestSol((problem->blocks[b]).subscip)), 0.0) ) ) 1629 { 1630 SCIPdebugMsg(scip, "Block is optimal or reached gaplimit or nodelimit.\n"); 1631 1632 if( status == SCIP_STATUS_TIMELIMIT ) 1633 { 1634 SCIPdebugMsg(scip, "Block reached time limit with at least one feasible solution.\n"); 1635 istimeleft = FALSE; 1636 } 1637 1638 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1639 { 1640 int linkvaridx; 1641 linkvaridx = blocktolinkvars[b].indexes[i]; 1642 1643 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1644 { 1645 int b2; 1646 b2 = linkvartoblocks[linkvaridx].indexes[k]; 1647 1648 if( b2 != b ) 1649 { 1650 SCIP_SOL* sol; 1651 BLOCKINFO binfo; 1652 BLOCKINFO* binfoout; 1653 SCIP_VAR* var; 1654 SCIP_Real val; 1655 1656 binfo.block = b; 1657 binfo.otherblock = b2; 1658 binfo.linkvaridx = linkvaridx; 1659 binfoout = (BLOCKINFO *)SCIPhashtableRetrieve(htable, (void *)&binfo); 1660 assert(binfoout != NULL); 1661 1662 sol = SCIPgetBestSol((problem->blocks[b]).subscip); 1663 assert(sol != NULL); 1664 var = binfoout->linkvar; 1665 val = SCIPgetSolVal((problem->blocks[b]).subscip, sol, var); 1666 1667 if( !EPSEQ(binfoout->linkvarval, val, SCIP_DEFAULT_EPSILON) ) 1668 solutionsdiffer = TRUE; 1669 1670 binfoout->linkvarval = val; 1671 } 1672 } 1673 } 1674 } 1675 else if( status == SCIP_STATUS_UNBOUNDED ) 1676 { 1677 SCIPdebugMsg(scip, "Block is unbounded.\n"); 1678 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1679 { 1680 int linkvaridx; 1681 linkvaridx = blocktolinkvars[b].indexes[i]; 1682 1683 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1684 { 1685 int b2; 1686 b2 = linkvartoblocks[linkvaridx].indexes[k]; 1687 1688 if( b2 != b ) 1689 { 1690 BLOCKINFO binfo; 1691 BLOCKINFO* binfoout; 1692 1693 binfo.block = b; 1694 binfo.otherblock = b2; 1695 binfo.linkvaridx = linkvaridx; 1696 binfoout = (BLOCKINFO*) SCIPhashtableRetrieve(htable, (void*) &binfo); 1697 assert(binfoout != NULL); 1698 1699 /* increase penalty coefficients to obtain a bounded subproblem */ 1700 binfoout->slackposobjcoeff *= 10.0; 1701 binfoout->slacknegobjcoeff *= 10.0; 1702 SCIP_CALL( SCIPchgVarObj((problem->blocks[b]).subscip, binfoout->slackposvar, binfoout->slackposobjcoeff) ); 1703 SCIP_CALL( SCIPchgVarObj((problem->blocks[b]).subscip, binfoout->slacknegvar, binfoout->slacknegobjcoeff) ); 1704 } 1705 } 1706 } 1707 } 1708 else if( status == SCIP_STATUS_TIMELIMIT ) 1709 { 1710 SCIPdebugMsg(scip, "Block reached time limit. No optimal solution available.\n"); 1711 goto TERMINATE; 1712 } 1713 else 1714 { 1715 SCIPdebugMsg(scip, "Block solving status %d not supported\n", status); 1716 goto TERMINATE; 1717 } 1718 1719 /* free solving data in order to change problem */ 1720 SCIP_CALL( SCIPfreeTransform((problem->blocks[b]).subscip) ); 1721 } 1722 while( status != SCIP_STATUS_OPTIMAL && status != SCIP_STATUS_GAPLIMIT && 1723 !(status == SCIP_STATUS_NODELIMIT && SCIPgetNSols((problem->blocks[b]).subscip) > 0) && 1724 !(status == SCIP_STATUS_TIMELIMIT && SCIPgetNSols((problem->blocks[b]).subscip) > 0 && 1725 SCIPisEQ(scip, SCIPgetSolOrigObj((problem->blocks[b]).subscip, SCIPgetBestSol((problem->blocks[b]).subscip)), 0.0) ) ); 1726 } 1727 } 1728 1729 /* check wether problem has been solved and if not update penalty coeffs */ 1730 doscaling = FALSE; 1731 solved = TRUE; 1732 increasedslacks = 0; 1733 maxpenalty = SCIP_REAL_MIN; 1734 for( b = 0; b < problem->nblocks; b++ ) 1735 { 1736 for( i = 0; i < blocktolinkvars[b].size; i++ ) 1737 { 1738 int linkvaridx; 1739 linkvaridx = blocktolinkvars[b].indexes[i]; 1740 1741 for( k = 0; k < linkvartoblocks[linkvaridx].size; k++ ) 1742 { 1743 int b2; 1744 b2 = linkvartoblocks[linkvaridx].indexes[k]; 1745 1746 if( b2 != b ) 1747 { 1748 SCIP_SOL* sol; 1749 BLOCKINFO binfo; 1750 BLOCKINFO* binfoout; 1751 SCIP_Real slackposval; 1752 SCIP_Real slacknegval; 1753 1754 binfo.block = b; 1755 binfo.otherblock = b2; 1756 binfo.linkvaridx = linkvaridx; 1757 binfoout = (BLOCKINFO*) SCIPhashtableRetrieve(htable, (void*) &binfo); 1758 assert(binfoout != NULL); 1759 1760 sol = SCIPgetBestSol((problem->blocks[b]).subscip); 1761 slackposval = SCIPgetSolVal((problem->blocks[b]).subscip, sol, binfoout->slackposvar); 1762 slacknegval = SCIPgetSolVal((problem->blocks[b]).subscip, sol, binfoout->slacknegvar); 1763 1764 /* increase penalty coefficient of positive slack variable */ 1765 if( SCIPisGT(scip, slackposval, 0.0) ) 1766 { 1767 binfoout->slackposobjcoeff *= 10.0; 1768 1769 if( binfoout->slackposobjcoeff > slackthreshold ) 1770 doscaling = TRUE; 1771 1772 if( binfoout->slackposobjcoeff > maxpenalty ) 1773 maxpenalty = binfoout->slackposobjcoeff; 1774 1775 solved = FALSE; 1776 increasedslacks++; 1777 } 1778 1779 /* increase penalty coefficient of negative slack variable */ 1780 if( SCIPisGT(scip, slacknegval, 0.0) ) 1781 { 1782 binfoout->slacknegobjcoeff *= 10.0; 1783 1784 if( binfoout->slacknegobjcoeff > slackthreshold ) 1785 doscaling = TRUE; 1786 1787 if( binfoout->slacknegobjcoeff > maxpenalty ) 1788 maxpenalty = binfoout->slacknegobjcoeff; 1789 1790 solved = FALSE; 1791 increasedslacks++; 1792 } 1793 } 1794 } 1795 } 1796 } 1797 1798 /* should sigmoid scaling be applied to the penalty parameters? */ 1799 if( doscaling && heurdata->scaling ) 1800 { 1801 SCIPdebugMsg(scip, "rescale penalty parameters\n"); 1802 1803 /* reset counter */ 1804 increasedslacks = 0; 1805 1806 /* rescale penalty parameters */ 1807 SCIP_CALL( scalePenalties(problem, linkvartoblocks, blocktolinkvars, htable, maxpenalty) ); 1808 } 1809 1810 /* adapt in some cases the gap parameter */ 1811 if( (aiter == 1 && solutionsdiffer == FALSE) || (doscaling && heurdata->scaling) ) 1812 { 1813 SCIP_Real mingap = 0.001; //todo 1814 SCIP_Real newgap = MAX(gap * 0.5, mingap); 1815 1816 if( newgap >= mingap ) 1817 { 1818 if( doscaling && heurdata->scaling ) 1819 (void) SCIPsnprintf(info, SCIP_MAXSTRLEN, "scale, %f", newgap); 1820 else 1821 (void) SCIPsnprintf(info, SCIP_MAXSTRLEN, "%f", newgap); 1822 1823 gap = newgap; 1824 } 1825 } 1826 1827 /* free solution process data */ 1828 if( !solved ) 1829 for( b = 0; b < problem->nblocks; b++ ) 1830 { 1831 SCIP_CALL( SCIPfreeTransform((problem->blocks[b]).subscip) ); 1832 } 1833 } 1834 1835 /* copy solution if present */ 1836 if( solved ) 1837 { 1838 SCIP_SOL* newsol; 1839 SCIP_Real* blocksolvals; 1840 1841 assert(increasedslacks == 0); 1842 1843 SCIP_CALL( SCIPallocBufferArray(scip, &blocksolvals, nvars) ); 1844 SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) ); 1845 1846 for( b = 0; b < problem->nblocks; b++ ) 1847 { 1848 SCIP_SOL* blocksol; 1849 SCIP_VAR** blockvars; 1850 int nblockvars; 1851 1852 /* get solution of block variables (without slack variables) */ 1853 blocksol = SCIPgetBestSol((problem->blocks[b]).subscip); 1854 assert(blocksol != NULL); 1855 blockvars = (problem->blocks[b]).subvars; 1856 nblockvars = (problem->blocks[b]).nsubvars; 1857 SCIP_CALL( SCIPgetSolVals((problem->blocks[b]).subscip, blocksol, nblockvars, blockvars, blocksolvals) ); 1858 1859 for( i = 0; i < nblockvars; i++ ) 1860 { 1861 SCIP_VAR* origvar; 1862 SCIP_Real solval; 1863 1864 origvar = SCIPfindVar(scip, SCIPvarGetName(blockvars[i])); 1865 solval = blocksolvals[i]; 1866 SCIP_CALL_ABORT( SCIPsetSolVal(scip, newsol, origvar, solval) ); 1867 } 1868 } 1869 1870 /* treat variables with no constraints; set value of variable to bound */ 1871 for( i = 0; i < numlinkvars; i++ ) 1872 { 1873 if( varonlyobj[i] ) 1874 { 1875 SCIP_Real fixedvalue; 1876 if( SCIPvarGetObj(linkvars[i]) < 0 ) 1877 { 1878 fixedvalue = SCIPvarGetUbLocal(linkvars[i]); 1879 if( SCIPisInfinity(scip, fixedvalue) ) 1880 break; // todo: maybe we should return the status UNBOUNDED instead 1881 } 1882 else 1883 { 1884 fixedvalue = SCIPvarGetLbLocal(linkvars[i]); 1885 if( SCIPisInfinity(scip, fixedvalue) ) 1886 break; // todo: maybe we should return the status UNBOUNDED instead 1887 } 1888 SCIP_CALL_ABORT( SCIPsetSolVal(scip, newsol, linkvars[i], fixedvalue) ); 1889 } 1890 } 1891 1892 SCIPdebugMsg(scip, "Objective value %.2f\n", SCIPgetSolOrigObj(scip, newsol)); 1893 1894 /* fix linking variables and reoptimize with original objective function */ 1895 if( heurdata->reoptimize ) 1896 { 1897 SCIP_SOL* improvedsol = NULL; 1898 SCIP_CALL( reoptimize(scip, heur, newsol, vars, nvars, linkvars, numlinkvars, &improvedsol, &success) ); 1899 assert(improvedsol != NULL || success == FALSE); 1900 1901 if( success ) 1902 { 1903 SCIP_CALL( SCIPtrySolFree(scip, &improvedsol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) ); 1904 if( !success ) 1905 { 1906 SCIPdebugMsg(scip, "Reoptimizing solution failed\n"); 1907 } 1908 else 1909 { 1910 SCIPdebugMsg(scip, "Reoptimizing solution successful\n"); 1911 *result = SCIP_FOUNDSOL; 1912 } 1913 } 1914 } 1915 1916 /* if reoptimization is turned off or reoptimization found no solution, try initial solution */ 1917 if( *result != SCIP_FOUNDSOL ) 1918 { 1919 SCIP_CALL( SCIPtrySolFree(scip, &newsol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) ); 1920 if( !success ) 1921 { 1922 SCIPdebugMsg(scip, "Solution copy failed\n"); 1923 } 1924 else 1925 { 1926 SCIPdebugMsg(scip, "Solution copy successful\n"); 1927 *result = SCIP_FOUNDSOL; 1928 } 1929 } 1930 else 1931 { 1932 SCIP_CALL( SCIPfreeSol(scip, &newsol) ); 1933 } 1934 1935 SCIPfreeBufferArray(scip, &blocksolvals); 1936 } 1937 else 1938 { 1939 SCIPdebugMsg(scip, "maximum number of penalty loops reached\n"); 1940 *result = SCIP_DIDNOTFIND; 1941 } 1942 1943 TERMINATE: 1944 /* release variables, constraints and free memory */ 1945 if( problem != NULL ) 1946 { 1947 for( b = 0; b < problem->nblocks; b++ ) 1948 { 1949 BLOCK curr_block = problem->blocks[b]; 1950 for( i = 0; i < (problem->blocks[b]).ncoupling; i++ ) 1951 { 1952 SCIP_CALL( SCIPreleaseCons(curr_block.subscip, &curr_block.couplingcons[i]) ); 1953 SCIP_CALL( SCIPreleaseVar(curr_block.subscip, &curr_block.slackspos[i]) ); 1954 SCIP_CALL( SCIPreleaseVar(curr_block.subscip, &curr_block.slacksneg[i]) ); 1955 } 1956 } 1957 } 1958 1959 if( htable != NULL ) 1960 SCIPhashtableFree(&htable); 1961 1962 if( blockinfolist != NULL ) 1963 SCIPfreeBufferArray(scip, &blockinfolist); 1964 1965 if( tmpcouplingcoef != NULL ) 1966 SCIPfreeBufferArray(scip, &tmpcouplingcoef); 1967 1968 if( tmpcouplingvars != NULL ) 1969 SCIPfreeBufferArray(scip, &tmpcouplingvars); 1970 1971 if( problem != NULL ) 1972 { 1973 for( b = problem->nblocks - 1; b >= 0; b-- ) 1974 { 1975 if( problem->blocks[b].couplingcons != NULL ) 1976 { 1977 SCIPfreeBufferArray(scip, &problem->blocks[b].couplingcons); 1978 SCIPfreeBufferArray(scip, &problem->blocks[b].slacksneg); 1979 SCIPfreeBufferArray(scip, &problem->blocks[b].slackspos); 1980 } 1981 } 1982 } 1983 1984 if( varonlyobj != NULL ) 1985 SCIPfreeBufferArray(scip, &varonlyobj); 1986 1987 if( problem != NULL && blocktolinkvars != NULL ) 1988 { 1989 for( b = problem->nblocks -1; b >= 0; b-- ) 1990 { 1991 if( blocktolinkvars[b].indexes != NULL ) 1992 SCIPfreeBufferArray(scip, &(blocktolinkvars[b].indexes)); 1993 } 1994 } 1995 1996 if( linkvars != NULL ) 1997 SCIPfreeBufferArray(scip, &linkvars); 1998 1999 if( alllinkvartoblocks != NULL ) 2000 SCIPfreeBufferArray(scip, &alllinkvartoblocks); 2001 2002 if( blocktolinkvars != NULL ) 2003 SCIPfreeBufferArray(scip, &blocktolinkvars); 2004 2005 if( linkvartoblocks != NULL ) 2006 SCIPfreeBufferArray(scip, &linkvartoblocks); 2007 2008 if( assigneddecomp != NULL ) 2009 SCIPfreeDecomp(scip, &assigneddecomp); 2010 2011 if( consssize != NULL ) 2012 SCIPfreeBufferArray(scip, &consssize); 2013 2014 if( conslabels != NULL ) 2015 SCIPfreeBufferArray(scip, &conslabels); 2016 2017 if( varlabels != NULL ) 2018 SCIPfreeBufferArray(scip, &varlabels); 2019 2020 if( sortedconss != NULL ) 2021 SCIPfreeBufferArray(scip, &sortedconss); 2022 2023 if( problem != NULL ) 2024 { 2025 SCIP_CALL( freeProblem(&problem, nblocks) ); 2026 } 2027 2028 SCIPdebugMsg(scip, "Leave padm heuristic\n"); 2029 return SCIP_OKAY; 2030 } 2031 2032 /* 2033 * primal heuristic specific interface methods 2034 */ 2035 2036 /** creates the PADM primal heuristic and includes it in SCIP */ 2037 SCIP_RETCODE SCIPincludeHeurPADM( 2038 SCIP* scip /**< SCIP data structure */ 2039 ) 2040 { 2041 SCIP_HEURDATA* heurdata; 2042 SCIP_HEUR* heur = NULL; 2043 2044 /* create PADM primal heuristic data */ 2045 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 2046 2047 /* include primal heuristic */ 2048 2049 /* use SCIPincludeHeurBasic() plus setter functions if you want to set callbacks one-by-one and your code should 2050 * compile independent of new callbacks being added in future SCIP versions 2051 */ 2052 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 2053 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 2054 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecPADM, heurdata) ); 2055 2056 assert(heur != NULL); 2057 2058 /* set non fundamental callbacks via setter functions */ 2059 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyPADM) ); 2060 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreePADM) ); 2061 2062 /* add padm primal heuristic parameters */ 2063 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/maxnodes", 2064 "maximum number of nodes to regard in all subproblems", 2065 &heurdata->maxnodes, TRUE, DEFAULT_MAXNODES, 0LL, SCIP_LONGINT_MAX, NULL, NULL) ); 2066 2067 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/minnodes", 2068 "minimum number of nodes to regard in one subproblem", 2069 &heurdata->minnodes, TRUE, DEFAULT_MINNODES, 0LL, SCIP_LONGINT_MAX, NULL, NULL) ); 2070 2071 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/nodefac", 2072 "factor to control nodelimits of subproblems", &heurdata->nodefac, TRUE, DEFAULT_NODEFAC, 0.0, 0.99, NULL, NULL) ); 2073 2074 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/admiterations", 2075 "maximal number of ADM iterations in each penalty loop", &heurdata->admiterations, TRUE, DEFAULT_ADMIT, 1, 100, NULL, NULL) ); 2076 2077 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/penaltyiterations", 2078 "maximal number of penalty iterations", &heurdata->penaltyiterations, TRUE, DEFAULT_PENALTYIT, 1, 100000, NULL, NULL) ); 2079 2080 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/gap", 2081 "mipgap at start", &heurdata->gap, TRUE, DEFAULT_GAP, 0.0, 16.0, NULL, NULL) ); 2082 2083 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/reoptimize", 2084 "should the problem get reoptimized with the original objective function?", &heurdata->reoptimize, FALSE, TRUE, NULL, NULL) ); 2085 2086 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/scaling", 2087 "enable sigmoid rescaling of penalty parameters", &heurdata->scaling, TRUE, TRUE, NULL, NULL) ); 2088 2089 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/assignlinking", 2090 "should linking constraints be assigned?", &heurdata->assignlinking, FALSE, TRUE, NULL, NULL) ); 2091 2092 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/original", 2093 "should the original problem be used? This is only for testing and not recommended!", &heurdata->original, TRUE, FALSE, NULL, NULL) ); 2094 2095 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/timing", 2096 "should the heuristic run before or after the processing of the node? (0: before, 1: after, 2: both)", 2097 &heurdata->timing, FALSE, 0, 0, 2, NULL, NULL) ); 2098 2099 return SCIP_OKAY; 2100 } 2101