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