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_multistart.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief multistart heuristic for convex and nonconvex MINLPs 28 * @author Benjamin Mueller 29 */ 30 31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 32 33 #include "blockmemshell/memory.h" 34 #include "scip/scip_expr.h" 35 #include "scip/pub_expr.h" 36 #include "scip/heur_multistart.h" 37 #include "scip/heur_subnlp.h" 38 #include "scip/pub_heur.h" 39 #include "scip/pub_message.h" 40 #include "scip/pub_misc.h" 41 #include "scip/pub_misc_sort.h" 42 #include "scip/pub_nlp.h" 43 #include "scip/pub_var.h" 44 #include "scip/scip_general.h" 45 #include "scip/scip_heur.h" 46 #include "scip/scip_mem.h" 47 #include "scip/scip_message.h" 48 #include "scip/scip_nlp.h" 49 #include "scip/scip_nlpi.h" 50 #include "scip/scip_numerics.h" 51 #include "scip/scip_param.h" 52 #include "scip/scip_prob.h" 53 #include "scip/scip_randnumgen.h" 54 #include "scip/scip_sol.h" 55 #include "scip/scip_timing.h" 56 #include <string.h> 57 58 59 #define HEUR_NAME "multistart" 60 #define HEUR_DESC "multistart heuristic for convex and nonconvex MINLPs" 61 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_LNS 62 #define HEUR_PRIORITY -2100000 63 #define HEUR_FREQ 0 64 #define HEUR_FREQOFS 0 65 #define HEUR_MAXDEPTH -1 66 #define HEUR_TIMING SCIP_HEURTIMING_AFTERNODE 67 #define HEUR_USESSUBSCIP TRUE /**< does the heuristic use a secondary SCIP instance? */ 68 69 #define DEFAULT_RANDSEED 131 /**< initial random seed */ 70 #define DEFAULT_NRNDPOINTS 100 /**< default number of generated random points per call */ 71 #define DEFAULT_MAXBOUNDSIZE 2e+4 /**< default maximum variable domain size for unbounded variables */ 72 #define DEFAULT_MAXITER 300 /**< default number of iterations to reduce the violation of a point */ 73 #define DEFAULT_MINIMPRFAC 0.05 /**< default minimum required improving factor to proceed in improvement of a point */ 74 #define DEFAULT_MINIMPRITER 10 /**< default number of iteration when checking the minimum improvement */ 75 #define DEFAULT_MAXRELDIST 0.15 /**< default maximum distance between two points in the same cluster */ 76 #define DEFAULT_GRADLIMIT 5e+6 /**< default limit for gradient computations for all improvePoint() calls */ 77 #define DEFAULT_MAXNCLUSTER 3 /**< default maximum number of considered clusters per heuristic call */ 78 #define DEFAULT_ONLYNLPS TRUE /**< should the heuristic run only on continuous problems? */ 79 80 #define MINFEAS -1e+4 /**< minimum feasibility for a point; used for filtering and improving 81 * feasibility */ 82 #define MINIMPRFAC 0.95 /**< improvement factor used to discard randomly generated points with a 83 * too large objective value */ 84 #define GRADCOSTFAC_LINEAR 1.0 /**< gradient cost factor for the number of linear variables */ 85 #define GRADCOSTFAC_NONLINEAR 3.0 /**< gradient cost factor for the number of nodes in nonlinear expression */ 86 87 /* 88 * Data structures 89 */ 90 91 /** primal heuristic data */ 92 struct SCIP_HeurData 93 { 94 int nrndpoints; /**< number of random points generated per execution call */ 95 SCIP_Real maxboundsize; /**< maximum variable domain size for unbounded variables */ 96 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */ 97 SCIP_HEUR* heursubnlp; /**< sub-NLP heuristic */ 98 99 int maxiter; /**< number of iterations to reduce the maximum violation of a point */ 100 SCIP_Real minimprfac; /**< minimum required improving factor to proceed in the improvement of a single point */ 101 int minimpriter; /**< number of iteration when checking the minimum improvement */ 102 103 SCIP_Real maxreldist; /**< maximum distance between two points in the same cluster */ 104 SCIP_Real gradlimit; /**< limit for gradient computations for all improvePoint() calls (0 for no limit) */ 105 int maxncluster; /**< maximum number of considered clusters per heuristic call */ 106 SCIP_Bool onlynlps; /**< should the heuristic run only on continuous problems? */ 107 }; 108 109 110 /* 111 * Local methods 112 */ 113 114 115 /** returns an unique index of a variable in the range of 0,..,SCIPgetNVars(scip)-1 */ 116 #ifndef NDEBUG 117 static 118 int getVarIndex( 119 SCIP_HASHMAP* varindex, /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 */ 120 SCIP_VAR* var /**< variable */ 121 ) 122 { 123 assert(varindex != NULL); 124 assert(var != NULL); 125 assert(SCIPhashmapExists(varindex, (void*)var)); 126 127 return SCIPhashmapGetImageInt(varindex, (void*)var); 128 } 129 #else 130 #define getVarIndex(varindex,var) (SCIPhashmapGetImageInt((varindex), (void*)(var))) 131 #endif 132 133 /** samples and stores random points; stores points which have a better objective value than the current incumbent 134 * solution 135 */ 136 static 137 SCIP_RETCODE sampleRandomPoints( 138 SCIP* scip, /**< SCIP data structure */ 139 SCIP_SOL** rndpoints, /**< array to store all random points */ 140 int nmaxrndpoints, /**< maximum number of random points to compute */ 141 SCIP_Real maxboundsize, /**< maximum variable domain size for unbounded variables */ 142 SCIP_RANDNUMGEN* randnumgen, /**< random number generator */ 143 SCIP_Real bestobj, /**< objective value in the transformed space of the current incumbent */ 144 int* nstored /**< pointer to store the number of randomly generated points */ 145 ) 146 { 147 SCIP_VAR** vars; 148 SCIP_SOL* sol; 149 SCIP_Real val; 150 SCIP_Real lb; 151 SCIP_Real ub; 152 int nvars; 153 int niter; 154 int i; 155 156 assert(scip != NULL); 157 assert(rndpoints != NULL); 158 assert(nmaxrndpoints > 0); 159 assert(maxboundsize > 0.0); 160 assert(randnumgen != NULL); 161 assert(nstored != NULL); 162 163 vars = SCIPgetVars(scip); 164 nvars = SCIPgetNVars(scip); 165 *nstored = 0; 166 167 SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) ); 168 169 for( niter = 0; niter < 3 * nmaxrndpoints && *nstored < nmaxrndpoints; ++niter ) 170 { 171 /* reset solution, in case the old one had infinite objective, which can give difficulties in updating the obj value */ 172 SCIP_CALL( SCIPclearSol(scip, sol) ); 173 174 for( i = 0; i < nvars; ++i ) 175 { 176 lb = MIN(SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])); /*lint !e666*/ 177 ub = MAX(SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])); /*lint !e666*/ 178 179 if( SCIPisFeasEQ(scip, lb, ub) ) 180 val = (lb + ub) / 2.0; 181 /* use a smaller domain for unbounded variables */ 182 else if( !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) ) 183 val = SCIPrandomGetReal(randnumgen, lb, ub); 184 else if( !SCIPisInfinity(scip, -lb) ) 185 val = lb + SCIPrandomGetReal(randnumgen, 0.0, maxboundsize); 186 else if( !SCIPisInfinity(scip, ub) ) 187 val = ub - SCIPrandomGetReal(randnumgen, 0.0, maxboundsize); 188 else 189 { 190 assert(SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub)); 191 val = SCIPrandomGetReal(randnumgen, -0.5*maxboundsize, 0.5*maxboundsize); 192 } 193 assert(SCIPisFeasGE(scip, val, lb) && SCIPisFeasLE(scip, val, ub)); 194 195 /* set solution value; round the sampled point for integer variables */ 196 if( SCIPvarGetType(vars[i]) < SCIP_VARTYPE_CONTINUOUS ) 197 val = SCIPfeasRound(scip, val); 198 SCIP_CALL( SCIPsetSolVal(scip, sol, vars[i], val) ); 199 } 200 201 /* add solution if it is good enough */ 202 if( SCIPisLE(scip, SCIPgetSolTransObj(scip, sol), bestobj) ) 203 { 204 SCIP_CALL( SCIPcreateSolCopy(scip, &rndpoints[*nstored], sol) ); 205 ++(*nstored); 206 } 207 } 208 assert(*nstored <= nmaxrndpoints); 209 SCIPdebugMsg(scip, "found %d randomly generated points\n", *nstored); 210 211 SCIP_CALL( SCIPfreeSol(scip, &sol) ); 212 213 return SCIP_OKAY; 214 } 215 216 /** computes the minimum feasibility of a given point; a negative value means that there is an infeasibility */ 217 static 218 SCIP_RETCODE getMinFeas( 219 SCIP* scip, /**< SCIP data structure */ 220 SCIP_NLROW** nlrows, /**< array containing all nlrows */ 221 int nnlrows, /**< total number of nlrows */ 222 SCIP_SOL* sol, /**< solution */ 223 SCIP_Real* minfeas /**< buffer to store the minimum feasibility */ 224 ) 225 { 226 SCIP_Real tmp; 227 int i; 228 229 assert(scip != NULL); 230 assert(sol != NULL); 231 assert(minfeas != NULL); 232 assert(nlrows != NULL); 233 assert(nnlrows > 0); 234 235 *minfeas = SCIPinfinity(scip); 236 237 for( i = 0; i < nnlrows; ++i ) 238 { 239 assert(nlrows[i] != NULL); 240 241 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrows[i], sol, &tmp) ); 242 *minfeas = MIN(*minfeas, tmp); 243 } 244 245 return SCIP_OKAY; 246 } 247 248 /** computes the gradient for a given point and nonlinear row */ 249 static 250 SCIP_RETCODE computeGradient( 251 SCIP* scip, /**< SCIP data structure */ 252 SCIP_NLROW* nlrow, /**< nonlinear row */ 253 SCIP_SOL* sol, /**< solution to compute the gradient for */ 254 SCIP_HASHMAP* varindex, /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 uniquely */ 255 SCIP_EXPRITER* exprit, /**< expression iterator that can be used */ 256 SCIP_Real* grad, /**< buffer to store the gradient; grad[varindex(i)] corresponds to SCIPgetVars(scip)[i] */ 257 SCIP_Real* norm /**< buffer to store ||grad||^2 */ 258 ) 259 { 260 SCIP_EXPR* expr; 261 SCIP_VAR* var; 262 int i; 263 264 assert(scip != NULL); 265 assert(nlrow != NULL); 266 assert(varindex != NULL); 267 assert(sol != NULL); 268 assert(norm != NULL); 269 270 BMSclearMemoryArray(grad, SCIPgetNVars(scip)); 271 *norm = 0.0; 272 273 /* linear part */ 274 for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ ) 275 { 276 var = SCIPnlrowGetLinearVars(nlrow)[i]; 277 assert(var != NULL); 278 assert(getVarIndex(varindex, var) >= 0 && getVarIndex(varindex, var) < SCIPgetNVars(scip)); 279 280 grad[getVarIndex(varindex, var)] += SCIPnlrowGetLinearCoefs(nlrow)[i]; 281 } 282 283 /* expression part */ 284 expr = SCIPnlrowGetExpr(nlrow); 285 286 if( expr != NULL ) 287 { 288 assert(exprit != NULL); 289 290 SCIP_CALL( SCIPevalExprGradient(scip, expr, sol, 0L) ); 291 292 /* TODO: change this when nlrows store the vars */ 293 SCIP_CALL( SCIPexpriterInit(exprit, expr, SCIP_EXPRITER_DFS, FALSE) ); 294 for( ; !SCIPexpriterIsEnd(exprit); expr = SCIPexpriterGetNext(exprit) ) /*lint !e441*/ /*lint !e440*/ 295 { 296 if( !SCIPisExprVar(scip, expr) ) 297 continue; 298 299 var = SCIPgetVarExprVar(expr); 300 assert(var != NULL); 301 assert(getVarIndex(varindex, var) >= 0 && getVarIndex(varindex, var) < SCIPgetNVars(scip)); 302 303 grad[getVarIndex(varindex, var)] += SCIPexprGetDerivative(expr); 304 } 305 } 306 307 /* compute ||grad||^2 */ 308 for( i = 0; i < SCIPgetNVars(scip); ++i ) 309 *norm += SQR(grad[i]); 310 311 return SCIP_OKAY; 312 } 313 314 /** use consensus vectors to improve feasibility for a given starting point */ 315 static 316 SCIP_RETCODE improvePoint( 317 SCIP* scip, /**< SCIP data structure */ 318 SCIP_NLROW** nlrows, /**< array containing all nlrows */ 319 int nnlrows, /**< total number of nlrows */ 320 SCIP_HASHMAP* varindex, /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 */ 321 SCIP_SOL* point, /**< random generated point */ 322 int maxiter, /**< maximum number of iterations */ 323 SCIP_Real minimprfac, /**< minimum required improving factor to proceed in the improvement of a single point */ 324 int minimpriter, /**< number of iteration when checking the minimum improvement */ 325 SCIP_Real* minfeas, /**< pointer to store the minimum feasibility */ 326 SCIP_Real* nlrowgradcosts, /**< estimated costs for each gradient computation */ 327 SCIP_Real* gradcosts /**< pointer to store the estimated gradient costs */ 328 ) 329 { 330 SCIP_VAR** vars; 331 SCIP_EXPRITER* exprit; 332 SCIP_Real* grad; 333 SCIP_Real* updatevec; 334 SCIP_Real lastminfeas; 335 int nvars; 336 int r; 337 int i; 338 339 assert(varindex != NULL); 340 assert(point != NULL); 341 assert(maxiter > 0); 342 assert(minfeas != NULL); 343 assert(nlrows != NULL); 344 assert(nnlrows > 0); 345 assert(nlrowgradcosts != NULL); 346 assert(gradcosts != NULL); 347 348 *gradcosts = 0.0; 349 350 SCIP_CALL( getMinFeas(scip, nlrows, nnlrows, point, minfeas) ); 351 #ifdef SCIP_DEBUG_IMPROVEPOINT 352 printf("start minfeas = %e\n", *minfeas); 353 #endif 354 355 /* stop since start point is feasible */ 356 if( !SCIPisFeasLT(scip, *minfeas, 0.0) ) 357 { 358 #ifdef SCIP_DEBUG_IMPROVEPOINT 359 printf("start point is feasible"); 360 #endif 361 return SCIP_OKAY; 362 } 363 364 lastminfeas = *minfeas; 365 vars = SCIPgetVars(scip); 366 nvars = SCIPgetNVars(scip); 367 368 SCIP_CALL( SCIPallocBufferArray(scip, &grad, nvars) ); 369 SCIP_CALL( SCIPallocBufferArray(scip, &updatevec, nvars) ); 370 SCIP_CALL( SCIPcreateExpriter(scip, &exprit) ); 371 372 /* main loop */ 373 for( r = 0; r < maxiter && SCIPisFeasLT(scip, *minfeas, 0.0); ++r ) 374 { 375 SCIP_Real feasibility; 376 SCIP_Real activity; 377 SCIP_Real nlrownorm; 378 SCIP_Real scale; 379 int nviolnlrows; 380 381 BMSclearMemoryArray(updatevec, nvars); 382 nviolnlrows = 0; 383 384 for( i = 0; i < nnlrows; ++i ) 385 { 386 int j; 387 388 SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrows[i], point, &feasibility) ); 389 390 /* do not consider non-violated constraints */ 391 if( SCIPisFeasGE(scip, feasibility, 0.0) ) 392 continue; 393 394 /* increase number of violated nlrows */ 395 ++nviolnlrows; 396 397 SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrows[i], point, &activity) ); 398 SCIP_CALL( computeGradient(scip, nlrows[i], point, varindex, exprit, grad, &nlrownorm) ); 399 400 /* update estimated costs for computing gradients */ 401 *gradcosts += nlrowgradcosts[i]; 402 403 /* stop if the gradient disappears at the current point */ 404 if( SCIPisZero(scip, nlrownorm) ) 405 { 406 #ifdef SCIP_DEBUG_IMPROVEPOINT 407 printf("gradient vanished at current point -> stop\n"); 408 #endif 409 goto TERMINATE; 410 } 411 412 /* compute -g(x_k) / ||grad(g)(x_k)||^2 for a constraint g(x_k) <= 0 */ 413 scale = -feasibility / nlrownorm; 414 if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) && SCIPisGT(scip, activity, SCIPnlrowGetRhs(nlrows[i])) ) 415 scale *= -1.0; 416 417 /* skip nonliner row if the scaler is too small or too large */ 418 if( SCIPisEQ(scip, scale, 0.0) || SCIPisHugeValue(scip, REALABS(scale)) ) 419 continue; 420 421 for( j = 0; j < nvars; ++j ) 422 updatevec[j] += scale * grad[j]; 423 } 424 425 /* if there are no violated rows, stop since start point is feasible */ 426 if( nviolnlrows == 0 ) 427 { 428 assert(updatevec[i] == 0.0); 429 return SCIP_OKAY; 430 } 431 432 for( i = 0; i < nvars; ++i ) 433 { 434 /* adjust point */ 435 updatevec[i] = SCIPgetSolVal(scip, point, vars[i]) + updatevec[i] / nviolnlrows; 436 updatevec[i] = MIN(updatevec[i], SCIPvarGetUbLocal(vars[i])); /*lint !e666*/ 437 updatevec[i] = MAX(updatevec[i], SCIPvarGetLbLocal(vars[i])); /*lint !e666*/ 438 439 SCIP_CALL( SCIPsetSolVal(scip, point, vars[i], updatevec[i]) ); 440 } 441 442 /* update feasibility */ 443 SCIP_CALL( getMinFeas(scip, nlrows, nnlrows, point, minfeas) ); 444 445 /* check stopping criterion */ 446 if( r % minimpriter == 0 && r > 0 ) 447 { 448 if( *minfeas <= MINFEAS 449 || (*minfeas-lastminfeas) / MAX(REALABS(*minfeas), REALABS(lastminfeas)) < minimprfac ) /*lint !e666*/ 450 break; 451 lastminfeas = *minfeas; 452 } 453 } 454 455 TERMINATE: 456 #ifdef SCIP_DEBUG_IMPROVEPOINT 457 printf("niter=%d minfeas=%e\n", r, *minfeas); 458 #endif 459 460 SCIPfreeExpriter(&exprit); 461 462 SCIPfreeBufferArray(scip, &updatevec); 463 SCIPfreeBufferArray(scip, &grad); 464 465 return SCIP_OKAY; 466 } 467 468 /** sorts points w.r.t their feasibilities; points with a feasibility which is too small (w.r.t. the geometric mean of 469 * all feasibilities) will be filtered out 470 */ 471 static 472 SCIP_RETCODE filterPoints( 473 SCIP* scip, /**< SCIP data structure */ 474 SCIP_SOL** points, /**< array containing improved points */ 475 SCIP_Real* feasibilities, /**< array containing feasibility for each point (sorted) */ 476 int npoints, /**< total number of points */ 477 int* nusefulpoints /**< pointer to store the total number of useful points */ 478 ) 479 { 480 SCIP_Real minfeas; 481 SCIP_Real meanfeas; 482 int i; 483 484 assert(points != NULL); 485 assert(feasibilities != NULL); 486 assert(npoints > 0); 487 assert(nusefulpoints != NULL); 488 489 /* sort points w.r.t their feasibilities; non-negative feasibility correspond to feasible points for the NLP */ 490 SCIPsortDownRealPtr(feasibilities, (void**)points, npoints); 491 minfeas = feasibilities[npoints - 1]; 492 493 /* check if all points are feasible */ 494 if( SCIPisFeasGE(scip, minfeas, 0.0) ) 495 { 496 *nusefulpoints = npoints; 497 return SCIP_OKAY; 498 } 499 500 *nusefulpoints = 0; 501 502 /* compute shifted geometric mean of feasibilities (shift value = 1 - minfeas) */ 503 meanfeas = 1.0; 504 for( i = 0; i < npoints; ++i ) 505 { 506 assert(feasibilities[i] - minfeas + 1.0 > 0.0); 507 meanfeas *= pow(feasibilities[i] - minfeas + 1.0, 1.0 / npoints); 508 } 509 meanfeas += minfeas - 1.0; 510 SCIPdebugMsg(scip, "meanfeas = %e\n", meanfeas); 511 512 /* keep all points with which have a feasibility not much below the geometric mean of infeasibilities */ 513 for( i = 0; i < npoints; ++i ) 514 { 515 if( SCIPisFeasLT(scip, feasibilities[i], 0.0) 516 && (feasibilities[i] <= 1.05 * meanfeas || SCIPisLE(scip, feasibilities[i], MINFEAS)) ) 517 break; 518 519 ++(*nusefulpoints); 520 } 521 522 return SCIP_OKAY; 523 } 524 525 /** returns the relative distance between two points; considers a smaller bounded domain for unbounded variables */ 526 static 527 SCIP_Real getRelDistance( 528 SCIP* scip, /**< SCIP data structure */ 529 SCIP_SOL* x, /**< first point */ 530 SCIP_SOL* y, /**< second point */ 531 SCIP_Real maxboundsize /**< maximum variable domain size for unbounded variables */ 532 ) 533 { 534 SCIP_VAR** vars; 535 int nvars; 536 SCIP_Real distance; 537 SCIP_Real solx; 538 SCIP_Real soly; 539 SCIP_Real lb; 540 SCIP_Real ub; 541 int i; 542 543 assert(x != NULL); 544 assert(y != NULL); 545 546 vars = SCIPgetVars(scip); 547 nvars = SCIPgetNVars(scip); 548 distance = 0.0; 549 550 if( nvars == 0 ) 551 return 0.0; 552 553 for( i = 0; i < nvars; ++i ) 554 { 555 lb = SCIPvarGetLbLocal(vars[i]); 556 ub = SCIPvarGetUbLocal(vars[i]); 557 solx = SCIPgetSolVal(scip, x, vars[i]); 558 soly = SCIPgetSolVal(scip, y, vars[i]); 559 560 /* adjust lower and upper bounds for unbounded variables*/ 561 if( SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub) ) 562 { 563 lb = -maxboundsize / 2.0; 564 ub = +maxboundsize / 2.0; 565 } 566 else if( SCIPisInfinity(scip, -lb) ) 567 { 568 lb = ub - maxboundsize; 569 } 570 else if( SCIPisInfinity(scip, ub) ) 571 { 572 ub = lb + maxboundsize; 573 } 574 575 /* project solution values to the variable domain */ 576 solx = MIN(MAX(solx, lb), ub); 577 soly = MIN(MAX(soly, lb), ub); 578 579 distance += REALABS(solx - soly) / MAX(1.0, ub - lb); 580 } 581 582 return distance / nvars; 583 } 584 585 /** cluster useful points with a greedy algorithm */ 586 static 587 SCIP_RETCODE clusterPointsGreedy( 588 SCIP* scip, /**< SCIP data structure */ 589 SCIP_SOL** points, /**< array containing improved points */ 590 int npoints, /**< total number of points */ 591 int* clusteridx, /**< array to store for each point the index of the cluster */ 592 int* ncluster, /**< pointer to store the total number of cluster */ 593 SCIP_Real maxboundsize, /**< maximum variable domain size for unbounded variables */ 594 SCIP_Real maxreldist, /**< maximum relative distance between any two points of the same cluster */ 595 int maxncluster /**< maximum number of clusters to compute */ 596 ) 597 { 598 int i; 599 600 assert(points != NULL); 601 assert(npoints > 0); 602 assert(clusteridx != NULL); 603 assert(ncluster != NULL); 604 assert(maxreldist >= 0.0); 605 assert(maxncluster >= 0); 606 607 /* initialize cluster indices */ 608 for( i = 0; i < npoints; ++i ) 609 clusteridx[i] = INT_MAX; 610 611 *ncluster = 0; 612 613 for( i = 0; i < npoints && (*ncluster < maxncluster); ++i ) 614 { 615 int j; 616 617 /* point is already assigned to a cluster */ 618 if( clusteridx[i] != INT_MAX ) 619 continue; 620 621 /* create a new cluster for i */ 622 clusteridx[i] = *ncluster; 623 624 for( j = i + 1; j < npoints; ++j ) 625 { 626 if( clusteridx[j] == INT_MAX && getRelDistance(scip, points[i], points[j], maxboundsize) <= maxreldist ) 627 clusteridx[j] = *ncluster; 628 } 629 630 ++(*ncluster); 631 } 632 633 #ifndef NDEBUG 634 for( i = 0; i < npoints; ++i ) 635 { 636 assert(clusteridx[i] >= 0); 637 assert(clusteridx[i] < *ncluster || clusteridx[i] == INT_MAX); 638 } 639 #endif 640 641 return SCIP_OKAY; 642 } 643 644 /** calls the sub-NLP heuristic for a given cluster */ 645 static 646 SCIP_RETCODE solveNLP( 647 SCIP* scip, /**< SCIP data structure */ 648 SCIP_HEUR* heur, /**< multi-start heuristic */ 649 SCIP_HEUR* nlpheur, /**< pointer to NLP local search heuristics */ 650 SCIP_SOL** points, /**< array containing improved points */ 651 int npoints, /**< total number of points */ 652 SCIP_Bool* success /**< pointer to store if we could find a solution */ 653 ) 654 { 655 SCIP_VAR** vars; 656 SCIP_SOL* refpoint; 657 SCIP_RESULT nlpresult; 658 SCIP_Real val; 659 int nbinvars; 660 int nintvars; 661 int nvars; 662 int i; 663 664 assert(points != NULL); 665 assert(npoints > 0); 666 667 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) ); 668 *success = FALSE; 669 670 SCIP_CALL( SCIPcreateSol(scip, &refpoint, heur) ); 671 672 /* compute reference point */ 673 for( i = 0; i < nvars; ++i ) 674 { 675 int p; 676 677 val = 0.0; 678 679 for( p = 0; p < npoints; ++p ) 680 { 681 assert(points[p] != NULL); 682 val += SCIPgetSolVal(scip, points[p], vars[i]); 683 } 684 685 SCIP_CALL( SCIPsetSolVal(scip, refpoint, vars[i], val / npoints) ); 686 } 687 688 /* round point for sub-NLP heuristic */ 689 SCIP_CALL( SCIProundSol(scip, refpoint, success) ); 690 SCIPdebugMsg(scip, "rounding of refpoint successfully? %u\n", *success); 691 692 /* round variables manually if the locks did not allow us to round them */ 693 if( !(*success) ) 694 { 695 for( i = 0; i < nbinvars + nintvars; ++i ) 696 { 697 val = SCIPgetSolVal(scip, refpoint, vars[i]); 698 699 if( !SCIPisFeasIntegral(scip, val) ) 700 { 701 assert(SCIPisFeasIntegral(scip, SCIPvarGetLbLocal(vars[i]))); 702 assert(SCIPisFeasIntegral(scip, SCIPvarGetUbLocal(vars[i]))); 703 704 /* round and adjust value */ 705 val = SCIPround(scip, val); 706 val = MIN(val, SCIPvarGetUbLocal(vars[i])); /*lint !e666*/ 707 val = MAX(val, SCIPvarGetLbLocal(vars[i])); /*lint !e666*/ 708 assert(SCIPisFeasIntegral(scip, val)); 709 710 SCIP_CALL( SCIPsetSolVal(scip, refpoint, vars[i], val) ); 711 } 712 } 713 } 714 715 /* call sub-NLP heuristic */ 716 SCIP_CALL( SCIPapplyHeurSubNlp(scip, nlpheur, &nlpresult, refpoint, NULL) ); 717 SCIP_CALL( SCIPfreeSol(scip, &refpoint) ); 718 719 /* let sub-NLP heuristic decide whether the solution is feasible or not */ 720 *success = nlpresult == SCIP_FOUNDSOL; 721 722 return SCIP_OKAY; 723 } 724 725 /** recursive helper function to count the number of nodes in a sub-expr */ 726 static 727 int getExprSize( 728 SCIP_EXPR* expr /**< expression */ 729 ) 730 { 731 int sum; 732 int i; 733 734 assert(expr != NULL); 735 736 sum = 0; 737 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 738 { 739 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i]; 740 sum += getExprSize(child); 741 } 742 return 1 + sum; 743 } 744 745 /** main function of the multi-start heuristic (see @ref heur_multistart.h for more details); it consists of the 746 * following four steps: 747 * 748 * 1. sampling points in the current domain; for unbounded variables we use a bounded box 749 * 750 * 2. reduce infeasibility by using a gradient descent method 751 * 752 * 3. cluster points; filter points with a too large infeasibility 753 * 754 * 4. compute start point for each cluster and use it in the sub-NLP heuristic (@ref heur_subnlp.h) 755 */ 756 static 757 SCIP_RETCODE applyHeur( 758 SCIP* scip, /**< SCIP data structure */ 759 SCIP_HEUR* heur, /**< heuristic */ 760 SCIP_HEURDATA* heurdata, /**< heuristic data */ 761 SCIP_RESULT* result /**< pointer to store the result */ 762 ) 763 { 764 SCIP_NLROW** nlrows; 765 SCIP_SOL** points; 766 SCIP_HASHMAP* varindex; 767 SCIP_Real* feasibilities; 768 SCIP_Real* nlrowgradcosts; 769 int* clusteridx; 770 SCIP_Real gradlimit; 771 SCIP_Real bestobj; 772 int nusefulpoints; 773 int nrndpoints; 774 int ncluster; 775 int nnlrows; 776 int npoints; 777 int start; 778 int i; 779 780 assert(scip != NULL); 781 assert(heur != NULL); 782 assert(result != NULL); 783 assert(heurdata != NULL); 784 785 SCIPdebugMsg(scip, "call applyHeur()\n"); 786 787 nlrows = SCIPgetNLPNlRows(scip); 788 nnlrows = SCIPgetNNLPNlRows(scip); 789 bestobj = SCIPgetNSols(scip) > 0 ? MINIMPRFAC * SCIPgetSolTransObj(scip, SCIPgetBestSol(scip)) : SCIPinfinity(scip); 790 791 SCIP_CALL( SCIPallocBufferArray(scip, &points, heurdata->nrndpoints) ); 792 SCIP_CALL( SCIPallocBufferArray(scip, &nlrowgradcosts, nnlrows) ); 793 SCIP_CALL( SCIPallocBufferArray(scip, &feasibilities, heurdata->nrndpoints) ); 794 SCIP_CALL( SCIPallocBufferArray(scip, &clusteridx, heurdata->nrndpoints) ); 795 SCIP_CALL( SCIPhashmapCreate(&varindex, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 796 797 /* create an unique mapping of all variables to 0,..,SCIPgetNVars(scip)-1 */ 798 for( i = 0; i < SCIPgetNVars(scip); ++i ) 799 { 800 SCIP_CALL( SCIPhashmapInsertInt(varindex, (void*)SCIPgetVars(scip)[i], i) ); 801 } 802 803 /* compute estimated costs of computing a gradient for each nlrow */ 804 for( i = 0; i < nnlrows; ++i ) 805 { 806 nlrowgradcosts[i] = GRADCOSTFAC_LINEAR * SCIPnlrowGetNLinearVars(nlrows[i]); 807 if( SCIPnlrowGetExpr(nlrows[i]) != NULL ) 808 nlrowgradcosts[i] += GRADCOSTFAC_NONLINEAR * getExprSize(SCIPnlrowGetExpr(nlrows[i])); 809 } 810 811 /* 812 * 1. sampling points in the current domain; for unbounded variables we use a bounded box 813 */ 814 SCIP_CALL( sampleRandomPoints(scip, points, heurdata->nrndpoints, heurdata->maxboundsize, heurdata->randnumgen, 815 bestobj, &nrndpoints) ); 816 assert(nrndpoints >= 0); 817 818 if( nrndpoints == 0 ) 819 goto TERMINATE; 820 821 /* 822 * 2. improve points via consensus vectors 823 */ 824 gradlimit = heurdata->gradlimit == 0.0 ? SCIPinfinity(scip) : heurdata->gradlimit; 825 for( npoints = 0; npoints < nrndpoints && gradlimit >= 0 && !SCIPisStopped(scip); ++npoints ) 826 { 827 SCIP_Real gradcosts; 828 829 SCIP_CALL( improvePoint(scip, nlrows, nnlrows, varindex, points[npoints], 830 heurdata->maxiter, heurdata->minimprfac, heurdata->minimpriter, &feasibilities[npoints], nlrowgradcosts, 831 &gradcosts) ); 832 833 gradlimit -= gradcosts; 834 SCIPdebugMsg(scip, "improve point %d / %d gradlimit = %g\n", npoints, nrndpoints, gradlimit); 835 } 836 assert(npoints >= 0 && npoints <= nrndpoints); 837 838 if( npoints == 0 ) 839 goto TERMINATE; 840 841 /* 842 * 3. filter and cluster points 843 */ 844 SCIP_CALL( filterPoints(scip, points, feasibilities, npoints, &nusefulpoints) ); 845 assert(nusefulpoints >= 0); 846 SCIPdebugMsg(scip, "nusefulpoints = %d\n", nusefulpoints); 847 848 if( nusefulpoints == 0 ) 849 goto TERMINATE; 850 851 SCIP_CALL( clusterPointsGreedy(scip, points, nusefulpoints, clusteridx, &ncluster, heurdata->maxboundsize, 852 heurdata->maxreldist, heurdata->maxncluster) ); 853 assert(ncluster >= 0 && ncluster <= heurdata->maxncluster); 854 SCIPdebugMsg(scip, "ncluster = %d\n", ncluster); 855 856 SCIPsortIntPtr(clusteridx, (void**)points, nusefulpoints); 857 858 /* 859 * 4. compute start point for each cluster and use it in the sub-NLP heuristic (@ref heur_subnlp.h) 860 */ 861 start = 0; 862 while( start < nusefulpoints && clusteridx[start] != INT_MAX && !SCIPisStopped(scip) ) 863 { 864 SCIP_Bool success; 865 int end; 866 867 end = start; 868 while( end < nusefulpoints && clusteridx[start] == clusteridx[end] ) 869 ++end; 870 871 assert(end - start > 0); 872 873 /* call sub-NLP heuristic */ 874 SCIP_CALL( solveNLP(scip, heur, heurdata->heursubnlp, &points[start], end - start, &success) ); 875 SCIPdebugMsg(scip, "solveNLP result = %u\n", success); 876 877 if( success ) 878 *result = SCIP_FOUNDSOL; 879 880 /* go to the next cluster */ 881 start = end; 882 } 883 884 TERMINATE: 885 /* free memory */ 886 for( i = nrndpoints - 1; i >= 0 ; --i ) 887 { 888 assert(points[i] != NULL); 889 SCIP_CALL( SCIPfreeSol(scip, &points[i]) ); 890 } 891 892 SCIPhashmapFree(&varindex); 893 SCIPfreeBufferArray(scip, &clusteridx); 894 SCIPfreeBufferArray(scip, &feasibilities); 895 SCIPfreeBufferArray(scip, &nlrowgradcosts); 896 SCIPfreeBufferArray(scip, &points); 897 898 return SCIP_OKAY; 899 } 900 901 /* 902 * Callback methods of primal heuristic 903 */ 904 905 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 906 static 907 SCIP_DECL_HEURCOPY(heurCopyMultistart) 908 { /*lint --e{715}*/ 909 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 910 911 /* call inclusion method of primal heuristic */ 912 SCIP_CALL( SCIPincludeHeurMultistart(scip) ); 913 914 return SCIP_OKAY; 915 } 916 917 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 918 static 919 SCIP_DECL_HEURFREE(heurFreeMultistart) 920 { /*lint --e{715}*/ 921 SCIP_HEURDATA* heurdata; 922 923 /* free heuristic data */ 924 heurdata = SCIPheurGetData(heur); 925 926 SCIPfreeBlockMemory(scip, &heurdata); 927 SCIPheurSetData(heur, NULL); 928 929 return SCIP_OKAY; 930 } 931 932 /** initialization method of primal heuristic (called after problem was transformed) */ 933 static 934 SCIP_DECL_HEURINIT(heurInitMultistart) 935 { /*lint --e{715}*/ 936 SCIP_HEURDATA* heurdata; 937 938 assert( heur != NULL ); 939 940 heurdata = SCIPheurGetData(heur); 941 assert(heurdata != NULL); 942 943 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, 944 DEFAULT_RANDSEED, TRUE) ); 945 946 /* try to find sub-NLP heuristic */ 947 heurdata->heursubnlp = SCIPfindHeur(scip, "subnlp"); 948 949 return SCIP_OKAY; 950 } 951 952 /** deinitialization method of primal heuristic (called before transformed problem is freed) */ 953 static 954 SCIP_DECL_HEUREXIT(heurExitMultistart) 955 { /*lint --e{715}*/ 956 SCIP_HEURDATA* heurdata; 957 958 assert( heur != NULL ); 959 960 heurdata = SCIPheurGetData(heur); 961 assert(heurdata != NULL); 962 assert(heurdata->randnumgen != NULL); 963 964 SCIPfreeRandom(scip, &heurdata->randnumgen); 965 966 return SCIP_OKAY; 967 } 968 969 /** execution method of primal heuristic */ 970 static 971 SCIP_DECL_HEUREXEC(heurExecMultistart) 972 { /*lint --e{715}*/ 973 SCIP_HEURDATA* heurdata; 974 975 assert( heur != NULL ); 976 977 heurdata = SCIPheurGetData(heur); 978 assert(heurdata != NULL); 979 980 *result = SCIP_DIDNOTRUN; 981 982 /* check cases for which the heuristic is not applicable */ 983 if( !SCIPisNLPConstructed(scip) || heurdata->heursubnlp == NULL || SCIPgetNNlpis(scip) <= 0 ) 984 return SCIP_OKAY; 985 986 /* check whether the heuristic should be applied for a problem containing integer variables */ 987 if( heurdata->onlynlps && (SCIPgetNBinVars(scip) > 0 || SCIPgetNIntVars(scip) > 0) ) 988 return SCIP_OKAY; 989 990 *result = SCIP_DIDNOTFIND; 991 992 SCIP_CALL( applyHeur(scip, heur, heurdata, result) ); 993 994 return SCIP_OKAY; 995 } 996 997 /* 998 * primal heuristic specific interface methods 999 */ 1000 1001 /** creates the multistart primal heuristic and includes it in SCIP */ 1002 SCIP_RETCODE SCIPincludeHeurMultistart( 1003 SCIP* scip /**< SCIP data structure */ 1004 ) 1005 { 1006 SCIP_HEURDATA* heurdata; 1007 SCIP_HEUR* heur; 1008 1009 /* create multistart primal heuristic data */ 1010 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 1011 BMSclearMemory(heurdata); 1012 1013 /* include primal heuristic */ 1014 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 1015 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 1016 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecMultistart, heurdata) ); 1017 1018 assert(heur != NULL); 1019 1020 /* set non fundamental callbacks via setter functions */ 1021 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyMultistart) ); 1022 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeMultistart) ); 1023 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitMultistart) ); 1024 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitMultistart) ); 1025 1026 /* add multistart primal heuristic parameters */ 1027 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/nrndpoints", 1028 "number of random points generated per execution call", 1029 &heurdata->nrndpoints, FALSE, DEFAULT_NRNDPOINTS, 0, INT_MAX, NULL, NULL) ); 1030 1031 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxboundsize", 1032 "maximum variable domain size for unbounded variables", 1033 &heurdata->maxboundsize, FALSE, DEFAULT_MAXBOUNDSIZE, 0.0, SCIPinfinity(scip), NULL, NULL) ); 1034 1035 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiter", 1036 "number of iterations to reduce the maximum violation of a point", 1037 &heurdata->maxiter, FALSE, DEFAULT_MAXITER, 0, INT_MAX, NULL, NULL) ); 1038 1039 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprfac", 1040 "minimum required improving factor to proceed in improvement of a single point", 1041 &heurdata->minimprfac, FALSE, DEFAULT_MINIMPRFAC, -SCIPinfinity(scip), SCIPinfinity(scip), NULL, NULL) ); 1042 1043 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/minimpriter", 1044 "number of iteration when checking the minimum improvement", 1045 &heurdata->minimpriter, FALSE, DEFAULT_MINIMPRITER, 1, INT_MAX, NULL, NULL) ); 1046 1047 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxreldist", 1048 "maximum distance between two points in the same cluster", 1049 &heurdata->maxreldist, FALSE, DEFAULT_MAXRELDIST, 0.0, SCIPinfinity(scip), NULL, NULL) ); 1050 1051 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/gradlimit", 1052 "limit for gradient computations for all improvePoint() calls (0 for no limit)", 1053 &heurdata->gradlimit, FALSE, DEFAULT_GRADLIMIT, 0.0, SCIPinfinity(scip), NULL, NULL) ); 1054 1055 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxncluster", 1056 "maximum number of considered clusters per heuristic call", 1057 &heurdata->maxncluster, FALSE, DEFAULT_MAXNCLUSTER, 0, INT_MAX, NULL, NULL) ); 1058 1059 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/onlynlps", 1060 "should the heuristic run only on continuous problems?", 1061 &heurdata->onlynlps, FALSE, DEFAULT_ONLYNLPS, NULL, NULL) ); 1062 1063 return SCIP_OKAY; 1064 } 1065