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_mpec.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief mpec primal heuristic 28 * @author Felipe Serrano 29 * @author Benjamin Mueller 30 */ 31 32 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 33 34 #include "blockmemshell/memory.h" 35 #include "scip/pub_expr.h" 36 #include "scip/expr_var.h" 37 #include "scip/expr_sum.h" 38 #include "scip/expr_pow.h" 39 #include "scip/heur_mpec.h" 40 #include "scip/heur_subnlp.h" 41 #include "scip/pub_cons.h" 42 #include "scip/pub_heur.h" 43 #include "scip/pub_message.h" 44 #include "scip/pub_misc.h" 45 #include "scip/pub_nlp.h" 46 #include "scip/pub_var.h" 47 #include "scip/scip_cons.h" 48 #include "scip/scip_general.h" 49 #include "scip/scip_heur.h" 50 #include "scip/scip_mem.h" 51 #include "scip/scip_message.h" 52 #include "scip/scip_nlp.h" 53 #include "scip/scip_nlpi.h" 54 #include "scip/scip_numerics.h" 55 #include "scip/scip_param.h" 56 #include "scip/scip_prob.h" 57 #include "scip/scip_sol.h" 58 #include "scip/scip_solvingstats.h" 59 #include "scip/scip_timing.h" 60 #include <string.h> 61 62 63 #define HEUR_NAME "mpec" 64 #define HEUR_DESC "regularization heuristic for convex and nonconvex MINLPs" 65 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_DIVING 66 #define HEUR_PRIORITY -2050000 67 #define HEUR_FREQ 50 68 #define HEUR_FREQOFS 0 69 #define HEUR_MAXDEPTH -1 70 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPNODE 71 #define HEUR_USESSUBSCIP TRUE /**< disable the heuristic in sub-SCIPs, even though it does not use any */ 72 73 #define DEFAULT_INITTHETA 0.125 /**< default initial regularization right-hand side value (< 0.25) */ 74 #define DEFAULT_SIGMA 0.5 /**< default regularization update factor (< 1) */ 75 #define DEFAULT_MAXITER 100 /**< default maximum number of iterations of the MPEC loop */ 76 #define DEFAULT_MAXNLPITER 500 /**< default maximum number of NLP iterations per solve */ 77 #define DEFAULT_MINGAPLEFT 0.05 /**< default minimum amount of gap left in order to call the heuristic */ 78 #define DEFAULT_SUBNLPTRIGGER 1e-3 /**< default maximum integrality violation before triggering a sub-NLP call */ 79 #define DEFAULT_MAXNLPCOST 1e+8 /**< default maximum cost available for solving NLPs per call of the heuristic */ 80 #define DEFAULT_MINIMPROVE 0.01 /**< default factor by which heuristic should at least improve the incumbent */ 81 #define DEFAULT_MAXNUNSUCC 10 /**< default maximum number of consecutive calls for which the heuristic did not find an improving solution */ 82 83 /* 84 * Data structures 85 */ 86 87 /** primal heuristic data */ 88 struct SCIP_HeurData 89 { 90 SCIP_NLPI* nlpi; /**< nlpi used to create the nlpi problem */ 91 SCIP_NLPIPROBLEM* nlpiprob; /**< nlpi problem representing the NLP relaxation */ 92 SCIP_HASHMAP* var2idx; /**< mapping between variables and nlpi indices */ 93 SCIP_HEUR* subnlp; /**< sub-NLP heuristic */ 94 95 SCIP_Real inittheta; /**< initial regularization right-hand side value */ 96 SCIP_Real sigma; /**< regularization update factor */ 97 SCIP_Real subnlptrigger; /**< maximum number of NLP iterations per solve */ 98 SCIP_Real maxnlpcost; /**< maximum cost available for solving NLPs per call of the heuristic */ 99 SCIP_Real minimprove; /**< factor by which heuristic should at least improve the incumbent */ 100 SCIP_Real mingapleft; /**< minimum amount of gap left in order to call the heuristic */ 101 int maxiter; /**< maximum number of iterations of the MPEC loop */ 102 int maxnlpiter; /**< maximum number of NLP iterations per solve */ 103 int nunsucc; /**< number of consecutive calls for which the heuristic did not find an 104 * improving solution */ 105 int maxnunsucc; /**< maximum number of consecutive calls for which the heuristic did not 106 * find an improving solution */ 107 }; 108 109 110 /* 111 * Local methods 112 */ 113 114 /** creates the data structure for generating the current NLP relaxation */ 115 static 116 SCIP_RETCODE createNLP( 117 SCIP* scip, /**< SCIP data structure */ 118 SCIP_HEURDATA* heurdata /**< heuristic data */ 119 ) 120 { 121 SCIP_Real cutoff = SCIPinfinity(scip); 122 123 assert(heurdata != NULL); 124 assert(heurdata->nlpi != NULL); 125 126 /* NLP has been already created */ 127 if( heurdata->nlpiprob != NULL ) 128 return SCIP_OKAY; 129 130 /* compute cutoff value to ensure minimum improvement */ 131 if( SCIPgetNSols(scip) > 0 ) 132 { 133 SCIP_Real upperbound = SCIPgetUpperbound(scip) - SCIPsumepsilon(scip); 134 135 assert( !SCIPisInfinity(scip, SCIPgetUpperbound(scip)) ); 136 137 if( !SCIPisInfinity(scip, -SCIPgetLowerbound(scip)) ) 138 { 139 cutoff = (1.0 - heurdata->minimprove) * SCIPgetUpperbound(scip) 140 + heurdata->minimprove * SCIPgetLowerbound(scip); 141 } 142 else 143 { 144 if( SCIPgetUpperbound(scip) >= 0.0 ) 145 cutoff = ( 1.0 - heurdata->minimprove ) * SCIPgetUpperbound(scip); 146 else 147 cutoff = ( 1.0 + heurdata->minimprove ) * SCIPgetUpperbound(scip); 148 } 149 cutoff = MIN(upperbound, cutoff); 150 SCIPdebugMsg(scip, "set objective limit %g in [%g,%g]\n", cutoff, SCIPgetLowerbound(scip), 151 SCIPgetUpperbound(scip)); 152 } 153 154 SCIP_CALL( SCIPhashmapCreate(&heurdata->var2idx, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 155 SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, heurdata->nlpi, &heurdata->nlpiprob, "MPEC-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip), 156 heurdata->var2idx, NULL, NULL, cutoff, TRUE, FALSE) ); 157 158 return SCIP_OKAY; 159 } 160 161 /** frees the data structures for the NLP relaxation */ 162 static 163 SCIP_RETCODE freeNLP( 164 SCIP* scip, /**< SCIP data structure */ 165 SCIP_HEURDATA* heurdata /**< heuristic data */ 166 ) 167 { 168 assert(heurdata != NULL); 169 170 /* NLP has not been created yet */ 171 if( heurdata->nlpiprob == NULL ) 172 return SCIP_OKAY; 173 174 assert(heurdata->nlpi != NULL); 175 assert(heurdata->var2idx != NULL); 176 177 SCIPhashmapFree(&heurdata->var2idx); 178 SCIP_CALL( SCIPfreeNlpiProblem(scip, heurdata->nlpi, &heurdata->nlpiprob) ); 179 180 return SCIP_OKAY; 181 } 182 183 /** adds or updates the regularization constraints to the NLP; for a given parameter theta we add for each non-fixed 184 * binary variable z the constraint z*(1-z) <= theta; if these constraint are already present we update the theta on 185 * the right-hand side 186 */ 187 static 188 SCIP_RETCODE addRegularScholtes( 189 SCIP* scip, /**< SCIP data structure */ 190 SCIP_HEURDATA* heurdata, /**< heuristic data */ 191 SCIP_VAR** binvars, /**< array containing all non-fixed binary variables */ 192 int nbinvars, /**< total number of non-fixed binary variables */ 193 SCIP_Real theta, /**< regularization parameter */ 194 SCIP_Bool update /**< should the regularization constraints be added or updated? */ 195 ) 196 { 197 int i; 198 199 assert(binvars != NULL); 200 assert(nbinvars > 0); 201 202 /* add or update regularization for each non-fixed binary variables */ 203 if( !update ) 204 { 205 SCIP_NLROW** nlrows; 206 207 SCIP_CALL( SCIPallocBufferArray(scip, &nlrows, nbinvars) ); 208 209 for( i = 0; i < nbinvars; ++i ) 210 { 211 SCIP_Real one = 1.0; 212 SCIP_Real minusone = -1.0; 213 SCIP_EXPR* varexpr; 214 SCIP_EXPR* powexpr; 215 SCIP_EXPR* sumexpr; 216 char name[SCIP_MAXSTRLEN]; 217 218 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_reg", SCIPvarGetName(binvars[i])); 219 220 /* -binvars[i]^2 */ 221 SCIP_CALL( SCIPcreateExprVar(scip, &varexpr, binvars[i], NULL, NULL) ); 222 SCIP_CALL( SCIPcreateExprPow(scip, &powexpr, varexpr, 2.0, NULL, NULL) ); 223 SCIP_CALL( SCIPcreateExprSum(scip, &sumexpr, 1, &powexpr, &minusone, 0.0, NULL, NULL) ); 224 225 /* binvars[i] - binvars[i]^2 <= theta */ 226 SCIP_CALL( SCIPcreateNlRow(scip, &nlrows[i], name, 0.0, 1, &binvars[i], &one, sumexpr, -SCIPinfinity(scip), theta, SCIP_EXPRCURV_CONCAVE) ); 227 228 SCIP_CALL( SCIPreleaseExpr(scip, &sumexpr) ); 229 SCIP_CALL( SCIPreleaseExpr(scip, &powexpr) ); 230 SCIP_CALL( SCIPreleaseExpr(scip, &varexpr) ); 231 } 232 233 SCIP_CALL( SCIPaddNlpiProblemNlRows(scip, heurdata->nlpi, heurdata->nlpiprob, heurdata->var2idx, nlrows, nbinvars) ); 234 235 for( i = nbinvars-1; i >= 0; --i ) 236 { 237 SCIP_CALL( SCIPreleaseNlRow(scip, &nlrows[i]) ); 238 } 239 240 SCIPfreeBufferArray(scip, &nlrows); 241 } 242 else 243 { 244 int startidx = SCIPgetNNLPNlRows(scip) + 1; /* the cutoff is a separate constraint */ 245 SCIP_Real* lhss; 246 SCIP_Real* rhss; 247 int* indices; 248 249 SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nbinvars) ); 250 SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nbinvars) ); 251 SCIP_CALL( SCIPallocBufferArray(scip, &indices, nbinvars) ); 252 253 for( i = 0; i < nbinvars; ++i ) 254 { 255 lhss[i] = -SCIPinfinity(scip); 256 rhss[i] = theta; 257 indices[i] = startidx + i; 258 } 259 260 SCIP_CALL( SCIPchgNlpiConsSides(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lhss, rhss) ); 261 262 SCIPfreeBufferArray(scip, &indices); 263 SCIPfreeBufferArray(scip, &rhss); 264 SCIPfreeBufferArray(scip, &lhss); 265 } 266 267 return SCIP_OKAY; 268 } 269 270 /** recursive helper function to count the number of nodes in a sub-expr */ 271 static 272 int getExprSize( 273 SCIP_EXPR* expr /**< expression */ 274 ) 275 { 276 int sum; 277 int i; 278 279 assert(expr != NULL); 280 281 sum = 0; 282 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 283 { 284 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i]; 285 sum += getExprSize(child); 286 } 287 return 1 + sum; 288 } 289 290 /** main execution function of the MPEC heuristic */ 291 static 292 SCIP_RETCODE heurExec( 293 SCIP* scip, /**< SCIP data structure */ 294 SCIP_HEUR* heur, /**< MPEC heuristic */ 295 SCIP_HEURDATA* heurdata, /**< heuristic data */ 296 SCIP_RESULT* result /**< pointer to store the result */ 297 ) 298 { 299 SCIP_NLPSTATISTICS nlpstatistics; 300 SCIP_NLPPARAM nlpparam = SCIP_NLPPARAM_DEFAULT(scip); /*lint !e446*/ 301 SCIP_VAR** binvars = NULL; 302 SCIP_Real* initguess = NULL; 303 SCIP_Real* ubs = NULL; 304 SCIP_Real* lbs = NULL; 305 int* indices = NULL; 306 SCIP_Real theta = heurdata->inittheta; 307 SCIP_Real nlpcostperiter = 0.0; 308 SCIP_Real nlpcostleft = heurdata->maxnlpcost; 309 SCIP_Bool reinit = TRUE; 310 SCIP_Bool fixed = FALSE; 311 SCIP_Bool subnlpcalled = FALSE; 312 int nbinvars = 0; 313 int i; 314 315 assert(heurdata->nlpiprob != NULL); 316 assert(heurdata->var2idx != NULL); 317 assert(heurdata->nlpi != NULL); 318 assert(result != NULL); 319 320 SCIP_CALL( SCIPallocBufferArray(scip, &binvars, SCIPgetNBinVars(scip)) ); 321 322 /* collect all non-fixed binary variables */ 323 for( i = 0; i < SCIPgetNBinVars(scip); ++i ) 324 { 325 SCIP_VAR* var = SCIPgetVars(scip)[i]; 326 assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY); 327 328 if( !SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) ) 329 binvars[nbinvars++] = var; 330 } 331 332 /* all binary variables are fixed */ 333 SCIPdebugMsg(scip, "nbinvars %d\n", nbinvars); 334 if( nbinvars == 0 ) 335 goto TERMINATE; 336 337 SCIP_CALL( SCIPallocBufferArray(scip, &initguess, SCIPgetNVars(scip)) ); 338 SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nbinvars) ); 339 SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nbinvars) ); 340 SCIP_CALL( SCIPallocBufferArray(scip, &indices, nbinvars) ); 341 342 /* compute estimate cost for each NLP iteration */ 343 for( i = 0; i < SCIPgetNNLPNlRows(scip); ++i ) 344 { 345 SCIP_NLROW* nlrow = SCIPgetNLPNlRows(scip)[i]; 346 assert(nlrow != NULL); 347 348 nlpcostperiter += 1.0 * SCIPnlrowGetNLinearVars(nlrow); 349 350 if( SCIPnlrowGetExpr(nlrow) != NULL ) 351 nlpcostperiter += 3.0 * getExprSize(SCIPnlrowGetExpr(nlrow)); 352 } 353 354 /* set initial guess */ 355 for( i = 0; i < SCIPgetNVars(scip); ++i ) 356 { 357 SCIP_VAR* var = SCIPgetVars(scip)[i]; 358 initguess[i] = SCIPgetSolVal(scip, NULL, var); 359 /* SCIPdebugMsg(scip, "set initial value for %s to %g\n", SCIPvarGetName(var), initguess[i]); */ 360 } 361 SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, initguess, NULL, NULL, NULL) ); 362 363 /* set parameters of NLP solver */ 364 nlpparam.feastol /= 10.0; 365 nlpparam.opttol /= 10.0; 366 nlpparam.iterlimit = heurdata->maxnlpiter; 367 368 /* main loop */ 369 for( i = 0; i < heurdata->maxiter && *result != SCIP_FOUNDSOL && nlpcostleft > 0.0 && !SCIPisStopped(scip); ++i ) 370 { 371 SCIP_Real* primal = NULL; 372 SCIP_Bool binaryfeasible; 373 SCIP_Bool regularfeasible; 374 SCIP_NLPSOLSTAT solstat; 375 SCIP_Real maxviolbin = 0.0; 376 SCIP_Real maxviolreg = 0.0; 377 int j; 378 379 /* add or update regularization */ 380 SCIP_CALL( addRegularScholtes(scip, heurdata, binvars, nbinvars, theta, i > 0) ); 381 382 /* solve NLP */ 383 SCIP_CALL( SCIPsolveNlpiParam(scip, heurdata->nlpi, heurdata->nlpiprob, nlpparam) ); 384 solstat = SCIPgetNlpiSolstat(scip, heurdata->nlpi, heurdata->nlpiprob); 385 386 /* give up if an error occurred or no primal values are accessible */ 387 if( solstat > SCIP_NLPSOLSTAT_LOCINFEASIBLE ) 388 { 389 SCIPdebugMsg(scip, "error occured during NLP solve -> stop!\n"); 390 break; 391 } 392 393 /* update nlpcostleft */ 394 SCIP_CALL( SCIPgetNlpiStatistics(scip, heurdata->nlpi, heurdata->nlpiprob, &nlpstatistics) ); 395 nlpcostleft -= nlpstatistics.niterations * nlpcostperiter * nbinvars; 396 SCIPdebugMsg(scip, "nlpcostleft = %e\n", nlpcostleft); 397 398 SCIP_CALL( SCIPgetNlpiSolution(scip, heurdata->nlpi, heurdata->nlpiprob, &primal, NULL, NULL, NULL, NULL) ); 399 assert(primal != NULL); 400 401 /* check for binary feasibility */ 402 binaryfeasible = TRUE; 403 regularfeasible = TRUE; 404 for( j = 0; j < nbinvars; ++j ) 405 { 406 int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]); 407 binaryfeasible = binaryfeasible && SCIPisFeasIntegral(scip, primal[idx]); 408 regularfeasible = regularfeasible && SCIPisLE(scip, primal[idx] - SQR(primal[idx]), theta); 409 410 maxviolreg = MAX(maxviolreg, primal[idx] - SQR(primal[idx]) - theta); 411 maxviolbin = MAX(maxviolbin, MIN(primal[idx], 1.0-primal[idx])); 412 } 413 SCIPdebugMsg(scip, "maxviol-regularization %g maxviol-integrality %g\n", maxviolreg, maxviolbin); 414 415 /* call sub-NLP heuristic when the maximum binary infeasibility is small enough (or this is the last iteration 416 * because we reached the nlpcost limit) 417 */ 418 if( !subnlpcalled && heurdata->subnlp != NULL 419 && (SCIPisLE(scip, maxviolbin, heurdata->subnlptrigger) || nlpcostleft <= 0.0) 420 && !SCIPisStopped(scip) ) 421 { 422 SCIP_SOL* refpoint; 423 SCIP_RESULT subnlpresult; 424 425 SCIPdebugMsg(scip, "call sub-NLP heuristic because binary infeasibility is small enough\n"); 426 SCIP_CALL( SCIPcreateSol(scip, &refpoint, heur) ); 427 428 for( j = 0; j < SCIPgetNVars(scip); ++j ) 429 { 430 SCIP_VAR* var = SCIPgetVars(scip)[j]; 431 SCIP_Real val = SCIPvarIsBinary(var) ? SCIPfeasRound(scip, primal[j]) : primal[j]; 432 SCIP_CALL( SCIPsetSolVal(scip, refpoint, var, val) ); 433 } 434 435 SCIP_CALL( SCIPapplyHeurSubNlp(scip, heurdata->subnlp, &subnlpresult, refpoint, NULL) ); 436 SCIP_CALL( SCIPfreeSol(scip, &refpoint) ); 437 SCIPdebugMsg(scip, "result of sub-NLP call: %d\n", subnlpresult); 438 439 /* stop MPEC heuristic when the sub-NLP heuristic has found a feasible solution */ 440 if( subnlpresult == SCIP_FOUNDSOL ) 441 { 442 SCIPdebugMsg(scip, "sub-NLP found a feasible solution -> stop!\n"); 443 break; 444 } 445 446 subnlpcalled = TRUE; 447 } 448 449 /* NLP feasible + binary feasible -> add solution and stop */ 450 if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE && binaryfeasible ) 451 { 452 SCIP_SOL* sol; 453 SCIP_Bool stored; 454 455 SCIP_CALL( SCIPcreateSol(scip, &sol, heur) ); 456 457 for( j = 0; j < SCIPgetNVars(scip); ++j ) 458 { 459 SCIP_VAR* var = SCIPgetVars(scip)[j]; 460 assert(j == SCIPhashmapGetImageInt(heurdata->var2idx, (void*)var)); 461 SCIP_CALL( SCIPsetSolVal(scip, sol, var, primal[j]) ); 462 } 463 464 #ifdef SCIP_DEBUG 465 SCIP_CALL( SCIPtrySolFree(scip, &sol, TRUE, TRUE, TRUE, TRUE, FALSE, &stored) ); 466 #else 467 SCIP_CALL( SCIPtrySolFree(scip, &sol, FALSE, FALSE, TRUE, TRUE, FALSE, &stored) ); 468 #endif 469 SCIPdebugMsg(scip, "found a solution (stored = %u)\n", stored); 470 471 if( stored ) 472 *result = SCIP_FOUNDSOL; 473 break; 474 } 475 476 /* NLP feasible + binary infeasible -> reduce theta */ 477 else if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE && !binaryfeasible ) 478 { 479 BMScopyMemoryArray(initguess, primal, SCIPgetNVars(scip)); 480 SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, primal, NULL, NULL, NULL) ); 481 SCIPdebugMsg(scip, "update theta from %g -> %g\n", theta, theta*heurdata->sigma); 482 483 if( !reinit ) 484 { 485 SCIPdebugMsg(scip, "reinit fixed the infeasibility\n"); 486 reinit = TRUE; 487 } 488 489 theta *= heurdata->sigma; 490 491 /* unfix binary variables */ 492 if( fixed ) 493 { 494 SCIPdebugMsg(scip, "unfixing binary variables\n"); 495 for( j = 0; j < nbinvars; ++j ) 496 { 497 lbs[j] = 0.0; 498 ubs[j] = 1.0; 499 indices[j] = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]); 500 } 501 SCIP_CALL( SCIPchgNlpiVarBounds(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lbs, ubs) ); 502 fixed = FALSE; 503 } 504 } 505 506 /* NLP infeasible + regularization feasible -> stop (give up) */ 507 else if( solstat > SCIP_NLPSOLSTAT_FEASIBLE && regularfeasible ) 508 { 509 SCIPdebugMsg(scip, "NLP is infeasible but regularization constraints are satisfied -> stop!\n"); 510 break; 511 } 512 513 /* NLP infeasible + binary infeasible -> set initial point / fix binary variables */ 514 else 515 { 516 assert(solstat > SCIP_NLPSOLSTAT_FEASIBLE && !regularfeasible); 517 518 SCIPdebugMsg(scip, "NLP solution is not feasible for the NLP and the binary variables\n"); 519 520 /* stop if fixing did not resolve the infeasibility */ 521 if( fixed ) 522 { 523 SCIPdebugMsg(scip, "fixing variables did not resolve infeasibility -> stop!\n"); 524 break; 525 } 526 527 /* fix variables if reinit is FALSE; otherwise set another initial point */ 528 if( !reinit ) 529 { 530 int nfixedvars = 0; 531 532 /* fix binary variables */ 533 for( j = 0; j < nbinvars; ++j ) 534 { 535 int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]); 536 indices[j] = idx; 537 538 if( SCIPisFeasLE(scip, primal[idx] - SQR(primal[idx]), theta) ) 539 { 540 lbs[j] = 0.0; 541 ubs[j] = 1.0; 542 } 543 else 544 { 545 lbs[j] = primal[idx] >= 0.5 ? 0.0 : 1.0; 546 ubs[j] = primal[idx] >= 0.5 ? 0.0 : 1.0; 547 ++nfixedvars; 548 /* SCIPdebugMsg(scip, "fix binary variable %s = %g\n", SCIPvarGetName(binvars[j]), ubs[j]); */ 549 } 550 } 551 SCIPdebugMsg(scip, "fixed %d binary variables\n", nfixedvars); 552 SCIP_CALL( SCIPchgNlpiVarBounds(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lbs, ubs) ); 553 fixed = TRUE; 554 } 555 else 556 { 557 SCIPdebugMsg(scip, "update initial guess\n"); 558 559 /* set initial point */ 560 for( j = 0; j < nbinvars; ++j ) 561 { 562 int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]); 563 initguess[idx] = primal[idx] >= 0.5 ? 0.0 : 1.0; 564 /* SCIPdebugMsg(scip, "update init guess for %s to %g\n", SCIPvarGetName(binvars[j]), initguess[idx]); */ 565 } 566 SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, initguess, NULL, NULL, NULL) ); 567 reinit = FALSE; 568 } 569 } 570 } 571 572 TERMINATE: 573 SCIPfreeBufferArrayNull(scip, &indices); 574 SCIPfreeBufferArrayNull(scip, &ubs); 575 SCIPfreeBufferArrayNull(scip, &lbs); 576 SCIPfreeBufferArrayNull(scip, &initguess); 577 SCIPfreeBufferArray(scip, &binvars); 578 579 return SCIP_OKAY; 580 } 581 582 583 /* 584 * Callback methods of primal heuristic 585 */ 586 587 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 588 static 589 SCIP_DECL_HEURCOPY(heurCopyMpec) 590 { /*lint --e{715}*/ 591 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 592 593 /* call inclusion method of primal heuristic */ 594 SCIP_CALL( SCIPincludeHeurMpec(scip) ); 595 596 return SCIP_OKAY; 597 } 598 599 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 600 static 601 SCIP_DECL_HEURFREE(heurFreeMpec) 602 { /*lint --e{715}*/ 603 SCIP_HEURDATA* heurdata = SCIPheurGetData(heur); 604 assert(heurdata != NULL); 605 606 SCIPfreeBlockMemory(scip, &heurdata); 607 SCIPheurSetData(heur, NULL); 608 609 return SCIP_OKAY; 610 } 611 612 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */ 613 static 614 SCIP_DECL_HEURINITSOL(heurInitsolMpec) 615 { /*lint --e{715}*/ 616 SCIP_HEURDATA* heurdata = SCIPheurGetData(heur); 617 618 assert(heurdata != NULL); 619 assert(heurdata->nlpi == NULL); 620 621 if( SCIPgetNNlpis(scip) > 0 ) 622 { 623 heurdata->nlpi = SCIPgetNlpis(scip)[0]; 624 heurdata->subnlp = SCIPfindHeur(scip, "subnlp"); 625 } 626 627 return SCIP_OKAY; 628 } 629 630 631 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */ 632 static 633 SCIP_DECL_HEUREXITSOL(heurExitsolMpec) 634 { /*lint --e{715}*/ 635 SCIP_HEURDATA* heurdata = SCIPheurGetData(heur); 636 637 assert(heurdata != NULL); 638 heurdata->nlpi = NULL; 639 640 return SCIP_OKAY; 641 } 642 643 /** execution method of primal heuristic */ 644 static 645 SCIP_DECL_HEUREXEC(heurExecMpec) 646 { /*lint --e{715}*/ 647 SCIP_HEURDATA* heurdata = SCIPheurGetData(heur); 648 SCIP_CONSHDLR* sosonehdlr = SCIPfindConshdlr(scip, "SOS1"); 649 SCIP_CONSHDLR* sostwohdlr = SCIPfindConshdlr(scip, "SOS2"); 650 651 assert(heurdata != NULL); 652 653 *result = SCIP_DIDNOTRUN; 654 655 if( SCIPgetNIntVars(scip) > 0 || SCIPgetNBinVars(scip) == 0 656 || heurdata->nlpi == NULL || !SCIPisNLPConstructed(scip) 657 || heurdata->mingapleft > SCIPgetGap(scip) 658 || heurdata->nunsucc > heurdata->maxnunsucc ) 659 return SCIP_OKAY; 660 661 /* skip heuristic if constraints without a nonlinear representation are present */ 662 if( (sosonehdlr != NULL && SCIPconshdlrGetNConss(sosonehdlr) > 0) || 663 (sostwohdlr != NULL && SCIPconshdlrGetNConss(sostwohdlr) > 0) ) 664 { 665 return SCIP_OKAY; 666 } 667 668 *result = SCIP_DIDNOTFIND; 669 670 /* call MPEC method */ 671 SCIP_CALL( createNLP(scip, heurdata) ); 672 SCIP_CALL( heurExec(scip, heur, heurdata, result) ); 673 SCIP_CALL( freeNLP(scip, heurdata) ); 674 675 /* update number of unsuccessful calls */ 676 heurdata->nunsucc = (*result == SCIP_FOUNDSOL) ? 0 : heurdata->nunsucc + 1; 677 678 return SCIP_OKAY; 679 } 680 681 682 /* 683 * primal heuristic specific interface methods 684 */ 685 686 /** creates the mpec primal heuristic and includes it in SCIP */ 687 SCIP_RETCODE SCIPincludeHeurMpec( 688 SCIP* scip /**< SCIP data structure */ 689 ) 690 { 691 SCIP_HEURDATA* heurdata = NULL; 692 SCIP_HEUR* heur = NULL; 693 694 /* create mpec primal heuristic data */ 695 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 696 BMSclearMemory(heurdata); 697 698 /* include primal heuristic */ 699 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 700 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 701 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecMpec, heurdata) ); 702 703 assert(heur != NULL); 704 705 /* set non fundamental callbacks via setter functions */ 706 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyMpec) ); 707 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeMpec) ); 708 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolMpec) ); 709 SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolMpec) ); 710 711 /* add mpec primal heuristic parameters */ 712 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/inittheta", 713 "initial regularization right-hand side value", 714 &heurdata->inittheta, FALSE, DEFAULT_INITTHETA, 0.0, 0.25, NULL, NULL) ); 715 716 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/sigma", 717 "regularization update factor", 718 &heurdata->sigma, FALSE, DEFAULT_SIGMA, 0.0, 1.0, NULL, NULL) ); 719 720 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/subnlptrigger", 721 "maximum number of NLP iterations per solve", 722 &heurdata->subnlptrigger, FALSE, DEFAULT_SUBNLPTRIGGER, 0.0, 1.0, NULL, NULL) ); 723 724 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxnlpcost", 725 "maximum cost available for solving NLPs per call of the heuristic", 726 &heurdata->maxnlpcost, FALSE, DEFAULT_MAXNLPCOST, 0.0, SCIPinfinity(scip), NULL, NULL) ); 727 728 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprove", 729 "factor by which heuristic should at least improve the incumbent", 730 &heurdata->minimprove, FALSE, DEFAULT_MINIMPROVE, 0.0, 1.0, NULL, NULL) ); 731 732 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/mingapleft", 733 "minimum amount of gap left in order to call the heuristic", 734 &heurdata->mingapleft, FALSE, DEFAULT_MINGAPLEFT, 0.0, SCIPinfinity(scip), NULL, NULL) ); 735 736 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiter", 737 "maximum number of iterations of the MPEC loop", 738 &heurdata->maxiter, FALSE, DEFAULT_MAXITER, 0, INT_MAX, NULL, NULL) ); 739 740 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxnlpiter", 741 "maximum number of NLP iterations per solve", 742 &heurdata->maxnlpiter, FALSE, DEFAULT_MAXNLPITER, 0, INT_MAX, NULL, NULL) ); 743 744 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxnunsucc", 745 "maximum number of consecutive calls for which the heuristic did not find an improving solution", 746 &heurdata->maxnunsucc, FALSE, DEFAULT_MAXNUNSUCC, 0, INT_MAX, NULL, NULL) ); 747 748 return SCIP_OKAY; 749 } 750