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 benderscut_opt.c 26 * @ingroup OTHER_CFILES 27 * @brief Generates a standard Benders' decomposition optimality cut 28 * @author Stephen J. Maher 29 */ 30 31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 32 33 #include "scip/pub_expr.h" 34 #include "scip/benderscut_opt.h" 35 #include "scip/cons_linear.h" 36 #include "scip/pub_benderscut.h" 37 #include "scip/pub_benders.h" 38 #include "scip/pub_lp.h" 39 #include "scip/pub_nlp.h" 40 #include "scip/pub_message.h" 41 #include "scip/pub_misc.h" 42 #include "scip/pub_misc_linear.h" 43 #include "scip/pub_var.h" 44 #include "scip/scip.h" 45 #include <string.h> 46 47 #define BENDERSCUT_NAME "optimality" 48 #define BENDERSCUT_DESC "Standard Benders' decomposition optimality cut" 49 #define BENDERSCUT_PRIORITY 5000 50 #define BENDERSCUT_LPCUT TRUE 51 52 #define SCIP_DEFAULT_ADDCUTS FALSE /** Should cuts be generated, instead of constraints */ 53 #define SCIP_DEFAULT_CALCMIR TRUE /** Should the mixed integer rounding procedure be used for the cut */ 54 55 /* 56 * Data structures 57 */ 58 59 /** Benders' decomposition cuts data */ 60 struct SCIP_BenderscutData 61 { 62 SCIP_Bool addcuts; /**< should cuts be generated instead of constraints */ 63 SCIP_Bool calcmir; /**< should the mixed integer rounding procedure be applied to cuts */ 64 }; 65 66 67 /* 68 * Local methods 69 */ 70 71 /** in the case of numerical troubles, the LP is resolved with solution polishing activated */ 72 static 73 SCIP_RETCODE polishSolution( 74 SCIP* subproblem, /**< the SCIP data structure */ 75 SCIP_Bool* success /**< TRUE is the resolving of the LP was successful */ 76 ) 77 { 78 int oldpolishing; 79 SCIP_Bool lperror; 80 SCIP_Bool cutoff; 81 82 assert(subproblem != NULL); 83 assert(SCIPinProbing(subproblem)); 84 85 (*success) = FALSE; 86 87 /* setting the solution polishing parameter */ 88 SCIP_CALL( SCIPgetIntParam(subproblem, "lp/solutionpolishing", &oldpolishing) ); 89 SCIP_CALL( SCIPsetIntParam(subproblem, "lp/solutionpolishing", 2) ); 90 91 /* resolving the probing LP */ 92 SCIP_CALL( SCIPsolveProbingLP(subproblem, -1, &lperror, &cutoff) ); 93 94 if( SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_OPTIMAL ) 95 (*success) = TRUE; 96 97 /* resetting the solution polishing parameter */ 98 SCIP_CALL( SCIPsetIntParam(subproblem, "lp/solutionpolishing", oldpolishing) ); 99 100 return SCIP_OKAY; 101 } 102 103 /** verifying the activity of the cut when master variables are within epsilon of their upper or lower bounds 104 * 105 * When setting up the Benders' decomposition subproblem, master variables taking values that are within epsilon 106 * greater than their upper bound or less than their lower bound are set to their upper and lower bounds respectively. 107 * As such, there can be a difference between the subproblem dual solution objective and the optimality cut activity, 108 * when computed using the master problem solution directly. This check is to verify whether this difference is an 109 * actual error or due to the violation of the upper and lower bounds when setting up the Benders' decomposition 110 * subproblem. 111 */ 112 static 113 SCIP_RETCODE checkSetupTolerances( 114 SCIP* masterprob, /**< the SCIP data structure */ 115 SCIP_SOL* sol, /**< the master problem solution */ 116 SCIP_VAR** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 117 SCIP_Real* vals, /**< pointer to array of coefficients of the variables in the generated cut */ 118 SCIP_Real lhs, /**< the left hand side of the cut */ 119 SCIP_Real checkobj, /**< the objective of the subproblem computed from the dual solution */ 120 int nvars, /**< the number of variables in the cut */ 121 SCIP_Bool* valid /**< returns true is the cut is valid */ 122 ) 123 { 124 SCIP_Real verifyobj; 125 int i; 126 127 assert(masterprob != NULL); 128 assert(vars != NULL); 129 assert(vals != NULL); 130 131 /* initialising the verify objective with the left hand side of the optimality cut */ 132 verifyobj = lhs; 133 134 /* computing the activity of the cut from the master solution and the constraint values */ 135 for( i = 0; i < nvars; i++ ) 136 { 137 SCIP_Real solval; 138 139 solval = SCIPgetSolVal(masterprob, sol, vars[i]); 140 141 /* checking whether the solution value is less than or greater than the variable bounds */ 142 if( !SCIPisLT(masterprob, solval, SCIPvarGetUbLocal(vars[i])) ) 143 solval = SCIPvarGetUbLocal(vars[i]); 144 else if( !SCIPisGT(masterprob, solval, SCIPvarGetLbLocal(vars[i])) ) 145 solval = SCIPvarGetLbLocal(vars[i]); 146 147 verifyobj -= solval*vals[i]; 148 } 149 150 (*valid) = SCIPisFeasEQ(masterprob, checkobj, verifyobj); 151 152 return SCIP_OKAY; 153 } 154 155 /** when solving NLP subproblems, numerical issues are addressed by tightening the feasibility tolerance */ 156 static 157 SCIP_RETCODE resolveNLPWithTighterFeastol( 158 SCIP* subproblem, /**< the SCIP data structure */ 159 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 160 SCIP_Real multiplier, /**< the amount by which to decrease the tolerance */ 161 SCIP_Bool* success /**< TRUE is the resolving of the LP was successful */ 162 ) 163 { 164 SCIP_NLPSOLSTAT nlpsolstat; 165 #ifdef SCIP_DEBUG 166 SCIP_NLPTERMSTAT nlptermstat; 167 #endif 168 SCIP_NLPPARAM nlpparam = SCIPbendersGetNLPParam(benders); 169 #ifdef SCIP_MOREDEBUG 170 SCIP_SOL* nlpsol; 171 #endif 172 173 assert(subproblem != NULL); 174 assert(SCIPinProbing(subproblem)); 175 176 (*success) = FALSE; 177 178 /* reduce the default feasibility and optimality tolerance by given factor (typically 0.01) */ 179 nlpparam.feastol *= multiplier; 180 nlpparam.opttol *= multiplier; 181 182 SCIP_CALL( SCIPsolveNLPParam(subproblem, nlpparam) ); 183 184 nlpsolstat = SCIPgetNLPSolstat(subproblem); 185 #ifdef SCIP_DEBUG 186 nlptermstat = SCIPgetNLPTermstat(subproblem); 187 SCIPdebugMsg(subproblem, "NLP solstat %d termstat %d\n", nlpsolstat, nlptermstat); 188 #endif 189 190 if( nlpsolstat == SCIP_NLPSOLSTAT_LOCOPT || nlpsolstat == SCIP_NLPSOLSTAT_GLOBOPT 191 || nlpsolstat == SCIP_NLPSOLSTAT_FEASIBLE ) 192 { 193 #ifdef SCIP_MOREDEBUG 194 SCIP_CALL( SCIPcreateNLPSol(subproblem, &nlpsol, NULL) ); 195 SCIP_CALL( SCIPprintSol(subproblem, nlpsol, NULL, FALSE) ); 196 SCIP_CALL( SCIPfreeSol(subproblem, &nlpsol) ); 197 #endif 198 199 (*success) = TRUE; 200 } 201 202 return SCIP_OKAY; 203 } 204 205 /** adds a variable and value to the constraint/row arrays */ 206 static 207 SCIP_RETCODE addVariableToArray( 208 SCIP* masterprob, /**< the SCIP instance of the master problem */ 209 SCIP_VAR*** vars, /**< pointer to the array of variables in the generated cut with non-zero coefficient */ 210 SCIP_Real** vals, /**< pointer to the array of coefficients of the variables in the generated cut */ 211 SCIP_VAR* addvar, /**< the variable that will be added to the array */ 212 SCIP_Real addval, /**< the value that will be added to the array */ 213 int* nvars, /**< the number of variables in the variable array */ 214 int* varssize /**< the length of the variable size */ 215 ) 216 { 217 assert(masterprob != NULL); 218 assert(vars != NULL); 219 assert(*vars != NULL); 220 assert(vals != NULL); 221 assert(*vals != NULL); 222 assert(addvar != NULL); 223 assert(nvars != NULL); 224 assert(varssize != NULL); 225 226 if( *nvars >= *varssize ) 227 { 228 *varssize = SCIPcalcMemGrowSize(masterprob, *varssize + 1); 229 SCIP_CALL( SCIPreallocBufferArray(masterprob, vars, *varssize) ); 230 SCIP_CALL( SCIPreallocBufferArray(masterprob, vals, *varssize) ); 231 } 232 assert(*nvars < *varssize); 233 234 (*vars)[*nvars] = addvar; 235 (*vals)[*nvars] = addval; 236 (*nvars)++; 237 238 return SCIP_OKAY; 239 } 240 241 /** returns the variable solution either from the NLP or from the primal vals array */ 242 static 243 SCIP_Real getNlpVarSol( 244 SCIP_VAR* var, /**< the variable for which the solution is requested */ 245 SCIP_Real* primalvals, /**< the primal solutions for the NLP, can be NULL */ 246 SCIP_HASHMAP* var2idx /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */ 247 ) 248 { 249 SCIP_Real varsol; 250 int idx; 251 252 assert(var != NULL); 253 assert((primalvals == NULL && var2idx == NULL) || (primalvals != NULL && var2idx != NULL)); 254 255 if( var2idx != NULL && primalvals != NULL ) 256 { 257 assert(SCIPhashmapExists(var2idx, (void*)var) ); 258 idx = SCIPhashmapGetImageInt(var2idx, (void*)var); 259 varsol = primalvals[idx]; 260 } 261 else 262 varsol = SCIPvarGetNLPSol(var); 263 264 return varsol; 265 } 266 267 /** calculates a MIR cut from the coefficients of the standard optimality cut */ 268 static 269 SCIP_RETCODE computeMIRForOptimalityCut( 270 SCIP* masterprob, /**< the SCIP instance of the master problem */ 271 SCIP_SOL* sol, /**< primal CIP solution */ 272 SCIP_VAR** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 273 SCIP_Real* vals, /**< pointer to array of coefficients of the variables in the generated cut */ 274 SCIP_Real lhs, /**< the left hand side of the cut */ 275 SCIP_Real rhs, /**< the right hand side of the cut */ 276 int nvars, /**< the number of variables in the cut */ 277 SCIP_Real* cutcoefs, /**< the coefficients of the MIR cut */ 278 int* cutinds, /**< the variable indices of the MIR cut */ 279 SCIP_Real* cutrhs, /**< the RHS of the MIR cut */ 280 int* cutnnz, /**< the number of non-zeros in the cut */ 281 SCIP_Bool* success /**< was the MIR cut successfully computed? */ 282 ) 283 { 284 SCIP_AGGRROW* aggrrow; 285 SCIP_Real* rowvals; 286 int* rowinds; 287 288 SCIP_Real cutefficacy; 289 int cutrank; 290 SCIP_Bool cutislocal; 291 292 SCIP_Bool cutsuccess; 293 294 int i; 295 296 /* creating the aggregation row. There will be only a single row in this aggregation, since it is only used to 297 * compute the MIR coefficients 298 */ 299 SCIP_CALL( SCIPaggrRowCreate(masterprob, &aggrrow) ); 300 301 /* retrieving the indices for the variables in the optimality cut. All of the values must be negated, since the 302 * aggregation row requires a RHS, where the optimality cut is computed with an LHS 303 */ 304 SCIP_CALL( SCIPallocBufferArray(masterprob, &rowvals, nvars) ); 305 SCIP_CALL( SCIPallocBufferArray(masterprob, &rowinds, nvars) ); 306 307 assert(SCIPisInfinity(masterprob, rhs)); 308 assert(!SCIPisInfinity(masterprob, lhs)); 309 for( i = 0; i < nvars; i++ ) 310 { 311 rowinds[i] = SCIPvarGetProbindex(vars[i]); 312 rowvals[i] = -vals[i]; 313 } 314 315 /* adding the optimality cut to the aggregation row */ 316 SCIP_CALL( SCIPaggrRowAddCustomCons(masterprob, aggrrow, rowinds, rowvals, nvars, -lhs, 1.0, 1, FALSE) ); 317 318 /* calculating a flow cover for the optimality cut */ 319 SCIP_CALL( SCIPcalcFlowCover(masterprob, sol, TRUE, 0.9999, FALSE, aggrrow, cutcoefs, cutrhs, cutinds, cutnnz, 320 &cutefficacy, NULL, &cutislocal, &cutsuccess) ); 321 (*success) = cutsuccess; 322 323 /* calculating the MIR coefficients for the optimality cut */ 324 SCIP_CALL( SCIPcalcMIR(masterprob, sol, TRUE, 0.9999, TRUE, FALSE, FALSE, NULL, NULL, 0.001, 0.999, 1.0, aggrrow, 325 cutcoefs, cutrhs, cutinds, cutnnz, &cutefficacy, &cutrank, &cutislocal, &cutsuccess) ); 326 (*success) = ((*success) || cutsuccess); 327 328 /* the cut is only successful if the efficacy is high enough */ 329 (*success) = (*success) && SCIPisEfficacious(masterprob, cutefficacy); 330 331 /* try to tighten the coefficients of the cut */ 332 if( (*success) ) 333 { 334 SCIP_Bool redundant; 335 int nchgcoefs; 336 337 redundant = SCIPcutsTightenCoefficients(masterprob, FALSE, cutcoefs, cutrhs, cutinds, cutnnz, &nchgcoefs); 338 339 (*success) = !redundant; 340 } 341 342 /* freeing the local memory */ 343 SCIPfreeBufferArray(masterprob, &rowinds); 344 SCIPfreeBufferArray(masterprob, &rowvals); 345 SCIPaggrRowFree(masterprob, &aggrrow); 346 347 return SCIP_OKAY; 348 } 349 350 /** computes a standard Benders' optimality cut from the dual solutions of the LP */ 351 static 352 SCIP_RETCODE computeStandardLPOptimalityCut( 353 SCIP* masterprob, /**< the SCIP instance of the master problem */ 354 SCIP* subproblem, /**< the SCIP instance of the subproblem */ 355 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 356 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 357 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 358 SCIP_Real* lhs, /**< the left hand side of the cut */ 359 SCIP_Real* rhs, /**< the right hand side of the cut */ 360 int* nvars, /**< the number of variables in the cut */ 361 int* varssize, /**< the number of variables in the array */ 362 SCIP_Real* checkobj, /**< stores the objective function computed from the dual solution */ 363 SCIP_Bool* success /**< was the cut generation successful? */ 364 ) 365 { 366 SCIP_VAR** subvars; 367 SCIP_VAR** fixedvars; 368 int nsubvars; 369 int nfixedvars; 370 SCIP_Real dualsol; 371 SCIP_Real addval; 372 int nrows; 373 int i; 374 375 (*checkobj) = 0; 376 377 assert(masterprob != NULL); 378 assert(subproblem != NULL); 379 assert(benders != NULL); 380 assert(vars != NULL); 381 assert(*vars != NULL); 382 assert(vals != NULL); 383 assert(*vals != NULL); 384 385 (*success) = FALSE; 386 387 /* looping over all LP rows and setting the coefficients of the cut */ 388 nrows = SCIPgetNLPRows(subproblem); 389 for( i = 0; i < nrows; i++ ) 390 { 391 SCIP_ROW* lprow; 392 393 lprow = SCIPgetLPRows(subproblem)[i]; 394 assert(lprow != NULL); 395 396 dualsol = SCIProwGetDualsol(lprow); 397 assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) ); 398 399 if( SCIPisZero(subproblem, dualsol) ) 400 continue; 401 402 if( dualsol > 0.0 ) 403 addval = dualsol*SCIProwGetLhs(lprow); 404 else 405 addval = dualsol*SCIProwGetRhs(lprow); 406 407 (*lhs) += addval; 408 409 /* if the bound becomes infinite, then the cut generation terminates. */ 410 if( SCIPisInfinity(masterprob, (*lhs)) || SCIPisInfinity(masterprob, -(*lhs)) 411 || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval)) 412 { 413 (*success) = FALSE; 414 SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut. lhs = %g addval = %g.\n", (*lhs), addval); 415 return SCIP_OKAY; 416 } 417 } 418 419 nsubvars = SCIPgetNVars(subproblem); 420 subvars = SCIPgetVars(subproblem); 421 nfixedvars = SCIPgetNFixedVars(subproblem); 422 fixedvars = SCIPgetFixedVars(subproblem); 423 424 /* looping over all variables to update the coefficients in the computed cut. */ 425 for( i = 0; i < nsubvars + nfixedvars; i++ ) 426 { 427 SCIP_VAR* var; 428 SCIP_VAR* mastervar; 429 SCIP_Real redcost; 430 431 if( i < nsubvars ) 432 var = subvars[i]; 433 else 434 var = fixedvars[i - nsubvars]; 435 436 /* retrieving the master problem variable for the given subproblem variable. */ 437 SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) ); 438 439 redcost = SCIPgetVarRedcost(subproblem, var); 440 441 (*checkobj) += SCIPvarGetUnchangedObj(var)*SCIPvarGetSol(var, TRUE); 442 443 /* checking whether the subproblem variable has a corresponding master variable. */ 444 if( mastervar != NULL ) 445 { 446 SCIP_Real coef; 447 448 coef = -1.0*(SCIPvarGetObj(var) + redcost); 449 450 if( !SCIPisZero(masterprob, coef) ) 451 { 452 /* adding the variable to the storage */ 453 SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) ); 454 } 455 } 456 else 457 { 458 if( !SCIPisZero(subproblem, redcost) ) 459 { 460 addval = 0; 461 462 if( SCIPisPositive(subproblem, redcost) ) 463 addval = redcost*SCIPvarGetLbLocal(var); 464 else if( SCIPisNegative(subproblem, redcost) ) 465 addval = redcost*SCIPvarGetUbLocal(var); 466 467 (*lhs) += addval; 468 469 /* if the bound becomes infinite, then the cut generation terminates. */ 470 if( SCIPisInfinity(masterprob, (*lhs)) || SCIPisInfinity(masterprob, -(*lhs)) 471 || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval)) 472 { 473 (*success) = FALSE; 474 SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut.\n"); 475 return SCIP_OKAY; 476 } 477 } 478 } 479 } 480 481 assert(SCIPisInfinity(masterprob, (*rhs))); 482 /* the rhs should be infinite. If it changes, then there is an error */ 483 if( !SCIPisInfinity(masterprob, (*rhs)) ) 484 { 485 (*success) = FALSE; 486 SCIPdebugMsg(masterprob, "RHS is not infinite. rhs = %g.\n", (*rhs)); 487 return SCIP_OKAY; 488 } 489 490 (*success) = TRUE; 491 492 return SCIP_OKAY; 493 } 494 495 /** computes a standard Benders' optimality cut from the dual solutions of the NLP */ 496 static 497 SCIP_RETCODE computeStandardNLPOptimalityCut( 498 SCIP* masterprob, /**< the SCIP instance of the master problem */ 499 SCIP* subproblem, /**< the SCIP instance of the subproblem */ 500 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 501 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 502 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 503 SCIP_Real* lhs, /**< the left hand side of the cut */ 504 SCIP_Real* rhs, /**< the right hand side of the cut */ 505 int* nvars, /**< the number of variables in the cut */ 506 int* varssize, /**< the number of variables in the array */ 507 SCIP_Real objective, /**< the objective function of the subproblem */ 508 SCIP_Real* primalvals, /**< the primal solutions for the NLP, can be NULL */ 509 SCIP_Real* consdualvals, /**< dual variables for the constraints, can be NULL */ 510 SCIP_Real* varlbdualvals, /**< the dual variables for the variable lower bounds, can be NULL */ 511 SCIP_Real* varubdualvals, /**< the dual variables for the variable upper bounds, can be NULL */ 512 SCIP_HASHMAP* row2idx, /**< mapping between the row in the subproblem to the index in the dual array, can be NULL */ 513 SCIP_HASHMAP* var2idx, /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */ 514 SCIP_Real* checkobj, /**< stores the objective function computed from the dual solution */ 515 SCIP_Bool* success /**< was the cut generation successful? */ 516 ) 517 { 518 SCIP_VAR** subvars; 519 SCIP_VAR** fixedvars; 520 int nsubvars; 521 int nfixedvars; 522 SCIP_Real dirderiv; 523 SCIP_Real dualsol; 524 int nrows; 525 int idx; 526 int i; 527 528 (*checkobj) = 0; 529 530 assert(masterprob != NULL); 531 assert(subproblem != NULL); 532 assert(benders != NULL); 533 assert(SCIPisNLPConstructed(subproblem)); 534 assert(SCIPgetNLPSolstat(subproblem) <= SCIP_NLPSOLSTAT_FEASIBLE || consdualvals != NULL); 535 assert(SCIPhasNLPSolution(subproblem) || consdualvals != NULL); 536 537 (*success) = FALSE; 538 539 if( !(primalvals == NULL && consdualvals == NULL && varlbdualvals == NULL && varubdualvals == NULL && row2idx == NULL && var2idx == NULL) 540 && !(primalvals != NULL && consdualvals != NULL && varlbdualvals != NULL && varubdualvals != NULL && row2idx != NULL && var2idx != NULL) ) /*lint !e845*/ 541 { 542 SCIPerrorMessage("The optimality cut must generated from either a SCIP instance or all of the dual solutions and indices must be supplied"); 543 (*success) = FALSE; 544 545 return SCIP_ERROR; 546 } 547 548 nsubvars = SCIPgetNNLPVars(subproblem); 549 subvars = SCIPgetNLPVars(subproblem); 550 nfixedvars = SCIPgetNFixedVars(subproblem); 551 fixedvars = SCIPgetFixedVars(subproblem); 552 553 /* our optimality cut implementation assumes that SCIP did not modify the objective function and sense, 554 * that is, that the objective function value of the NLP corresponds to the value of the auxiliary variable 555 * if that wouldn't be the case, then the scaling and offset may have to be considered when adding the 556 * auxiliary variable to the cut (cons/row)? 557 */ 558 assert(SCIPgetTransObjoffset(subproblem) == 0.0); 559 assert(SCIPgetTransObjscale(subproblem) == 1.0); 560 assert(SCIPgetObjsense(subproblem) == SCIP_OBJSENSE_MINIMIZE); 561 562 (*lhs) = objective; 563 assert(!SCIPisInfinity(subproblem, REALABS(*lhs))); 564 565 (*rhs) = SCIPinfinity(masterprob); 566 567 dirderiv = 0.0; 568 569 /* looping over all NLP rows and setting the corresponding coefficients of the cut */ 570 nrows = SCIPgetNNLPNlRows(subproblem); 571 for( i = 0; i < nrows; i++ ) 572 { 573 SCIP_NLROW* nlrow; 574 575 nlrow = SCIPgetNLPNlRows(subproblem)[i]; 576 assert(nlrow != NULL); 577 578 if( row2idx != NULL && consdualvals != NULL ) 579 { 580 assert(SCIPhashmapExists(row2idx, (void*)nlrow) ); 581 idx = SCIPhashmapGetImageInt(row2idx, (void*)nlrow); 582 dualsol = consdualvals[idx]; 583 } 584 else 585 dualsol = SCIPnlrowGetDualsol(nlrow); 586 assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) ); 587 588 if( SCIPisZero(subproblem, dualsol) ) 589 continue; 590 591 SCIP_CALL( SCIPaddNlRowGradientBenderscutOpt(masterprob, subproblem, benders, nlrow, 592 -dualsol, primalvals, var2idx, &dirderiv, vars, vals, nvars, varssize) ); 593 } 594 595 /* looping over sub- and fixed variables to compute checkobj */ 596 for( i = 0; i < nsubvars; i++ ) 597 (*checkobj) += SCIPvarGetObj(subvars[i]) * getNlpVarSol(subvars[i], primalvals, var2idx); 598 599 for( i = 0; i < nfixedvars; i++ ) 600 *checkobj += SCIPvarGetUnchangedObj(fixedvars[i]) * getNlpVarSol(fixedvars[i], primalvals, var2idx); 601 602 *lhs += dirderiv; 603 604 /* if the side became infinite or dirderiv was infinite, then the cut generation terminates. */ 605 if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs) 606 || SCIPisInfinity(masterprob, dirderiv) || SCIPisInfinity(masterprob, -dirderiv)) 607 { 608 (*success) = FALSE; 609 SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut. lhs = %g dirderiv = %g.\n", *lhs, dirderiv); 610 return SCIP_OKAY; 611 } 612 613 (*success) = TRUE; 614 615 return SCIP_OKAY; 616 } 617 618 619 /** Adds the auxiliary variable to the generated cut. If this is the first optimality cut for the subproblem, then the 620 * auxiliary variable is first created and added to the master problem. 621 */ 622 static 623 SCIP_RETCODE addAuxiliaryVariableToCut( 624 SCIP* masterprob, /**< the SCIP instance of the master problem */ 625 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 626 SCIP_VAR** vars, /**< the variables in the generated cut with non-zero coefficient */ 627 SCIP_Real* vals, /**< the coefficients of the variables in the generated cut */ 628 int* nvars, /**< the number of variables in the cut */ 629 int probnumber /**< the number of the pricing problem */ 630 ) 631 { 632 SCIP_VAR* auxiliaryvar; 633 634 assert(masterprob != NULL); 635 assert(benders != NULL); 636 assert(vars != NULL); 637 assert(vals != NULL); 638 639 auxiliaryvar = SCIPbendersGetAuxiliaryVar(benders, probnumber); 640 641 vars[(*nvars)] = auxiliaryvar; 642 vals[(*nvars)] = 1.0; 643 (*nvars)++; 644 645 return SCIP_OKAY; 646 } 647 648 649 /* 650 * Callback methods of Benders' decomposition cuts 651 */ 652 653 /** destructor of Benders' decomposition cuts to free user data (called when SCIP is exiting) */ 654 static 655 SCIP_DECL_BENDERSCUTFREE(benderscutFreeOpt) 656 { /*lint --e{715}*/ 657 SCIP_BENDERSCUTDATA* benderscutdata; 658 659 assert( benderscut != NULL ); 660 assert( strcmp(SCIPbenderscutGetName(benderscut), BENDERSCUT_NAME) == 0 ); 661 662 /* free Benders' cut data */ 663 benderscutdata = SCIPbenderscutGetData(benderscut); 664 assert( benderscutdata != NULL ); 665 666 SCIPfreeBlockMemory(scip, &benderscutdata); 667 668 SCIPbenderscutSetData(benderscut, NULL); 669 670 return SCIP_OKAY; 671 } 672 673 674 /** execution method of Benders' decomposition cuts */ 675 static 676 SCIP_DECL_BENDERSCUTEXEC(benderscutExecOpt) 677 { /*lint --e{715}*/ 678 SCIP* subproblem; 679 SCIP_BENDERSCUTDATA* benderscutdata; 680 SCIP_Bool nlprelaxation; 681 SCIP_Bool addcut; 682 char cutname[SCIP_MAXSTRLEN]; 683 684 assert(scip != NULL); 685 assert(benders != NULL); 686 assert(benderscut != NULL); 687 assert(result != NULL); 688 assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders)); 689 690 /* retrieving the Benders' cut data */ 691 benderscutdata = SCIPbenderscutGetData(benderscut); 692 693 /* if the cuts are generated prior to the solving stage, then rows can not be generated. So constraints must be 694 * added to the master problem. 695 */ 696 if( SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE ) 697 addcut = FALSE; 698 else 699 addcut = benderscutdata->addcuts; 700 701 /* setting the name of the generated cut */ 702 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "optimalitycut_%d_%" SCIP_LONGINT_FORMAT, probnumber, 703 SCIPbenderscutGetNFound(benderscut) ); 704 705 subproblem = SCIPbendersSubproblem(benders, probnumber); 706 707 if( subproblem == NULL ) 708 { 709 SCIPdebugMsg(scip, "The subproblem %d is set to NULL. The <%s> Benders' decomposition cut can not be executed.\n", 710 probnumber, BENDERSCUT_NAME); 711 712 (*result) = SCIP_DIDNOTRUN; 713 return SCIP_OKAY; 714 } 715 716 /* setting a flag to indicate whether the NLP relaxation should be used to generate cuts */ 717 nlprelaxation = SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem); 718 719 /* only generate optimality cuts if the subproblem LP or NLP is optimal, 720 * since we use the dual solution of the LP/NLP to construct the optimality cut 721 */ 722 if( SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING && 723 ((!nlprelaxation && SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_OPTIMAL) || 724 (nlprelaxation && SCIPgetNLPSolstat(subproblem) <= SCIP_NLPSOLSTAT_FEASIBLE)) ) 725 { 726 /* generating a cut for a given subproblem */ 727 SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname, 728 SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut, 729 FALSE, result) ); 730 731 /* if it was not possible to generate a cut, this could be due to numerical issues. So the solution to the LP is 732 * resolved and the generation of the cut is reattempted. For NLPs, we do not have such a polishing yet. 733 */ 734 if( (*result) == SCIP_DIDNOTFIND ) 735 { 736 SCIP_Bool success; 737 738 SCIPdebugMsg(scip, "Numerical trouble generating optimality cut for subproblem %d.\n", probnumber); 739 740 if( !nlprelaxation ) 741 { 742 SCIPdebugMsg(scip, "Attempting to polish the LP solution to find an alternative dual extreme point.\n"); 743 744 SCIP_CALL( polishSolution(subproblem, &success) ); 745 746 /* only attempt to generate a cut if the solution polishing was successful */ 747 if( success ) 748 { 749 SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname, 750 SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut, 751 FALSE, result) ); 752 } 753 } 754 else 755 { 756 SCIP_Real multiplier = 0.01; 757 758 SCIPdebugMsg(scip, "Attempting to resolve the NLP with a tighter feasibility tolerance to find an " 759 "alternative dual extreme point.\n"); 760 761 while( multiplier > 1e-06 && (*result) == SCIP_DIDNOTFIND ) 762 { 763 SCIP_CALL( resolveNLPWithTighterFeastol(subproblem, benders, multiplier, &success) ); 764 765 if( success ) 766 { 767 SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname, 768 SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut, 769 FALSE, result) ); 770 } 771 772 multiplier *= 0.1; 773 } 774 } 775 } 776 } 777 778 return SCIP_OKAY; 779 } 780 781 782 /* 783 * Benders' decomposition cuts specific interface methods 784 */ 785 786 /** creates the opt Benders' decomposition cuts and includes it in SCIP */ 787 SCIP_RETCODE SCIPincludeBenderscutOpt( 788 SCIP* scip, /**< SCIP data structure */ 789 SCIP_BENDERS* benders /**< Benders' decomposition */ 790 ) 791 { 792 SCIP_BENDERSCUTDATA* benderscutdata; 793 SCIP_BENDERSCUT* benderscut; 794 char paramname[SCIP_MAXSTRLEN]; 795 796 assert(benders != NULL); 797 798 /* create opt Benders' decomposition cuts data */ 799 SCIP_CALL( SCIPallocBlockMemory(scip, &benderscutdata) ); 800 801 benderscut = NULL; 802 803 /* include Benders' decomposition cuts */ 804 SCIP_CALL( SCIPincludeBenderscutBasic(scip, benders, &benderscut, BENDERSCUT_NAME, BENDERSCUT_DESC, 805 BENDERSCUT_PRIORITY, BENDERSCUT_LPCUT, benderscutExecOpt, benderscutdata) ); 806 807 assert(benderscut != NULL); 808 809 /* setting the non fundamental callbacks via setter functions */ 810 SCIP_CALL( SCIPsetBenderscutFree(scip, benderscut, benderscutFreeOpt) ); 811 812 /* add opt Benders' decomposition cuts parameters */ 813 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/benderscut/%s/addcuts", 814 SCIPbendersGetName(benders), BENDERSCUT_NAME); 815 SCIP_CALL( SCIPaddBoolParam(scip, paramname, 816 "should cuts be generated and added to the cutpool instead of global constraints directly added to the problem.", 817 &benderscutdata->addcuts, FALSE, SCIP_DEFAULT_ADDCUTS, NULL, NULL) ); 818 819 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/benderscut/%s/mir", 820 SCIPbendersGetName(benders), BENDERSCUT_NAME); 821 SCIP_CALL( SCIPaddBoolParam(scip, paramname, 822 "should the mixed integer rounding procedure be applied to cuts", 823 &benderscutdata->calcmir, FALSE, SCIP_DEFAULT_CALCMIR, NULL, NULL) ); 824 825 return SCIP_OKAY; 826 } 827 828 /** Generates a classical Benders' optimality cut using the dual solutions from the subproblem or the input arrays. If 829 * the dual solutions are input as arrays, then a mapping between the array indices and the rows/variables is required. 830 * As a cut strengthening approach, when an optimality cut is being generated (i.e. not for feasibility cuts) a MIR 831 * procedure is performed on the row. This procedure attempts to find a stronger constraint, if this doesn't happen, 832 * then the original constraint is added to SCIP. 833 * 834 * This method can also be used to generate a feasibility cut, if a problem to minimise the infeasibilities has been solved 835 * to generate the dual solutions 836 */ 837 SCIP_RETCODE SCIPgenerateAndApplyBendersOptCut( 838 SCIP* masterprob, /**< the SCIP instance of the master problem */ 839 SCIP* subproblem, /**< the SCIP instance of the pricing problem */ 840 SCIP_BENDERS* benders, /**< the benders' decomposition */ 841 SCIP_BENDERSCUT* benderscut, /**< the benders' decomposition cut method */ 842 SCIP_SOL* sol, /**< primal CIP solution */ 843 int probnumber, /**< the number of the pricing problem */ 844 char* cutname, /**< the name for the cut to be generated */ 845 SCIP_Real objective, /**< the objective function of the subproblem */ 846 SCIP_Real* primalvals, /**< the primal solutions for the NLP, can be NULL */ 847 SCIP_Real* consdualvals, /**< dual variables for the constraints, can be NULL */ 848 SCIP_Real* varlbdualvals, /**< the dual variables for the variable lower bounds, can be NULL */ 849 SCIP_Real* varubdualvals, /**< the dual variables for the variable upper bounds, can be NULL */ 850 SCIP_HASHMAP* row2idx, /**< mapping between the row in the subproblem to the index in the dual array, can be NULL */ 851 SCIP_HASHMAP* var2idx, /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */ 852 SCIP_BENDERSENFOTYPE type, /**< the enforcement type calling this function */ 853 SCIP_Bool addcut, /**< should the Benders' cut be added as a cut or constraint */ 854 SCIP_Bool feasibilitycut, /**< is this called for the generation of a feasibility cut */ 855 SCIP_RESULT* result /**< the result from solving the subproblems */ 856 ) 857 { 858 SCIP_CONSHDLR* consbenders; 859 SCIP_CONS* cons; 860 SCIP_ROW* row; 861 SCIP_VAR** vars; 862 SCIP_Real* vals; 863 SCIP_Real lhs; 864 SCIP_Real rhs; 865 int nvars; 866 int varssize; 867 int nmastervars; 868 SCIP_Bool calcmir; 869 SCIP_Bool optimal; 870 SCIP_Bool success; 871 SCIP_Bool mirsuccess; 872 873 SCIP_Real checkobj; 874 SCIP_Real verifyobj; 875 876 assert(masterprob != NULL); 877 assert(subproblem != NULL); 878 assert(benders != NULL); 879 assert(benderscut != NULL); 880 assert(result != NULL); 881 assert((primalvals == NULL && consdualvals == NULL && varlbdualvals == NULL && varubdualvals == NULL 882 && row2idx == NULL && var2idx == NULL) 883 || (primalvals != NULL && consdualvals != NULL && varlbdualvals != NULL && varubdualvals != NULL 884 && row2idx != NULL && var2idx != NULL)); 885 886 row = NULL; 887 cons = NULL; 888 889 calcmir = SCIPbenderscutGetData(benderscut)->calcmir && SCIPgetStage(masterprob) >= SCIP_STAGE_INITSOLVE && SCIPgetSubscipDepth(masterprob) == 0; 890 success = FALSE; 891 mirsuccess = FALSE; 892 893 /* retrieving the Benders' decomposition constraint handler */ 894 consbenders = SCIPfindConshdlr(masterprob, "benders"); 895 896 /* checking the optimality of the original problem with a comparison between the auxiliary variable and the 897 * objective value of the subproblem */ 898 if( feasibilitycut ) 899 optimal = FALSE; 900 else 901 { 902 SCIP_CALL( SCIPcheckBendersSubproblemOptimality(masterprob, benders, sol, probnumber, &optimal) ); 903 } 904 905 if( optimal ) 906 { 907 (*result) = SCIP_FEASIBLE; 908 SCIPdebugMsg(masterprob, "No cut added for subproblem %d\n", probnumber); 909 return SCIP_OKAY; 910 } 911 912 /* allocating memory for the variable and values arrays */ 913 nmastervars = SCIPgetNVars(masterprob) + SCIPgetNFixedVars(masterprob); 914 SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vars, nmastervars) ); 915 SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vals, nmastervars) ); 916 lhs = 0.0; 917 rhs = SCIPinfinity(masterprob); 918 nvars = 0; 919 varssize = nmastervars; 920 921 if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) ) 922 { 923 /* computing the coefficients of the optimality cut */ 924 SCIP_CALL( computeStandardNLPOptimalityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &rhs, &nvars, 925 &varssize, objective, primalvals, consdualvals, varlbdualvals, varubdualvals, row2idx, 926 var2idx, &checkobj, &success) ); 927 } 928 else 929 { 930 /* computing the coefficients of the optimality cut */ 931 SCIP_CALL( computeStandardLPOptimalityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &rhs, &nvars, 932 &varssize, &checkobj, &success) ); 933 } 934 935 /* if success is FALSE, then there was an error in generating the optimality cut. No cut will be added to the master 936 * problem. Otherwise, the constraint is added to the master problem. 937 */ 938 if( !success ) 939 { 940 (*result) = SCIP_DIDNOTFIND; 941 SCIPdebugMsg(masterprob, "Error in generating Benders' optimality cut for problem %d.\n", probnumber); 942 } 943 else 944 { 945 /* initially a row/constraint is created for the optimality cut using the master variables and coefficients 946 * computed in computeStandardLPOptimalityCut. At this stage, the auxiliary variable is not added since the 947 * activity of the row/constraint in its current form is used to determine the validity of the optimality cut. 948 */ 949 if( addcut ) 950 { 951 SCIP_CALL( SCIPcreateEmptyRowConshdlr(masterprob, &row, consbenders, cutname, lhs, rhs, FALSE, FALSE, TRUE) ); 952 SCIP_CALL( SCIPaddVarsToRow(masterprob, row, nvars, vars, vals) ); 953 } 954 else 955 { 956 SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cons, cutname, nvars, vars, vals, lhs, rhs) ); 957 SCIP_CALL( SCIPsetConsDynamic(masterprob, cons, TRUE) ); 958 SCIP_CALL( SCIPsetConsRemovable(masterprob, cons, TRUE) ); 959 } 960 961 /* computing the objective function from the cut activity to verify the accuracy of the constraint */ 962 verifyobj = 0.0; 963 if( addcut ) 964 { 965 verifyobj += SCIProwGetLhs(row) - SCIPgetRowSolActivity(masterprob, row, sol); 966 } 967 else 968 { 969 verifyobj += SCIPgetLhsLinear(masterprob, cons) - SCIPgetActivityLinear(masterprob, cons, sol); 970 } 971 972 if( feasibilitycut && verifyobj < SCIPfeastol(masterprob) ) 973 { 974 success = FALSE; 975 SCIPdebugMsg(masterprob, "The violation of the feasibility cut (%g) is too small. Skipping feasibility cut.\n", verifyobj); 976 } 977 978 /* it is possible that numerics will cause the generated cut to be invalid. This cut should not be added to the 979 * master problem, since its addition could cut off feasible solutions. The success flag is set of false, indicating 980 * that the Benders' cut could not find a valid cut. 981 */ 982 if( !feasibilitycut && !SCIPisFeasEQ(masterprob, checkobj, verifyobj) ) 983 { 984 SCIP_Bool valid; 985 986 /* the difference in the checkobj and verifyobj could be due to the setup tolerances. This is checked, and if 987 * so, then the generated cut is still valid 988 */ 989 SCIP_CALL( checkSetupTolerances(masterprob, sol, vars, vals, lhs, checkobj, nvars, &valid) ); 990 991 if( !valid ) 992 { 993 success = FALSE; 994 SCIPdebugMsg(masterprob, "The objective function and cut activity are not equal (%g != %g).\n", checkobj, 995 verifyobj); 996 997 #ifdef SCIP_DEBUG 998 /* we only need to abort if cut strengthen is not used. If cut strengthen has been used in this round and the 999 * cut could not be generated, then another subproblem solving round will be executed 1000 */ 1001 if( !SCIPbendersInStrengthenRound(benders) ) 1002 { 1003 #ifdef SCIP_MOREDEBUG 1004 int i; 1005 1006 for( i = 0; i < nvars; i++ ) 1007 printf("<%s> %g %g\n", SCIPvarGetName(vars[i]), vals[i], SCIPgetSolVal(masterprob, sol, vars[i])); 1008 #endif 1009 SCIPABORT(); 1010 } 1011 #endif 1012 } 1013 } 1014 1015 if( success ) 1016 { 1017 /* adding the auxiliary variable to the optimality cut. The auxiliary variable is added to the vars and vals 1018 * arrays prior to the execution of the MIR procedure. This is necessary because the MIR procedure must be 1019 * executed on the complete cut, not just the row/constraint without the auxiliary variable. 1020 */ 1021 if( !feasibilitycut ) 1022 { 1023 SCIP_CALL( addAuxiliaryVariableToCut(masterprob, benders, vars, vals, &nvars, probnumber) ); 1024 } 1025 1026 /* performing the MIR procedure. If the procedure is successful, then the vars and vals arrays are no longer 1027 * needed for creating the optimality cut. These are superseeded with the cutcoefs and cutinds arrays. In the 1028 * case that the MIR procedure is successful, the row/constraint that has been created previously is destroyed 1029 * and the MIR cut is added in its place 1030 */ 1031 if( calcmir ) 1032 { 1033 SCIP_Real* cutcoefs; 1034 int* cutinds; 1035 SCIP_Real cutrhs; 1036 int cutnnz; 1037 1038 /* allocating memory to compute the MIR cut */ 1039 SCIP_CALL( SCIPallocBufferArray(masterprob, &cutcoefs, nvars) ); 1040 SCIP_CALL( SCIPallocBufferArray(masterprob, &cutinds, nvars) ); 1041 1042 SCIP_CALL( computeMIRForOptimalityCut(masterprob, sol, vars, vals, lhs, rhs, nvars, cutcoefs, 1043 cutinds, &cutrhs, &cutnnz, &mirsuccess) ); 1044 1045 /* if the MIR cut was computed successfully, then the current row/constraint needs to be destroyed and 1046 * replaced with the updated coefficients 1047 */ 1048 if( mirsuccess ) 1049 { 1050 SCIP_VAR** mastervars; 1051 int i; 1052 1053 mastervars = SCIPgetVars(masterprob); 1054 1055 if( addcut ) 1056 { 1057 SCIP_CALL( SCIPreleaseRow(masterprob, &row) ); 1058 1059 SCIP_CALL( SCIPcreateEmptyRowConshdlr(masterprob, &row, consbenders, cutname, 1060 -SCIPinfinity(masterprob), cutrhs, FALSE, FALSE, TRUE) ); 1061 1062 for( i = 0; i < cutnnz; i++) 1063 { 1064 SCIP_CALL( SCIPaddVarToRow(masterprob, row, mastervars[cutinds[i]], cutcoefs[i]) ); 1065 } 1066 } 1067 else 1068 { 1069 SCIP_CALL( SCIPreleaseCons(masterprob, &cons) ); 1070 1071 SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cons, cutname, 0, NULL, NULL, 1072 -SCIPinfinity(masterprob), cutrhs) ); 1073 SCIP_CALL( SCIPsetConsDynamic(masterprob, cons, TRUE) ); 1074 SCIP_CALL( SCIPsetConsRemovable(masterprob, cons, TRUE) ); 1075 1076 for( i = 0; i < cutnnz; i++ ) 1077 { 1078 SCIP_CALL( SCIPaddCoefLinear(masterprob, cons, mastervars[cutinds[i]], cutcoefs[i]) ); 1079 } 1080 } 1081 } 1082 1083 /* freeing the memory required to compute the MIR cut */ 1084 SCIPfreeBufferArray(masterprob, &cutinds); 1085 SCIPfreeBufferArray(masterprob, &cutcoefs); 1086 } 1087 1088 /* adding the constraint to the master problem */ 1089 if( addcut ) 1090 { 1091 SCIP_Bool infeasible; 1092 1093 /* adding the auxiliary variable coefficient to the row. This is only added if the MIR procedure is not 1094 * successful. If the MIR procedure was successful, then the auxiliary variable is already included in the 1095 * row 1096 */ 1097 if( !feasibilitycut && !mirsuccess ) 1098 { 1099 SCIP_CALL( SCIPaddVarToRow(masterprob, row, vars[nvars - 1], vals[nvars - 1]) ); 1100 } 1101 1102 if( type == SCIP_BENDERSENFOTYPE_LP || type == SCIP_BENDERSENFOTYPE_RELAX ) 1103 { 1104 SCIP_CALL( SCIPaddRow(masterprob, row, FALSE, &infeasible) ); 1105 assert(!infeasible); 1106 } 1107 else 1108 { 1109 assert(type == SCIP_BENDERSENFOTYPE_CHECK || type == SCIP_BENDERSENFOTYPE_PSEUDO); 1110 SCIP_CALL( SCIPaddPoolCut(masterprob, row) ); 1111 } 1112 1113 (*result) = SCIP_SEPARATED; 1114 } 1115 else 1116 { 1117 /* adding the auxiliary variable coefficient to the row. This is only added if the MIR procedure is not 1118 * successful. If the MIR procedure was successful, then the auxiliary variable is already included in the 1119 * constraint. 1120 */ 1121 if( !feasibilitycut && !mirsuccess ) 1122 { 1123 SCIP_CALL( SCIPaddCoefLinear(masterprob, cons, vars[nvars - 1], vals[nvars - 1]) ); 1124 } 1125 1126 SCIPdebugPrintCons(masterprob, cons, NULL); 1127 1128 SCIP_CALL( SCIPaddCons(masterprob, cons) ); 1129 1130 (*result) = SCIP_CONSADDED; 1131 } 1132 1133 /* storing the data that is used to create the cut */ 1134 SCIP_CALL( SCIPstoreBendersCut(masterprob, benders, vars, vals, lhs, rhs, nvars) ); 1135 } 1136 else 1137 { 1138 (*result) = SCIP_DIDNOTFIND; 1139 SCIPdebugMsg(masterprob, "Error in generating Benders' %s cut for problem %d.\n", feasibilitycut ? "feasibility" : "optimality", probnumber); 1140 } 1141 1142 /* releasing the row or constraint */ 1143 if( addcut ) 1144 { 1145 /* release the row */ 1146 SCIP_CALL( SCIPreleaseRow(masterprob, &row) ); 1147 } 1148 else 1149 { 1150 /* release the constraint */ 1151 SCIP_CALL( SCIPreleaseCons(masterprob, &cons) ); 1152 } 1153 } 1154 1155 SCIPfreeBufferArray(masterprob, &vals); 1156 SCIPfreeBufferArray(masterprob, &vars); 1157 1158 return SCIP_OKAY; 1159 } 1160 1161 1162 /** adds the gradient of a nonlinear row in the current NLP solution of a subproblem to a linear row or constraint in the master problem 1163 * 1164 * Only computes gradient w.r.t. master problem variables. 1165 * Computes also the directional derivative, that is, mult times gradient times solution. 1166 */ 1167 SCIP_RETCODE SCIPaddNlRowGradientBenderscutOpt( 1168 SCIP* masterprob, /**< the SCIP instance of the master problem */ 1169 SCIP* subproblem, /**< the SCIP instance of the subproblem */ 1170 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 1171 SCIP_NLROW* nlrow, /**< nonlinear row */ 1172 SCIP_Real mult, /**< multiplier */ 1173 SCIP_Real* primalvals, /**< the primal solutions for the NLP, can be NULL */ 1174 SCIP_HASHMAP* var2idx, /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */ 1175 SCIP_Real* dirderiv, /**< storage to add directional derivative */ 1176 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 1177 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 1178 int* nvars, /**< the number of variables in the cut */ 1179 int* varssize /**< the number of variables in the array */ 1180 ) 1181 { 1182 SCIP_EXPR* expr; 1183 SCIP_VAR* var; 1184 SCIP_VAR* mastervar; 1185 SCIP_Real coef; 1186 int i; 1187 1188 assert(masterprob != NULL); 1189 assert(subproblem != NULL); 1190 assert(benders != NULL); 1191 assert(nlrow != NULL); 1192 assert((primalvals == NULL && var2idx == NULL) || (primalvals != NULL && var2idx != NULL)); 1193 assert(mult != 0.0); 1194 assert(dirderiv != NULL); 1195 assert(vars != NULL); 1196 assert(vals != NULL); 1197 1198 /* linear part */ 1199 for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ ) 1200 { 1201 var = SCIPnlrowGetLinearVars(nlrow)[i]; 1202 assert(var != NULL); 1203 1204 /* retrieving the master problem variable for the given subproblem variable. */ 1205 SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) ); 1206 if( mastervar == NULL ) 1207 continue; 1208 1209 coef = mult * SCIPnlrowGetLinearCoefs(nlrow)[i]; 1210 1211 /* adding the variable to the storage */ 1212 SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) ); 1213 1214 *dirderiv += coef * getNlpVarSol(var, primalvals, var2idx); 1215 } 1216 1217 /* expression part */ 1218 expr = SCIPnlrowGetExpr(nlrow); 1219 if( expr != NULL ) 1220 { 1221 SCIP_SOL* primalsol; 1222 SCIP_EXPRITER* it; 1223 1224 /* create primalsol, either from primalvals, or pointing to NLP solution */ 1225 if( primalvals != NULL ) 1226 { 1227 SCIP_CALL( SCIPcreateSol(subproblem, &primalsol, NULL) ); 1228 1229 /* TODO would be better to change primalvals to a SCIP_SOL and do this once for the whole NLP instead of repeating it for each expr */ 1230 for( i = 0; i < SCIPhashmapGetNEntries(var2idx); ++i ) 1231 { 1232 SCIP_HASHMAPENTRY* entry; 1233 entry = SCIPhashmapGetEntry(var2idx, i); 1234 if( entry == NULL ) 1235 continue; 1236 SCIP_CALL( SCIPsetSolVal(subproblem, primalsol, (SCIP_VAR*) SCIPhashmapEntryGetOrigin(entry), primalvals[SCIPhashmapEntryGetImageInt(entry)]) ); 1237 } 1238 } 1239 else 1240 { 1241 SCIP_CALL( SCIPcreateNLPSol(subproblem, &primalsol, NULL) ); 1242 } 1243 1244 /* eval gradient */ 1245 SCIP_CALL( SCIPevalExprGradient(subproblem, expr, primalsol, 0L) ); 1246 1247 assert(SCIPexprGetDerivative(expr) != SCIP_INVALID); /* TODO this should be a proper check&abort */ /*lint !e777*/ 1248 1249 SCIP_CALL( SCIPfreeSol(subproblem, &primalsol) ); 1250 1251 /* update corresponding gradient entry */ 1252 SCIP_CALL( SCIPcreateExpriter(subproblem, &it) ); 1253 SCIP_CALL( SCIPexpriterInit(it, expr, SCIP_EXPRITER_DFS, FALSE) ); 1254 for( ; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*/ /*lint !e440*/ 1255 { 1256 if( !SCIPisExprVar(subproblem, expr) ) 1257 continue; 1258 1259 var = SCIPgetVarExprVar(expr); 1260 assert(var != NULL); 1261 1262 /* retrieving the master problem variable for the given subproblem variable. */ 1263 SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) ); 1264 if( mastervar == NULL ) 1265 continue; 1266 1267 assert(SCIPexprGetDerivative(expr) != SCIP_INVALID); /*lint !e777*/ 1268 coef = mult * SCIPexprGetDerivative(expr); 1269 1270 /* adding the variable to the storage */ 1271 SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) ); 1272 1273 *dirderiv += coef * getNlpVarSol(var, primalvals, var2idx); 1274 } 1275 SCIPfreeExpriter(&it); 1276 } 1277 1278 return SCIP_OKAY; 1279 } 1280