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_feas.c 26 * @ingroup OTHER_CFILES 27 * @brief Standard feasibility cuts for Benders' decomposition 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_feas.h" 35 #include "scip/benderscut_opt.h" 36 #include "scip/cons_linear.h" 37 #include "scip/pub_benderscut.h" 38 #include "scip/pub_benders.h" 39 #include "scip/pub_lp.h" 40 #include "scip/pub_message.h" 41 #include "scip/pub_misc.h" 42 #include "scip/pub_misc_linear.h" 43 #include "scip/pub_nlp.h" 44 #include "scip/pub_var.h" 45 #include "scip/scip_benders.h" 46 #include "scip/scip_cons.h" 47 #include "scip/scip_general.h" 48 #include "scip/scip_lp.h" 49 #include "scip/scip_mem.h" 50 #include "scip/scip_message.h" 51 #include "scip/scip_nlp.h" 52 #include "scip/scip_nlpi.h" 53 #include "scip/scip_numerics.h" 54 #include "scip/scip_prob.h" 55 #include "scip/scip_solvingstats.h" 56 #include "scip/scip_var.h" 57 58 #define BENDERSCUT_NAME "feas" 59 #define BENDERSCUT_DESC "Standard feasibility cuts for Benders' decomposition" 60 #define BENDERSCUT_PRIORITY 10000 61 #define BENDERSCUT_LPCUT TRUE 62 63 /* 64 * Local methods 65 */ 66 67 /** adds a variable and value to the constraint/row arrays */ 68 static 69 SCIP_RETCODE addVariableToArray( 70 SCIP* masterprob, /**< the SCIP instance of the master problem */ 71 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 72 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 73 SCIP_VAR* addvar, /**< the variable that will be added to the array */ 74 SCIP_Real addval, /**< the value that will be added to the array */ 75 int* nvars, /**< the number of variables in the variable array */ 76 int* varssize /**< the length of the variable size */ 77 ) 78 { 79 assert(masterprob != NULL); 80 assert(vars != NULL); 81 assert(*vars != NULL); 82 assert(vals != NULL); 83 assert(*vals != NULL); 84 assert(addvar != NULL); 85 assert(nvars != NULL); 86 assert(varssize != NULL); 87 88 if( *nvars >= *varssize ) 89 { 90 *varssize = SCIPcalcMemGrowSize(masterprob, *varssize + 1); 91 SCIP_CALL( SCIPreallocBufferArray(masterprob, vars, *varssize) ); 92 SCIP_CALL( SCIPreallocBufferArray(masterprob, vals, *varssize) ); 93 } 94 assert(*nvars < *varssize); 95 96 (*vars)[*nvars] = addvar; 97 (*vals)[*nvars] = addval; 98 (*nvars)++; 99 100 return SCIP_OKAY; 101 } 102 103 /** computing as standard Benders' feasibility cut from the dual solutions of the LP */ 104 static 105 SCIP_RETCODE computeStandardLPFeasibilityCut( 106 SCIP* masterprob, /**< the SCIP instance of the master problem */ 107 SCIP* subproblem, /**< the SCIP instance of the pricing problem */ 108 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 109 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 110 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 111 SCIP_Real* lhs, /**< the left hand side of the cut */ 112 int* nvars, /**< the number of variables in the cut */ 113 int* varssize, /**< the number of variables in the array */ 114 SCIP_Bool* success /**< was the cut generation successful? */ 115 ) 116 { 117 SCIP_VAR** subvars; 118 int nsubvars; 119 int nrows; 120 SCIP_Real dualsol; 121 SCIP_Real addval; /* the value that must be added to the lhs */ 122 int i; 123 124 assert(masterprob != NULL); 125 assert(subproblem != NULL); 126 assert(benders != NULL); 127 assert(SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_INFEASIBLE); 128 129 (*success) = FALSE; 130 131 /* looping over all LP rows and setting the coefficients of the cut */ 132 nrows = SCIPgetNLPRows(subproblem); 133 for( i = 0; i < nrows; i++ ) 134 { 135 SCIP_ROW* lprow; 136 137 lprow = SCIPgetLPRows(subproblem)[i]; 138 assert(lprow != NULL); 139 140 dualsol = SCIProwGetDualfarkas(lprow); 141 assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) ); 142 143 if( SCIPisDualfeasZero(subproblem, dualsol) ) 144 continue; 145 146 if( dualsol > 0.0 ) 147 addval = dualsol*SCIProwGetLhs(lprow); 148 else 149 addval = dualsol*SCIProwGetRhs(lprow); 150 151 *lhs += addval; 152 153 /* if the bound becomes infinite, then the cut generation terminates. */ 154 if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs) 155 || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval)) 156 { 157 (*success) = FALSE; 158 SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut.\n"); 159 return SCIP_OKAY; 160 } 161 } 162 163 nsubvars = SCIPgetNVars(subproblem); 164 subvars = SCIPgetVars(subproblem); 165 166 /* looping over all variables to update the coefficients in the computed cut. */ 167 for( i = 0; i < nsubvars; i++ ) 168 { 169 SCIP_VAR* var; 170 SCIP_VAR* mastervar; 171 172 var = subvars[i]; 173 174 /* retrieving the master problem variable for the given subproblem variable. */ 175 SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) ); 176 177 dualsol = SCIPgetVarFarkasCoef(subproblem, var); 178 179 if( SCIPisZero(subproblem, dualsol) ) 180 continue; 181 182 /* checking whether the original variable is a linking variable. 183 * If this is the case, then the corresponding master variable is added to the generated cut. 184 * If the pricing variable is not a linking variable, then the farkas dual value is added to the lhs 185 */ 186 if( mastervar != NULL ) 187 { 188 SCIPdebugMsg(masterprob ,"Adding coeffs to feasibility cut: <%s> dualsol %g\n", SCIPvarGetName(mastervar), dualsol); 189 190 /* adding the variable to the storage */ 191 SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, dualsol, nvars, varssize) ); 192 } 193 else 194 { 195 addval = 0; 196 197 if( SCIPisPositive(subproblem, dualsol) ) 198 addval = dualsol*SCIPvarGetUbGlobal(var); 199 else if( SCIPisNegative(subproblem, dualsol) ) 200 addval = dualsol*SCIPvarGetLbGlobal(var); 201 202 *lhs -= addval; 203 204 /* if the bound becomes infinite, then the cut generation terminates. */ 205 if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs) 206 || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval)) 207 { 208 (*success) = FALSE; 209 SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut.\n"); 210 return SCIP_OKAY; 211 } 212 } 213 } 214 215 (*success) = TRUE; 216 217 return SCIP_OKAY; 218 } 219 220 221 /** computing as standard Benders' feasibility cut from the dual solutions of the NLP 222 * 223 * NOTE: The cut must be created before being passed to this function 224 */ 225 static 226 SCIP_RETCODE computeStandardNLPFeasibilityCut( 227 SCIP* masterprob, /**< the SCIP instance of the master problem */ 228 SCIP* subproblem, /**< the SCIP instance of the pricing problem */ 229 SCIP_BENDERS* benders, /**< the benders' decomposition structure */ 230 SCIP_VAR*** vars, /**< pointer to array of variables in the generated cut with non-zero coefficient */ 231 SCIP_Real** vals, /**< pointer to array of coefficients of the variables in the generated cut */ 232 SCIP_Real* lhs, /**< the left hand side of the cut */ 233 int* nvars, /**< the number of variables in the cut */ 234 int* varssize, /**< the number of variables in the array */ 235 SCIP_Bool* success /**< was the cut generation successful? */ 236 ) 237 { 238 int nrows; 239 SCIP_Real activity; 240 SCIP_Real dirderiv; 241 SCIP_Real dualsol; 242 int i; 243 244 assert(masterprob != NULL); 245 assert(subproblem != NULL); 246 assert(benders != NULL); 247 assert(SCIPisNLPConstructed(subproblem)); 248 assert(SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_LOCINFEASIBLE || SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_GLOBINFEASIBLE); 249 250 (*success) = FALSE; 251 252 *lhs = 0.0; 253 dirderiv = 0.0; 254 255 /* looping over all NLP rows and setting the corresponding coefficients of the cut */ 256 nrows = SCIPgetNNLPNlRows(subproblem); 257 for( i = 0; i < nrows; i++ ) 258 { 259 SCIP_NLROW* nlrow; 260 261 nlrow = SCIPgetNLPNlRows(subproblem)[i]; 262 assert(nlrow != NULL); 263 264 dualsol = SCIPnlrowGetDualsol(nlrow); 265 assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) ); 266 267 if( SCIPisZero(subproblem, dualsol) ) 268 continue; 269 270 SCIP_CALL( SCIPaddNlRowGradientBenderscutOpt(masterprob, subproblem, benders, nlrow, -dualsol, 271 NULL, NULL, &dirderiv, vars, vals, nvars, varssize) ); 272 273 SCIP_CALL( SCIPgetNlRowActivity(subproblem, nlrow, &activity) ); 274 275 if( dualsol > 0.0 ) 276 { 277 assert(!SCIPisInfinity(subproblem, SCIPnlrowGetRhs(nlrow))); 278 *lhs += dualsol * (activity - SCIPnlrowGetRhs(nlrow)); 279 } 280 else 281 { 282 assert(!SCIPisInfinity(subproblem, -SCIPnlrowGetLhs(nlrow))); 283 *lhs += dualsol * (activity - SCIPnlrowGetLhs(nlrow)); 284 } 285 } 286 287 *lhs += dirderiv; 288 289 /* if the side became infinite or dirderiv was infinite, then the cut generation terminates. */ 290 if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs) 291 || SCIPisInfinity(masterprob, dirderiv) || SCIPisInfinity(masterprob, -dirderiv)) 292 { 293 (*success) = FALSE; 294 SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut. lhs = %g dirderiv = %g.\n", *lhs, dirderiv); 295 return SCIP_OKAY; 296 } 297 298 (*success) = TRUE; 299 300 return SCIP_OKAY; 301 } 302 303 /** generates and applies Benders' cuts */ 304 static 305 SCIP_RETCODE generateAndApplyBendersCuts( 306 SCIP* masterprob, /**< the SCIP instance of the master problem */ 307 SCIP* subproblem, /**< the SCIP instance of the pricing problem */ 308 SCIP_BENDERS* benders, /**< the benders' decomposition */ 309 SCIP_BENDERSCUT* benderscut, /**< the benders' decomposition cut method */ 310 SCIP_SOL* sol, /**< primal CIP solution */ 311 int probnumber, /**< the number of the pricing problem */ 312 SCIP_RESULT* result /**< the result from solving the subproblems */ 313 ) 314 { 315 SCIP_CONS* cut; 316 SCIP_VAR** vars; 317 SCIP_Real* vals; 318 SCIP_Real lhs; 319 SCIP_Real activity; 320 int nvars; 321 int varssize; 322 int nmastervars; 323 char cutname[SCIP_MAXSTRLEN]; 324 SCIP_Bool success; 325 326 assert(masterprob != NULL); 327 assert(subproblem != NULL); 328 assert(benders != NULL); 329 assert(result != NULL); 330 331 /* allocating memory for the variable and values arrays */ 332 nmastervars = SCIPgetNVars(masterprob) + SCIPgetNFixedVars(masterprob); 333 SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vars, nmastervars) ); 334 SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vals, nmastervars) ); 335 lhs = 0.0; 336 nvars = 0; 337 varssize = nmastervars; 338 339 /* setting the name of the generated cut */ 340 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "feasibilitycut_%d_%" SCIP_LONGINT_FORMAT, probnumber, 341 SCIPbenderscutGetNFound(benderscut) ); 342 343 if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) ) 344 { 345 /* computing the coefficients of the feasibility cut from the NLP */ 346 SCIP_CALL( computeStandardNLPFeasibilityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &nvars, &varssize, 347 &success) ); 348 } 349 else 350 { 351 if( SCIPgetNLPIterations(subproblem) == 0 ) 352 { 353 SCIPverbMessage(masterprob, SCIP_VERBLEVEL_FULL, NULL, "There were no iterations in pricing problem %d. " 354 "A Benders' decomposition feasibility cut will be generated from the presolved LP data.\n", probnumber); 355 } 356 357 /* computing the coefficients of the feasibility cut from the LP */ 358 SCIP_CALL( computeStandardLPFeasibilityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &nvars, &varssize, 359 &success) ); 360 } 361 362 /* if success is FALSE, then there was an error in generating the feasibility cut. No cut will be added to the master 363 * problem. Otherwise, the constraint is added to the master problem. 364 */ 365 if( !success ) 366 { 367 (*result) = SCIP_DIDNOTFIND; 368 SCIPdebugMsg(masterprob, "Error in generating Benders' feasibility cut for problem %d.\n", probnumber); 369 } 370 else 371 { 372 /* creating a constraint with the variables and coefficients previously stored */ 373 SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cut, cutname, nvars, vars, vals, lhs, SCIPinfinity(masterprob)) ); 374 SCIP_CALL( SCIPsetConsDynamic(masterprob, cut, TRUE) ); 375 SCIP_CALL( SCIPsetConsRemovable(masterprob, cut, TRUE) ); 376 377 assert(SCIPisInfinity(masterprob, SCIPgetRhsLinear(masterprob, cut))); 378 379 /* the activity of the cut should be less than the lhs. This will ensure that the evaluated solution will be cut off. 380 * It is possible that the activity is greater than the lhs. This could be caused by numerical difficulties. In this 381 * case, no cut will be generated. 382 */ 383 lhs = SCIPgetLhsLinear(masterprob, cut); 384 activity = SCIPgetActivityLinear(masterprob, cut, sol); 385 if( SCIPisGE(masterprob, activity, lhs) ) 386 { 387 success = FALSE; 388 SCIPdebugMsg(masterprob ,"Invalid feasibility cut - activity is greater than lhs %g >= %g.\n", activity, lhs); 389 #ifdef SCIP_DEBUG 390 SCIPABORT(); 391 #endif 392 } 393 394 assert(cut != NULL); 395 396 if( success ) 397 { 398 /* adding the constraint to the master problem */ 399 SCIP_CALL( SCIPaddCons(masterprob, cut) ); 400 401 SCIPdebugPrintCons(masterprob, cut, NULL); 402 403 (*result) = SCIP_CONSADDED; 404 } 405 406 SCIP_CALL( SCIPreleaseCons(masterprob, &cut) ); 407 } 408 409 SCIPfreeBufferArray(masterprob, &vals); 410 SCIPfreeBufferArray(masterprob, &vars); 411 412 return SCIP_OKAY; 413 } 414 415 /* 416 * Callback methods of Benders' decomposition cuts 417 */ 418 419 /** execution method of Benders' decomposition cuts */ 420 static 421 SCIP_DECL_BENDERSCUTEXEC(benderscutExecFeas) 422 { /*lint --e{715}*/ 423 SCIP* subproblem; 424 SCIP_Bool nlprelaxation; 425 426 assert(scip != NULL); 427 assert(benders != NULL); 428 assert(benderscut != NULL); 429 assert(result != NULL); 430 assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders)); 431 432 subproblem = SCIPbendersSubproblem(benders, probnumber); 433 434 if( subproblem == NULL ) 435 { 436 SCIPdebugMsg(scip, "The subproblem %d is set to NULL. The <%s> Benders' decomposition cut can not be executed.\n", 437 probnumber, BENDERSCUT_NAME); 438 439 (*result) = SCIP_DIDNOTRUN; 440 return SCIP_OKAY; 441 } 442 443 /* setting a flag to indicate whether the NLP relaxation should be used to generate cuts */ 444 nlprelaxation = SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem); 445 446 /* only generate feasibility cuts if the subproblem LP or NLP is infeasible, 447 * since we use the farkas proof from the LP or the dual solution of the NLP to construct the feasibility cut 448 */ 449 if( SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING && 450 ((!nlprelaxation && SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_INFEASIBLE) || 451 (nlprelaxation && (SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_LOCINFEASIBLE || SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_GLOBINFEASIBLE))) ) 452 { 453 /* generating a cut for a given subproblem */ 454 SCIP_CALL( generateAndApplyBendersCuts(scip, subproblem, benders, benderscut, 455 sol, probnumber, result) ); 456 } 457 458 return SCIP_OKAY; 459 } 460 461 462 /* 463 * Benders' decomposition cuts specific interface methods 464 */ 465 466 /** creates the Standard Feasibility Benders' decomposition cuts and includes it in SCIP */ 467 SCIP_RETCODE SCIPincludeBenderscutFeas( 468 SCIP* scip, /**< SCIP data structure */ 469 SCIP_BENDERS* benders /**< Benders' decomposition */ 470 ) 471 { 472 SCIP_BENDERSCUT* benderscut; 473 474 assert(benders != NULL); 475 476 benderscut = NULL; 477 478 /* include Benders' decomposition cuts */ 479 SCIP_CALL( SCIPincludeBenderscutBasic(scip, benders, &benderscut, BENDERSCUT_NAME, BENDERSCUT_DESC, 480 BENDERSCUT_PRIORITY, BENDERSCUT_LPCUT, benderscutExecFeas, NULL) ); 481 482 assert(benderscut != NULL); 483 484 return SCIP_OKAY; 485 } 486