Scippy

SCIP

Solving Constraint Integer Programs

benders.c
Go to the documentation of this file.
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 benders.c
17  * @ingroup OTHER_CFILES
18  * @brief methods 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 <assert.h>
25 #include <string.h>
26 
27 #include "scip/def.h"
28 #include "scip/set.h"
29 #include "scip/clock.h"
30 #include "scip/dcmp.h"
31 #include "scip/paramset.h"
32 #include "scip/lp.h"
33 #include "scip/prob.h"
34 #include "scip/pricestore.h"
35 #include "scip/scip.h"
36 #include "scip/scipdefplugins.h"
37 #include "scip/benders.h"
38 #include "scip/pub_message.h"
39 #include "scip/pub_misc.h"
40 #include "scip/cons_linear.h"
41 #include "scip/cons_nonlinear.h"
42 
43 #include "scip/struct_benders.h"
44 #include "scip/struct_benderscut.h"
45 
46 #include "scip/benderscut.h"
47 
48 /* Defaults for parameters */
49 #define SCIP_DEFAULT_TRANSFERCUTS FALSE /** should Benders' cuts generated in LNS heuristics be transferred to the main SCIP instance? */
50 #define SCIP_DEFAULT_CUTSASCONSS TRUE /** should the transferred cuts be added as constraints? */
51 #define SCIP_DEFAULT_LNSCHECK TRUE /** should the Benders' decomposition be used in LNS heuristics */
52 #define SCIP_DEFAULT_LNSMAXDEPTH -1 /** maximum depth at which the LNS check is performed */
53 #define SCIP_DEFAULT_LNSMAXCALLS 10 /** the maximum number of Benders' decomposition calls in LNS heuristics */
54 #define SCIP_DEFAULT_LNSMAXCALLSROOT 0 /** the maximum number of root node Benders' decomposition calls in LNS heuristics */
55 #define SCIP_DEFAULT_SUBPROBFRAC 1.0 /** fraction of subproblems that are solved in each iteration */
56 #define SCIP_DEFAULT_UPDATEAUXVARBOUND FALSE /** should the auxiliary variable lower bound be updated by solving the subproblem */
57 #define SCIP_DEFAULT_AUXVARSIMPLINT FALSE /** set the auxiliary variables as implint if the subproblem objective is integer */
58 #define SCIP_DEFAULT_CUTCHECK TRUE /** should cuts be generated during the checking of solutions? */
59 #define SCIP_DEFAULT_STRENGTHENMULT 0.5 /** the convex combination multiplier for the cut strengthening */
60 #define SCIP_DEFAULT_NOIMPROVELIMIT 5 /** the maximum number of cut strengthening without improvement */
61 #define SCIP_DEFAULT_STRENGTHENPERTURB 1e-06 /** the amount by which the cut strengthening solution is perturbed */
62 #define SCIP_DEFAULT_STRENGTHENENABLED FALSE /** enable the core point cut strengthening approach */
63 #define SCIP_DEFAULT_STRENGTHENINTPOINT 'r' /** where should the strengthening interior point be sourced from ('l'p relaxation, 'f'irst solution, 'i'ncumbent solution, 'r'elative interior point, vector of 'o'nes, vector of 'z'eros) */
64 #define SCIP_DEFAULT_NUMTHREADS 1 /** the number of parallel threads to use when solving the subproblems */
65 #define SCIP_DEFAULT_EXECFEASPHASE FALSE /** should a feasibility phase be executed during the root node processing */
66 #define SCIP_DEFAULT_SLACKVARCOEF 1e+6 /** the initial objective coefficient of the slack variables in the subproblem */
67 #define SCIP_DEFAULT_MAXSLACKVARCOEF 1e+9 /** the maximal objective coefficient of the slack variables in the subproblem */
68 #define SCIP_DEFAULT_CHECKCONSCONVEXITY TRUE /** should the constraints of the subproblem be checked for convexity? */
69 
70 #define BENDERS_MAXPSEUDOSOLS 5 /** the maximum number of pseudo solutions checked before suggesting
71  * merge candidates */
72 
73 #define BENDERS_ARRAYSIZE 1000 /**< the initial size of the added constraints/cuts arrays */
74 
75 #define AUXILIARYVAR_NAME "##bendersauxiliaryvar" /** the name for the Benders' auxiliary variables in the master problem */
76 #define SLACKVAR_NAME "##bendersslackvar" /** the name for the Benders' slack variables added to each
77  * constraints in the subproblems */
78 #define NLINEARCONSHDLRS 5
79 
80 /* event handler properties */
81 #define NODEFOCUS_EVENTHDLR_NAME "bendersnodefocus"
82 #define NODEFOCUS_EVENTHDLR_DESC "node focus event handler for Benders' decomposition"
83 
84 #define MIPNODEFOCUS_EVENTHDLR_NAME "bendersmipsolvenodefocus"
85 #define MIPNODEFOCUS_EVENTHDLR_DESC "node focus event handler for the MIP solve method for Benders' decomposition"
86 
87 #define UPPERBOUND_EVENTHDLR_NAME "bendersupperbound"
88 #define UPPERBOUND_EVENTHDLR_DESC "found solution event handler to terminate subproblem solve for a given upper bound"
89 
90 #define NODESOLVED_EVENTHDLR_NAME "bendersnodesolved"
91 #define NODESOLVED_EVENTHDLR_DESC "node solved event handler for the Benders' integer cuts"
92 
93 
94 /** event handler data */
95 struct SCIP_EventhdlrData
96 {
97  int filterpos; /**< the event filter entry */
98  int numruns; /**< the number of times that the problem has been solved */
99  SCIP_Real upperbound; /**< an upper bound for the problem */
100  SCIP_Bool solvecip; /**< is the event called from a MIP subproblem solve*/
101 };
102 
103 
104 /* ---------------- Local methods for event handlers ---------------- */
105 
106 /** initialises the members of the eventhandler data */
107 static
109  SCIP* scip, /**< the SCIP data structure */
110  SCIP_EVENTHDLRDATA* eventhdlrdata /**< the event handler data */
111  )
112 {
113  assert(scip != NULL);
114  assert(eventhdlrdata != NULL);
115 
116  eventhdlrdata->filterpos = -1;
117  eventhdlrdata->numruns = 0;
118  eventhdlrdata->upperbound = -SCIPinfinity(scip);
119  eventhdlrdata->solvecip = FALSE;
120 
121  return SCIP_OKAY;
122 }
123 
124 /** initsol method for the event handlers */
125 static
127  SCIP* scip, /**< the SCIP data structure */
128  SCIP_EVENTHDLR* eventhdlr, /**< the event handlers data structure */
129  SCIP_EVENTTYPE eventtype /**< event type mask to select events to catch */
130  )
131 {
132  SCIP_EVENTHDLRDATA* eventhdlrdata;
133 
134  assert(scip != NULL);
135  assert(eventhdlr != NULL);
136 
137  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
138 
139  SCIP_CALL(SCIPcatchEvent(scip, eventtype, eventhdlr, NULL, &eventhdlrdata->filterpos));
140 
141  return SCIP_OKAY;
142 }
143 
144 /** the exit sol method for the event handlers */
145 static
147  SCIP* scip, /**< the SCIP data structure */
148  SCIP_EVENTHDLR* eventhdlr, /**< the event handlers data structure */
149  SCIP_EVENTTYPE eventtype /**< event type mask to select events to catch */
150  )
151 {
152  SCIP_EVENTHDLRDATA* eventhdlrdata;
153 
154  assert(scip != NULL);
155  assert(eventhdlr != NULL);
156 
157  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
158 
159  if( eventhdlrdata->filterpos >= 0 )
160  {
161  SCIP_CALL(SCIPdropEvent(scip, eventtype, eventhdlr, NULL, eventhdlrdata->filterpos));
162  eventhdlrdata->filterpos = -1;
163  }
164 
165  return SCIP_OKAY;
166 }
167 
168 /** the exit method for the event handlers */
169 static
171  SCIP* scip, /**< the SCIP data structure */
172  SCIP_EVENTHDLR* eventhdlr /**< the event handlers data structure */
173  )
174 {
175  SCIP_EVENTHDLRDATA* eventhdlrdata;
176 
177  assert(scip != NULL);
178  assert(eventhdlr != NULL);
179 
180  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
181 
182  /* reinitialise the event handler data */
183  SCIP_CALL( initEventhandlerData(scip, eventhdlrdata) );
184 
185  return SCIP_OKAY;
186 }
187 
188 /** free method for the event handler */
189 static
191  SCIP* scip, /**< the SCIP data structure */
192  SCIP_EVENTHDLR* eventhdlr /**< the event handlers data structure */
193  )
194 {
195  SCIP_EVENTHDLRDATA* eventhdlrdata;
196 
197  assert(scip != NULL);
198  assert(eventhdlr != NULL);
199 
200  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
201  assert(eventhdlrdata != NULL);
202 
203  SCIPfreeBlockMemory(scip, &eventhdlrdata);
204 
205  SCIPeventhdlrSetData(eventhdlr, NULL);
206 
207  return SCIP_OKAY;
208 }
209 
210 
211 
212 /* ---------------- Callback methods of node focus event handler ---------------- */
213 
214 /** exec the event handler */
215 static
216 SCIP_DECL_EVENTEXEC(eventExecBendersNodefocus)
217 { /*lint --e{715}*/
218  SCIP_EVENTHDLRDATA* eventhdlrdata;
219 
220  assert(scip != NULL);
221  assert(eventhdlr != NULL);
222  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODEFOCUS_EVENTHDLR_NAME) == 0);
223 
224  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
225 
226  /* sending an interrupt solve signal to return the control back to the Benders' decomposition plugin.
227  * This will ensure the SCIP stage is SCIP_STAGE_SOLVING, allowing the use of probing mode. */
229 
230  SCIP_CALL(SCIPdropEvent(scip, SCIP_EVENTTYPE_NODEFOCUSED, eventhdlr, NULL, eventhdlrdata->filterpos));
231  eventhdlrdata->filterpos = -1;
232 
233  return SCIP_OKAY;
234 }
235 
236 /** solving process initialization method of event handler (called when branch and bound process is about to begin) */
237 static
238 SCIP_DECL_EVENTINITSOL(eventInitsolBendersNodefocus)
239 {
240  assert(scip != NULL);
241  assert(eventhdlr != NULL);
242  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODEFOCUS_EVENTHDLR_NAME) == 0);
243 
245 
246  return SCIP_OKAY;
247 }
248 
249 /** solving process deinitialization method of event handler (called before branch and bound process data is freed) */
250 static
251 SCIP_DECL_EVENTEXITSOL(eventExitsolBendersNodefocus)
252 {
253  assert(scip != NULL);
254  assert(eventhdlr != NULL);
255  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODEFOCUS_EVENTHDLR_NAME) == 0);
256 
258 
259  return SCIP_OKAY;
260 }
261 
262 /** deinitialization method of event handler (called before transformed problem is freed) */
263 static
264 SCIP_DECL_EVENTEXIT(eventExitBendersNodefocus)
265 {
266  assert(scip != NULL);
267  assert(eventhdlr != NULL);
268  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODEFOCUS_EVENTHDLR_NAME) == 0);
269 
270  SCIP_CALL( exitEventhandler(scip, eventhdlr) );
271 
272  return SCIP_OKAY;
273 }
274 
275 /** deinitialization method of event handler (called before transformed problem is freed) */
276 static
277 SCIP_DECL_EVENTFREE(eventFreeBendersNodefocus)
278 {
279  assert(scip != NULL);
280  assert(eventhdlr != NULL);
281  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODEFOCUS_EVENTHDLR_NAME) == 0);
282 
283  SCIP_CALL( freeEventhandler(scip, eventhdlr) );
284 
285  return SCIP_OKAY;
286 }
287 
288 
289 /* ---------------- Callback methods of MIP solve node focus event handler ---------------- */
290 
291 /** exec the event handler */
292 static
293 SCIP_DECL_EVENTEXEC(eventExecBendersMipnodefocus)
294 { /*lint --e{715}*/
295  SCIP_EVENTHDLRDATA* eventhdlrdata;
296 
297  assert(scip != NULL);
298  assert(eventhdlr != NULL);
299  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), MIPNODEFOCUS_EVENTHDLR_NAME) == 0);
300 
301  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
302 
303  /* interrupting the solve so that the control is returned back to the Benders' core. */
304  if( eventhdlrdata->numruns == 0 && !eventhdlrdata->solvecip )
305  {
307  }
308 
309  SCIP_CALL(SCIPdropEvent(scip, SCIP_EVENTTYPE_NODEFOCUSED, eventhdlr, NULL, eventhdlrdata->filterpos));
310  eventhdlrdata->filterpos = -1;
311 
312  eventhdlrdata->numruns++;
313 
314  return SCIP_OKAY;
315 }
316 
317 /** solving process initialization method of event handler (called when branch and bound process is about to begin) */
318 static
319 SCIP_DECL_EVENTINITSOL(eventInitsolBendersMipnodefocus)
320 {
321  assert(scip != NULL);
322  assert(eventhdlr != NULL);
323  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), MIPNODEFOCUS_EVENTHDLR_NAME) == 0);
324 
326 
327  return SCIP_OKAY;
328 }
329 
330 /** solving process deinitialization method of event handler (called before branch and bound process data is freed) */
331 static
332 SCIP_DECL_EVENTEXITSOL(eventExitsolBendersMipnodefocus)
333 {
334  assert(scip != NULL);
335  assert(eventhdlr != NULL);
336  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), MIPNODEFOCUS_EVENTHDLR_NAME) == 0);
337 
339 
340  return SCIP_OKAY;
341 }
342 
343 /** deinitialization method of event handler (called before transformed problem is freed) */
344 static
345 SCIP_DECL_EVENTEXIT(eventExitBendersMipnodefocus)
346 {
347  assert(scip != NULL);
348  assert(eventhdlr != NULL);
349  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), MIPNODEFOCUS_EVENTHDLR_NAME) == 0);
350 
351  SCIP_CALL( exitEventhandler(scip, eventhdlr) );
352 
353  return SCIP_OKAY;
354 }
355 
356 /** deinitialization method of event handler (called before transformed problem is freed) */
357 static
358 SCIP_DECL_EVENTFREE(eventFreeBendersMipnodefocus)
359 {
360  assert(scip != NULL);
361  assert(eventhdlr != NULL);
362  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), MIPNODEFOCUS_EVENTHDLR_NAME) == 0);
363 
364  SCIP_CALL( freeEventhandler(scip, eventhdlr) );
365 
366  return SCIP_OKAY;
367 }
368 
369 /* ---------------- Callback methods of solution found event handler ---------------- */
370 
371 /** exec the event handler */
372 static
373 SCIP_DECL_EVENTEXEC(eventExecBendersUpperbound)
374 { /*lint --e{715}*/
375  SCIP_EVENTHDLRDATA* eventhdlrdata;
376  SCIP_SOL* bestsol;
377 
378  assert(scip != NULL);
379  assert(eventhdlr != NULL);
380  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), UPPERBOUND_EVENTHDLR_NAME) == 0);
381 
382  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
383  assert(eventhdlrdata != NULL);
384 
385  bestsol = SCIPgetBestSol(scip);
386 
387  if( SCIPisLT(scip, SCIPgetSolOrigObj(scip, bestsol)*(int)SCIPgetObjsense(scip), eventhdlrdata->upperbound) )
388  {
390  }
391 
392  return SCIP_OKAY;
393 }
394 
395 /** solving process initialization method of event handler (called when branch and bound process is about to begin) */
396 static
397 SCIP_DECL_EVENTINITSOL(eventInitsolBendersUpperbound)
398 {
399  assert(scip != NULL);
400  assert(eventhdlr != NULL);
401  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), UPPERBOUND_EVENTHDLR_NAME) == 0);
402 
404 
405  return SCIP_OKAY;
406 }
407 
408 /** solving process deinitialization method of event handler (called before branch and bound process data is freed) */
409 static
410 SCIP_DECL_EVENTEXITSOL(eventExitsolBendersUpperbound)
411 {
412  assert(scip != NULL);
413  assert(eventhdlr != NULL);
414  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), UPPERBOUND_EVENTHDLR_NAME) == 0);
415 
417 
418  return SCIP_OKAY;
419 }
420 
421 /** deinitialization method of event handler (called before transformed problem is freed) */
422 static
423 SCIP_DECL_EVENTEXIT(eventExitBendersUpperbound)
424 {
425  assert(scip != NULL);
426  assert(eventhdlr != NULL);
427  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), UPPERBOUND_EVENTHDLR_NAME) == 0);
428 
429  SCIP_CALL( exitEventhandler(scip, eventhdlr) );
430 
431  return SCIP_OKAY;
432 }
433 
434 /** deinitialization method of event handler (called before transformed problem is freed) */
435 static
436 SCIP_DECL_EVENTFREE(eventFreeBendersUpperbound)
437 {
438  assert(scip != NULL);
439  assert(eventhdlr != NULL);
440  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), UPPERBOUND_EVENTHDLR_NAME) == 0);
441 
442  SCIP_CALL( freeEventhandler(scip, eventhdlr) );
443 
444  return SCIP_OKAY;
445 }
446 
447 /** updates the upper bound in the event handler data */
448 static
450  SCIP_BENDERS* benders, /**< Benders' decomposition */
451  int probnumber, /**< the subproblem number */
452  SCIP_Real upperbound /**< the upper bound value */
453  )
454 {
455  SCIP_EVENTHDLR* eventhdlr;
456  SCIP_EVENTHDLRDATA* eventhdlrdata;
457 
458  assert(benders != NULL);
459  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
460 
461  eventhdlr = SCIPfindEventhdlr(SCIPbendersSubproblem(benders, probnumber), UPPERBOUND_EVENTHDLR_NAME);
462  assert(eventhdlr != NULL);
463 
464  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
465  assert(eventhdlrdata != NULL);
466 
467  eventhdlrdata->upperbound = upperbound;
468 
469  return SCIP_OKAY;
470 }
471 
472 /* ---------------- Callback methods of the node solved event handler ---------------- */
473 
474 /** Updates the cut constant of the Benders' cuts data.
475  * This function solves the master problem with only the auxiliary variables in the objective function.
476  */
477 static
479  SCIP* masterprob, /**< the SCIP instance of the master problem */
480  SCIP_BENDERS* benders /**< Benders' decomposition */
481  )
482 {
483  SCIP_VAR** vars;
484  int nvars;
485  int nsubproblems;
486  int i;
487  SCIP_Bool lperror;
488  SCIP_Bool cutoff;
489 
490  assert(masterprob != NULL);
491  assert(benders != NULL);
492 
493  /* don't run in probing or in repropagation */
494  if( SCIPinProbing(masterprob) || SCIPinRepropagation(masterprob) || SCIPinDive(masterprob) )
495  return SCIP_OKAY;
496 
497  nsubproblems = SCIPbendersGetNSubproblems(benders);
498 
499  SCIP_CALL( SCIPstartProbing(masterprob) );
500 
501  /* change the master problem variables to 0 */
502  nvars = SCIPgetNVars(masterprob);
503  vars = SCIPgetVars(masterprob);
504 
505  /* setting the objective function coefficient to 0 for all variables */
506  for( i = 0; i < nvars; i++ )
507  {
508  if( SCIPvarGetStatus(vars[i]) == SCIP_VARSTATUS_COLUMN )
509  {
510  SCIP_CALL( SCIPchgVarObjProbing(masterprob, vars[i], 0.0) );
511  }
512  }
513 
514  /* solving an LP for all subproblems to find the lower bound */
515  for( i = 0; i < nsubproblems; i++)
516  {
517  SCIP_VAR* auxiliaryvar;
518 
519  auxiliaryvar = SCIPbendersGetAuxiliaryVar(benders, i);
520 
521  if( SCIPvarGetStatus(auxiliaryvar) != SCIP_VARSTATUS_COLUMN )
522  continue;
523 
524  SCIP_CALL( SCIPchgVarObjProbing(masterprob, auxiliaryvar, 1.0) );
525 
526  /* solving the probing LP to get a lower bound on the auxiliary variables */
527  SCIP_CALL( SCIPsolveProbingLP(masterprob, -1, &lperror, &cutoff) );
528 
529  if( !SCIPisInfinity(masterprob, -SCIPgetSolTransObj(masterprob, NULL)) )
531 
532  SCIPdebugMsg(masterprob, "Cut constant for subproblem %d: %g\n", i,
534 
535  SCIP_CALL( SCIPchgVarObjProbing(masterprob, auxiliaryvar, 0.0) );
536  }
537 
538  SCIP_CALL( SCIPendProbing(masterprob) );
539 
540  return SCIP_OKAY;
541 }
542 
543 /** exec the event handler */
544 static
545 SCIP_DECL_EVENTEXEC(eventExecBendersNodesolved)
546 { /*lint --e{715}*/
547  SCIP_BENDERS* benders;
548 
549  assert(scip != NULL);
550  assert(eventhdlr != NULL);
551  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODESOLVED_EVENTHDLR_NAME) == 0);
552 
553  benders = (SCIP_BENDERS*)SCIPeventhdlrGetData(eventhdlr); /*lint !e826*/
554 
555  if( SCIPbendersGetNSubproblems(benders) > 0
557  {
559  }
560 
562 
563  return SCIP_OKAY;
564 }
565 
566 /** solving process initialization method of event handler (called when branch and bound process is about to begin) */
567 static
568 SCIP_DECL_EVENTINITSOL(eventInitsolBendersNodesolved)
569 {
570  SCIP_BENDERS* benders;
571 
572  assert(scip != NULL);
573  assert(eventhdlr != NULL);
574  assert(strcmp(SCIPeventhdlrGetName(eventhdlr), NODESOLVED_EVENTHDLR_NAME) == 0);
575 
576  /* getting the Benders' decomposition data structure */
577  benders = (SCIP_BENDERS*)SCIPeventhdlrGetData(eventhdlr); /*lint !e826*/
578 
579  /* The event is only caught if there is an active Benders' decomposition, the integer subproblem are solved and
580  * the Benders' decomposition has not been copied in thread safe mode
581  */
583  && !benders->threadsafe )
584  {
586  }
587 
588  return SCIP_OKAY;
589 }
590 
591 
592 /* ---------------- Methods for the parallelisation of Benders' decomposition ---------------- */
593 
594 /** comparison method for sorting the subproblems.
595  * The subproblem that has been called the least is prioritised
596  */
597 static
598 SCIP_DECL_SORTPTRCOMP(benderssubcompdefault)
599 {
600  SCIP_SUBPROBLEMSOLVESTAT* solvestat1;
601  SCIP_SUBPROBLEMSOLVESTAT* solvestat2;
602 
603  assert(elem1 != NULL);
604  assert(elem2 != NULL);
605 
606  solvestat1 = (SCIP_SUBPROBLEMSOLVESTAT*)elem1;
607  solvestat2 = (SCIP_SUBPROBLEMSOLVESTAT*)elem2;
608 
609  /* prefer subproblems with fewer calls, using the index as tie breaker */
610  if( MAX(solvestat1->ncalls, solvestat2->ncalls) == 0 )
611  return solvestat1->idx - solvestat2->idx;
612  else if( solvestat1->ncalls != solvestat2->ncalls )
613  return solvestat1->ncalls - solvestat2->ncalls;
614  else
615  {
616  /* prefer the harder problem (with more average iterations) */
617  int avgiterdiff = (int)solvestat2->avgiter - (int)solvestat1->avgiter;
618 
619  if( avgiterdiff != 0 )
620  return avgiterdiff;
621 
622  return solvestat1->idx - solvestat2->idx;
623  }
624 
625 /* the code below does not give a total order of the elements */
626 #ifdef SCIP_DISABLED_CODE
627  if( solvestat1->ncalls == 0 )
628  if( solvestat2->ncalls == 0 )
629  if( solvestat1->idx < solvestat2->idx )
630  return -1;
631  else
632  return 1;
633  else
634  return -1;
635  else if( solvestat2->ncalls == 0 )
636  return 1;
637  else
638  {
639  if( solvestat1->ncalls < solvestat2->ncalls )
640  return -1;
641  else if( solvestat2->ncalls < solvestat1->ncalls )
642  return 1;
643  else
644  {
645  /* we want to execute the hard subproblems first */
646  if( solvestat1->avgiter > solvestat2->avgiter )
647  return 1;
648  else
649  return -1;
650  }
651  }
652 #endif
653 }
654 
655 /* Local methods */
656 
657 /** A workaround for GCG. This is a temp vardata that is set for the auxiliary variables */
658 struct SCIP_VarData
659 {
660  int vartype; /**< the variable type. In GCG this indicates whether the variable is a
661  * master problem or subproblem variable. */
662 };
663 
664 /** adds the auxiliary variables to the Benders' decomposition master problem */
665 static
667  SCIP* scip, /**< SCIP data structure */
668  SCIP_BENDERS* benders /**< Benders' decomposition structure */
669  )
670 {
671  SCIP_BENDERS* topbenders; /* the highest priority Benders' decomposition */
672  SCIP_VAR* auxiliaryvar;
673  SCIP_VARDATA* vardata;
674  char varname[SCIP_MAXSTRLEN]; /* the name of the auxiliary variable */
675  SCIP_Bool shareauxvars;
676  int i;
677 
678  /* this is a workaround for GCG. GCG expects that the variable has vardata when added. So a dummy vardata is created */
679  SCIP_CALL( SCIPallocBlockMemory(scip, &vardata) );
680  vardata->vartype = -1;
681 
682  /* getting the highest priority Benders' decomposition */
683  topbenders = SCIPgetBenders(scip)[0];
684 
685  /* if the current Benders is the highest priority Benders, then we need to create the auxiliary variables.
686  * Otherwise, if the shareauxvars flag is set, then the auxiliary variables from the highest priority Benders' are
687  * stored with this Benders. */
688  shareauxvars = FALSE;
689  if( topbenders != benders && SCIPbendersShareAuxVars(benders) )
690  shareauxvars = TRUE;
691 
692  for( i = 0; i < SCIPbendersGetNSubproblems(benders); i++ )
693  {
694  /* if the auxiliary variables are shared, then a pointer to the variable is retrieved from topbenders,
695  * otherwise the auxiliaryvariable is created. */
696  if( shareauxvars )
697  {
698  auxiliaryvar = SCIPbendersGetAuxiliaryVar(topbenders, i);
699 
700  SCIP_CALL( SCIPcaptureVar(scip, auxiliaryvar) );
701  }
702  else
703  {
704  SCIP_VARTYPE vartype;
705 
706  /* set the variable type of the auxiliary variables to implied integer if the objective function of the
707  * subproblem is guaranteed to be integer. This behaviour is controlled through a user parameter.
708  * NOTE: It is only possible to determine if the objective function is integral if the subproblem is defined as
709  * a SCIP instance, i.e. not NULL.
710  */
711  if( benders->auxvarsimplint && SCIPbendersSubproblem(benders, i) != NULL
712  && SCIPisObjIntegral(SCIPbendersSubproblem(benders, i)) )
713  vartype = SCIP_VARTYPE_IMPLINT;
714  else
715  vartype = SCIP_VARTYPE_CONTINUOUS;
716 
717  (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "%s_%d_%s", AUXILIARYVAR_NAME, i, SCIPbendersGetName(benders) );
718  SCIP_CALL( SCIPcreateVarBasic(scip, &auxiliaryvar, varname, benders->subproblowerbound[i], SCIPinfinity(scip),
719  1.0, vartype) );
720 
721  SCIPvarSetData(auxiliaryvar, vardata);
722 
723  SCIP_CALL( SCIPaddVar(scip, auxiliaryvar) );
724 
725  /* adding the down lock for the Benders' decomposition constraint handler */
726  SCIP_CALL( SCIPaddVarLocksType(scip, auxiliaryvar, SCIP_LOCKTYPE_MODEL, 1, 0) );
727  }
728 
729  benders->auxiliaryvars[i] = auxiliaryvar;
730  }
731 
732  SCIPfreeBlockMemory(scip, &vardata);
733 
734  return SCIP_OKAY;
735 }
736 
737 /** assigns the copied auxiliary variables in the target SCIP to the target Benders' decomposition data */
738 static
740  SCIP* scip, /**< SCIP data structure, the target scip */
741  SCIP_BENDERS* benders /**< Benders' decomposition */
742  )
743 {
744  SCIP_BENDERS* topbenders; /* the highest priority Benders' decomposition */
745  SCIP_VAR* targetvar;
746  SCIP_VARDATA* vardata;
747  char varname[SCIP_MAXSTRLEN]; /* the name of the auxiliary variable */
748  SCIP_Bool shareauxvars;
749  int subscipdepth;
750  int i;
751  int j;
752 
753  assert(scip != NULL);
754  assert(benders != NULL);
755 
756  /* this is a workaround for GCG. GCG expects that the variable has vardata when added. So a dummy vardata is created */
757  SCIP_CALL( SCIPallocBlockMemory(scip, &vardata) );
758  vardata->vartype = -1;
759 
760  /* getting the highest priority Benders' decomposition */
761  topbenders = SCIPgetBenders(scip)[0];
762 
763  /* if the auxiliary variable are shared, then the variable name will have a suffix of the highest priority Benders'
764  * name. So the shareauxvars flag indicates how to search for the auxiliary variables */
765  shareauxvars = FALSE;
766  if( topbenders != benders && SCIPbendersShareAuxVars(benders) )
767  shareauxvars = TRUE;
768 
769  subscipdepth = SCIPgetSubscipDepth(scip);
770 
771  for( i = 0; i < SCIPbendersGetNSubproblems(benders); i++ )
772  {
773  char prefix[SCIP_MAXSTRLEN];
774  char tmpprefix[SCIP_MAXSTRLEN];
775  int len = 1;
776 
777  j = 0;
778  targetvar = NULL;
779 
780  /* the prefix for the variable names is required for UG, since we don't know how many copies have been made. To
781  * find the target variable, we start with an empty prefix. Then t_ is prepended until the target variable is
782  * found
783  */
784  prefix[0] = '\0';
785  while( targetvar == NULL && j <= subscipdepth )
786  {
787  if( shareauxvars )
788  (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "%s%s_%d_%s", prefix, AUXILIARYVAR_NAME, i, SCIPbendersGetName(topbenders));
789  else
790  (void) SCIPsnprintf(varname, SCIP_MAXSTRLEN, "%s%s_%d_%s", prefix, AUXILIARYVAR_NAME, i, SCIPbendersGetName(benders));
791 
792  /* finding the variable in the copied problem that has the same name as the auxiliary variable */
793  targetvar = SCIPfindVar(scip, varname);
794 
795  (void) SCIPsnprintf(tmpprefix, len, "t_%s", prefix);
796  len += 2;
797  (void) strncpy(prefix, tmpprefix, len); /*lint !e732*/
798 
799  j++;
800  }
801 
802  if( targetvar != NULL )
803  {
804  SCIPvarSetData(targetvar, vardata);
805 
806  benders->auxiliaryvars[i] = SCIPvarGetTransVar(targetvar);
807 
808  SCIP_CALL( SCIPcaptureVar(scip, benders->auxiliaryvars[i]) );
809  }
810  else
811  {
812  SCIPABORT();
813  }
814  }
815 
816  SCIPfreeBlockMemory(scip, &vardata);
817 
818  return SCIP_OKAY;
819 }
820 
821 /** sets the subproblem objective value array to -infinity */
822 static
824  SCIP_BENDERS* benders, /**< the Benders' decomposition structure */
825  SCIP_SET* set /**< global SCIP settings */
826  )
827 {
828  SCIP* subproblem;
829  SCIP_Real inf;
830  int nsubproblems;
831  int i;
832 
833  assert(benders != NULL);
834 
835  nsubproblems = SCIPbendersGetNSubproblems(benders);
836 
837  for( i = 0; i < nsubproblems; i++ )
838  {
839  subproblem = SCIPbendersSubproblem(benders, i);
840  if( subproblem != NULL )
841  inf = SCIPinfinity(subproblem);
842  else
843  inf = SCIPsetInfinity(set);
844 
845  SCIPbendersSetSubproblemObjval(benders, i, inf);
846  }
847 }
848 
849 /** compares two Benders' decompositions w. r. to their priority */
850 SCIP_DECL_SORTPTRCOMP(SCIPbendersComp)
851 { /*lint --e{715}*/
852  return ((SCIP_BENDERS*)elem2)->priority - ((SCIP_BENDERS*)elem1)->priority;
853 }
854 
855 /** comparison method for sorting Benders' decompositions w.r.t. to their name */
856 SCIP_DECL_SORTPTRCOMP(SCIPbendersCompName)
857 {
858  return strcmp(SCIPbendersGetName((SCIP_BENDERS*)elem1), SCIPbendersGetName((SCIP_BENDERS*)elem2));
859 }
860 
861 /** method to call, when the priority of a Benders' decomposition was changed */
862 static
863 SCIP_DECL_PARAMCHGD(paramChgdBendersPriority)
864 { /*lint --e{715}*/
865  SCIP_PARAMDATA* paramdata;
866 
867  paramdata = SCIPparamGetData(param);
868  assert(paramdata != NULL);
869 
870  /* use SCIPsetBendersPriority() to mark the Benders' decompositions as unsorted */
871  SCIPsetBendersPriority(scip, (SCIP_BENDERS*)paramdata, SCIPparamGetInt(param)); /*lint !e740*/
872 
873  return SCIP_OKAY;
874 }
875 
876 /** creates a variable mapping between the master problem variables of the source scip and the sub scip */
877 static
879  SCIP_BENDERS* benders, /**< Benders' decomposition of the target SCIP instance */
880  SCIP_SET* sourceset, /**< global SCIP settings from the source SCIP */
881  SCIP_HASHMAP* varmap /**< a hashmap to store the mapping of source variables corresponding
882  * target variables; must not be NULL */
883  )
884 {
885  SCIP_VAR** vars;
886  SCIP_VAR* targetvar;
887  int nvars;
888  int i;
889 
890  assert(benders != NULL);
891  assert(sourceset != NULL);
892  assert(benders->iscopy);
893  assert(benders->mastervarsmap == NULL);
894 
895  /* getting the master problem variable data */
896  vars = SCIPgetVars(sourceset->scip);
897  nvars = SCIPgetNVars(sourceset->scip);
898 
899  /* creating the hashmap for the mapping between the master variable of the target and source scip */
900  SCIP_CALL( SCIPhashmapCreate(&benders->mastervarsmap, SCIPblkmem(sourceset->scip), nvars) );
901 
902  for( i = 0; i < nvars; i++ )
903  {
904  /* getting the variable pointer for the target SCIP variables. The variable mapping returns the target SCIP
905  * varibale for a given source SCIP variable. */
906  targetvar = (SCIP_VAR*) SCIPhashmapGetImage(varmap, vars[i]);
907  if( targetvar != NULL )
908  {
909  SCIP_CALL( SCIPhashmapInsert(benders->mastervarsmap, targetvar, vars[i]) );
910  SCIP_CALL( SCIPcaptureVar(sourceset->scip, vars[i]) );
911  }
912  }
913 
914  return SCIP_OKAY;
915 }
916 
917 /** copies the given Benders' decomposition to a new SCIP */
919  SCIP_BENDERS* benders, /**< Benders' decomposition */
920  SCIP_SET* sourceset, /**< SCIP_SET of SCIP to copy from */
921  SCIP_SET* targetset, /**< SCIP_SET of SCIP to copy to */
922  SCIP_HASHMAP* varmap, /**< a hashmap to store the mapping of source variables corresponding
923  * target variables; if NULL, then the transfer of cuts is not possible */
924  SCIP_Bool threadsafe, /**< must the Benders' decomposition copy be thread safe */
925  SCIP_Bool* valid /**< was the copying process valid? */
926  )
927 {
928  SCIP_BENDERS* targetbenders; /* the copy of the Benders' decomposition struct in the target set */
929  int i;
930 
931  assert(benders != NULL);
932  assert(targetset != NULL);
933  assert(valid != NULL);
934  assert(targetset->scip != NULL);
935 
936  (*valid) = FALSE;
937 
938  if( benders->benderscopy != NULL && targetset->benders_copybenders && SCIPbendersIsActive(benders) )
939  {
940  SCIPsetDebugMsg(targetset, "including Benders' decomposition %s in subscip %p\n", SCIPbendersGetName(benders), (void*)targetset->scip);
941  SCIP_CALL( benders->benderscopy(targetset->scip, benders, threadsafe) );
942 
943  /* copying the Benders' cuts */
944  targetbenders = SCIPsetFindBenders(targetset, SCIPbendersGetName(benders));
945 
946  /* storing the pointer to the source scip instance */
947  targetbenders->sourcescip = sourceset->scip;
948 
949  /* the flag is set to indicate that the Benders' decomposition is a copy */
950  targetbenders->iscopy = TRUE;
951 
952  /* storing whether the lnscheck should be performed */
953  targetbenders->lnscheck = benders->lnscheck;
954  targetbenders->lnsmaxdepth = benders->lnsmaxdepth;
955  targetbenders->lnsmaxcalls = benders->lnsmaxcalls;
956  targetbenders->lnsmaxcallsroot = benders->lnsmaxcallsroot;
957 
958  /* storing whether the Benders' copy required thread safety */
959  targetbenders->threadsafe = threadsafe;
960 
961  /* calling the copy method for the Benders' cuts */
963  for( i = 0; i < benders->nbenderscuts; i++ )
964  {
965  SCIP_CALL( SCIPbenderscutCopyInclude(targetbenders, benders->benderscuts[i], targetset) );
966  }
967 
968  /* When the Benders' decomposition is copied then a variable mapping between the master problem variables is
969  * required. This variable mapping is used to transfer the cuts generated in the target SCIP to the source SCIP.
970  * The variable map is stored in the target Benders' decomposition. This will be freed when the sub-SCIP is freed.
971  */
972  if( varmap != NULL )
973  {
974  SCIP_CALL( createMasterVarMapping(targetbenders, sourceset, varmap) );
975  }
976 
977  assert((varmap != NULL && targetbenders->mastervarsmap != NULL)
978  || (varmap == NULL && targetbenders->mastervarsmap == NULL));
979  }
980 
981  /* if the Benders' decomposition is active, then copy is not valid. */
982  (*valid) = !SCIPbendersIsActive(benders);
983 
984  return SCIP_OKAY;
985 }
986 
987 /** internal method for creating a Benders' decomposition structure */
988 static
990  SCIP_BENDERS** benders, /**< pointer to Benders' decomposition data structure */
991  SCIP_SET* set, /**< global SCIP settings */
992  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
993  BMS_BLKMEM* blkmem, /**< block memory for parameter settings */
994  const char* name, /**< name of Benders' decomposition */
995  const char* desc, /**< description of Benders' decomposition */
996  int priority, /**< priority of the Benders' decomposition */
997  SCIP_Bool cutlp, /**< should Benders' cuts be generated for LP solutions */
998  SCIP_Bool cutpseudo, /**< should Benders' cuts be generated for pseudo solutions */
999  SCIP_Bool cutrelax, /**< should Benders' cuts be generated for relaxation solutions */
1000  SCIP_Bool shareauxvars, /**< should this Benders' use the highest priority Benders aux vars */
1001  SCIP_DECL_BENDERSCOPY ((*benderscopy)), /**< copy method of Benders' decomposition or NULL if you don't want to copy your plugin into sub-SCIPs */
1002  SCIP_DECL_BENDERSFREE ((*bendersfree)), /**< destructor of Benders' decomposition */
1003  SCIP_DECL_BENDERSINIT ((*bendersinit)), /**< initialize Benders' decomposition */
1004  SCIP_DECL_BENDERSEXIT ((*bendersexit)), /**< deinitialize Benders' decomposition */
1005  SCIP_DECL_BENDERSINITPRE((*bendersinitpre)),/**< presolving initialization method for Benders' decomposition */
1006  SCIP_DECL_BENDERSEXITPRE((*bendersexitpre)),/**< presolving deinitialization method for Benders' decomposition */
1007  SCIP_DECL_BENDERSINITSOL((*bendersinitsol)),/**< solving process initialization method of Benders' decomposition */
1008  SCIP_DECL_BENDERSEXITSOL((*bendersexitsol)),/**< solving process deinitialization method of Benders' decomposition */
1009  SCIP_DECL_BENDERSGETVAR((*bendersgetvar)),/**< returns the master variable for a given subproblem variable */
1010  SCIP_DECL_BENDERSCREATESUB((*benderscreatesub)),/**< creates a Benders' decomposition subproblem */
1011  SCIP_DECL_BENDERSPRESUBSOLVE((*benderspresubsolve)),/**< called prior to the subproblem solving loop */
1012  SCIP_DECL_BENDERSSOLVESUBCONVEX((*benderssolvesubconvex)),/**< the solving method for convex Benders' decomposition subproblems */
1013  SCIP_DECL_BENDERSSOLVESUB((*benderssolvesub)),/**< the solving method for the Benders' decomposition subproblems */
1014  SCIP_DECL_BENDERSPOSTSOLVE((*benderspostsolve)),/**< called after the subproblems are solved. */
1015  SCIP_DECL_BENDERSFREESUB((*bendersfreesub)),/**< the freeing method for the Benders' decomposition subproblems */
1016  SCIP_BENDERSDATA* bendersdata /**< Benders' decomposition data */
1017  )
1018 {
1019  char paramname[SCIP_MAXSTRLEN];
1020  char paramdesc[SCIP_MAXSTRLEN];
1021 
1022  assert(benders != NULL);
1023  assert(name != NULL);
1024  assert(desc != NULL);
1025 
1026  /* Checking whether the benderssolvesub and the bendersfreesub are both implemented or both are not implemented */
1027  if( (benderssolvesubconvex == NULL && benderssolvesub == NULL && bendersfreesub != NULL)
1028  || ((benderssolvesubconvex != NULL || benderssolvesub != NULL) && bendersfreesub == NULL) )
1029  {
1030  SCIPerrorMessage("Benders' decomposition <%s> requires that if bendersFreesub%s is implemented, then at least "
1031  "one of bendersSolvesubconvex%s or bendersSolvesub%s are implemented.\n", name, name, name, name);
1032  return SCIP_INVALIDCALL;
1033  }
1034 
1035  SCIP_ALLOC( BMSallocMemory(benders) );
1036  BMSclearMemory(*benders);
1037  SCIP_ALLOC( BMSduplicateMemoryArray(&(*benders)->name, name, strlen(name)+1) );
1038  SCIP_ALLOC( BMSduplicateMemoryArray(&(*benders)->desc, desc, strlen(desc)+1) );
1039  (*benders)->priority = priority;
1040  (*benders)->cutlp = cutlp;
1041  (*benders)->cutpseudo = cutpseudo;
1042  (*benders)->cutrelax = cutrelax;
1043  (*benders)->shareauxvars = shareauxvars;
1044  (*benders)->benderscopy = benderscopy;
1045  (*benders)->bendersfree = bendersfree;
1046  (*benders)->bendersinit = bendersinit;
1047  (*benders)->bendersexit = bendersexit;
1048  (*benders)->bendersinitpre = bendersinitpre;
1049  (*benders)->bendersexitpre = bendersexitpre;
1050  (*benders)->bendersinitsol = bendersinitsol;
1051  (*benders)->bendersexitsol = bendersexitsol;
1052  (*benders)->bendersgetvar = bendersgetvar;
1053  (*benders)->benderscreatesub = benderscreatesub;
1054  (*benders)->benderspresubsolve = benderspresubsolve;
1055  (*benders)->benderssolvesubconvex = benderssolvesubconvex;
1056  (*benders)->benderssolvesub = benderssolvesub;
1057  (*benders)->benderspostsolve = benderspostsolve;
1058  (*benders)->bendersfreesub = bendersfreesub;
1059  (*benders)->bendersdata = bendersdata;
1060  SCIP_CALL( SCIPclockCreate(&(*benders)->setuptime, SCIP_CLOCKTYPE_DEFAULT) );
1061  SCIP_CALL( SCIPclockCreate(&(*benders)->bendersclock, SCIP_CLOCKTYPE_DEFAULT) );
1062 
1063  /* add parameters */
1064  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/priority", name);
1065  (void) SCIPsnprintf(paramdesc, SCIP_MAXSTRLEN, "priority of Benders' decomposition <%s>", name);
1066  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname, paramdesc,
1067  &(*benders)->priority, FALSE, priority, INT_MIN/4, INT_MAX/4,
1068  paramChgdBendersPriority, (SCIP_PARAMDATA*)(*benders)) ); /*lint !e740*/
1069 
1070  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutlp", name);
1071  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1072  "should Benders' cuts be generated for LP solutions?", &(*benders)->cutlp, FALSE, cutlp, NULL, NULL) ); /*lint !e740*/
1073 
1074  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutpseudo", name);
1075  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1076  "should Benders' cuts be generated for pseudo solutions?", &(*benders)->cutpseudo, FALSE, cutpseudo, NULL, NULL) ); /*lint !e740*/
1077 
1078  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutrelax", name);
1079  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1080  "should Benders' cuts be generated for relaxation solutions?", &(*benders)->cutrelax, FALSE, cutrelax, NULL, NULL) ); /*lint !e740*/
1081 
1082  /* These parameters are left for the user to decide in a settings file. This departs from the usual SCIP convention
1083  * where the settings available at the creation of the plugin can be set in the function call.
1084  */
1085  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/transfercuts", name);
1086  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1087  "should Benders' cuts from LNS heuristics be transferred to the main SCIP instance?", &(*benders)->transfercuts,
1088  FALSE, SCIP_DEFAULT_TRANSFERCUTS, NULL, NULL) ); /*lint !e740*/
1089 
1090  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/lnscheck", name);
1091  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1092  "should Benders' decomposition be used in LNS heurisics?", &(*benders)->lnscheck, FALSE, SCIP_DEFAULT_LNSCHECK,
1093  NULL, NULL) ); /*lint !e740*/
1094 
1095  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/lnsmaxdepth", name);
1096  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname,
1097  "maximum depth at which the LNS check is performed (-1: no limit)", &(*benders)->lnsmaxdepth, TRUE,
1099 
1100  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/lnsmaxcalls", name);
1101  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname,
1102  "the maximum number of Benders' decomposition calls in LNS heuristics (-1: no limit)", &(*benders)->lnsmaxcalls,
1103  TRUE, SCIP_DEFAULT_LNSMAXCALLS, -1, INT_MAX, NULL, NULL) );
1104 
1105  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/lnsmaxcallsroot", name);
1106  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname,
1107  "the maximum number of root node Benders' decomposition calls in LNS heuristics (-1: no limit)",
1108  &(*benders)->lnsmaxcallsroot, TRUE, SCIP_DEFAULT_LNSMAXCALLSROOT, -1, INT_MAX, NULL, NULL) );
1109 
1110  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutsasconss", name);
1111  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1112  "should the transferred cuts be added as constraints?", &(*benders)->cutsasconss, FALSE,
1113  SCIP_DEFAULT_CUTSASCONSS, NULL, NULL) ); /*lint !e740*/
1114 
1115  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/subprobfrac", name);
1116  SCIP_CALL( SCIPsetAddRealParam(set, messagehdlr, blkmem, paramname,
1117  "fraction of subproblems that are solved in each iteration", &(*benders)->subprobfrac, FALSE,
1118  SCIP_DEFAULT_SUBPROBFRAC, 0.0, 1.0, NULL, NULL) ); /*lint !e740*/
1119 
1120  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/updateauxvarbound", name);
1121  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1122  "should the auxiliary variable bound be updated by solving the subproblem?", &(*benders)->updateauxvarbound,
1123  FALSE, SCIP_DEFAULT_UPDATEAUXVARBOUND, NULL, NULL) ); /*lint !e740*/
1124 
1125  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/auxvarsimplint", name);
1126  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1127  "if the subproblem objective is integer, then define the auxiliary variables as implied integers?",
1128  &(*benders)->auxvarsimplint, FALSE, SCIP_DEFAULT_AUXVARSIMPLINT, NULL, NULL) ); /*lint !e740*/
1129 
1130  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutcheck", name);
1131  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1132  "should Benders' cuts be generated while checking solutions?",
1133  &(*benders)->cutcheck, FALSE, SCIP_DEFAULT_CUTCHECK, NULL, NULL) ); /*lint !e740*/
1134 
1135  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutstrengthenmult", name);
1136  SCIP_CALL( SCIPsetAddRealParam(set, messagehdlr, blkmem, paramname,
1137  "the convex combination multiplier for the cut strengthening", &(*benders)->convexmult, FALSE,
1138  SCIP_DEFAULT_STRENGTHENMULT, 0.0, 1.0, NULL, NULL) ); /*lint !e740*/
1139 
1140  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/noimprovelimit", name);
1141  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname,
1142  "the maximum number of cut strengthening without improvement", &(*benders)->noimprovelimit, TRUE,
1143  SCIP_DEFAULT_NOIMPROVELIMIT, 0, INT_MAX, NULL, NULL) );
1144 
1145  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/corepointperturb", name);
1146  SCIP_CALL( SCIPsetAddRealParam(set, messagehdlr, blkmem, paramname,
1147  "the constant use to perturb the cut strengthening core point", &(*benders)->perturbeps, FALSE,
1148  SCIP_DEFAULT_STRENGTHENPERTURB, 0.0, 1.0, NULL, NULL) ); /*lint !e740*/
1149 
1150  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutstrengthenenabled", name);
1151  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1152  "should the core point cut strengthening be employed (only applied to fractional solutions or continuous subproblems)?",
1153  &(*benders)->strengthenenabled, FALSE, SCIP_DEFAULT_STRENGTHENENABLED, NULL, NULL) ); /*lint !e740*/
1154 
1155  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/cutstrengthenintpoint", name);
1156  SCIP_CALL( SCIPsetAddCharParam(set, messagehdlr, blkmem, paramname,
1157  "where should the strengthening interior point be sourced from ('l'p relaxation, 'f'irst solution, 'i'ncumbent solution, 'r'elative interior point, vector of 'o'nes, vector of 'z'eros)",
1158  &(*benders)->strengthenintpoint, FALSE, SCIP_DEFAULT_STRENGTHENINTPOINT, "lfiroz", NULL, NULL) ); /*lint !e740*/
1159 
1160  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/numthreads", name);
1161  SCIP_CALL( SCIPsetAddIntParam(set, messagehdlr, blkmem, paramname,
1162  "the number of threads to use when solving the subproblems", &(*benders)->numthreads, TRUE,
1163  SCIP_DEFAULT_NUMTHREADS, 1, INT_MAX, NULL, NULL) );
1164 
1165  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/execfeasphase", name);
1166  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1167  "should a feasibility phase be executed during the root node, i.e. adding slack variables to constraints to ensure feasibility",
1168  &(*benders)->execfeasphase, FALSE, SCIP_DEFAULT_EXECFEASPHASE, NULL, NULL) ); /*lint !e740*/
1169 
1170  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/slackvarcoef", name);
1171  SCIP_CALL( SCIPsetAddRealParam(set, messagehdlr, blkmem, paramname,
1172  "the initial objective coefficient of the slack variables in the subproblem", &(*benders)->slackvarcoef, FALSE,
1173  SCIP_DEFAULT_SLACKVARCOEF, 0.0, SCIPsetInfinity(set), NULL, NULL) ); /*lint !e740*/
1174 
1175  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/maxslackvarcoef", name);
1176  SCIP_CALL( SCIPsetAddRealParam(set, messagehdlr, blkmem, paramname,
1177  "the maximal objective coefficient of the slack variables in the subproblem", &(*benders)->maxslackvarcoef, FALSE,
1178  SCIP_DEFAULT_MAXSLACKVARCOEF, 0.0, SCIPsetInfinity(set), NULL, NULL) ); /*lint !e740*/
1179 
1180  (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/checkconsconvexity", name);
1181  SCIP_CALL( SCIPsetAddBoolParam(set, messagehdlr, blkmem, paramname,
1182  "should the constraints of the subproblems be checked for convexity?", &(*benders)->checkconsconvexity, FALSE,
1183  SCIP_DEFAULT_CHECKCONSCONVEXITY, NULL, NULL) ); /*lint !e740*/
1184 
1185  return SCIP_OKAY;
1186 }
1187 
1188 /** creates a Benders' decomposition structure
1189  *
1190  * To use the Benders' decomposition for solving a problem, it first has to be activated with a call to SCIPactivateBenders().
1191  */
1193  SCIP_BENDERS** benders, /**< pointer to Benders' decomposition data structure */
1194  SCIP_SET* set, /**< global SCIP settings */
1195  SCIP_MESSAGEHDLR* messagehdlr, /**< message handler */
1196  BMS_BLKMEM* blkmem, /**< block memory for parameter settings */
1197  const char* name, /**< name of Benders' decomposition */
1198  const char* desc, /**< description of Benders' decomposition */
1199  int priority, /**< priority of the Benders' decomposition */
1200  SCIP_Bool cutlp, /**< should Benders' cuts be generated for LP solutions */
1201  SCIP_Bool cutpseudo, /**< should Benders' cuts be generated for pseudo solutions */
1202  SCIP_Bool cutrelax, /**< should Benders' cuts be generated for relaxation solutions */
1203  SCIP_Bool shareauxvars, /**< should this Benders' use the highest priority Benders aux vars */
1204  SCIP_DECL_BENDERSCOPY ((*benderscopy)), /**< copy method of Benders' decomposition or NULL if you don't want to copy your plugin into sub-SCIPs */
1205  SCIP_DECL_BENDERSFREE ((*bendersfree)), /**< destructor of Benders' decomposition */
1206  SCIP_DECL_BENDERSINIT ((*bendersinit)), /**< initialize Benders' decomposition */
1207  SCIP_DECL_BENDERSEXIT ((*bendersexit)), /**< deinitialize Benders' decomposition */
1208  SCIP_DECL_BENDERSINITPRE((*bendersinitpre)),/**< presolving initialization method for Benders' decomposition */
1209  SCIP_DECL_BENDERSEXITPRE((*bendersexitpre)),/**< presolving deinitialization method for Benders' decomposition */
1210  SCIP_DECL_BENDERSINITSOL((*bendersinitsol)),/**< solving process initialization method of Benders' decomposition */
1211  SCIP_DECL_BENDERSEXITSOL((*bendersexitsol)),/**< solving process deinitialization method of Benders' decomposition */
1212  SCIP_DECL_BENDERSGETVAR((*bendersgetvar)),/**< returns the master variable for a given subproblem variable */
1213  SCIP_DECL_BENDERSCREATESUB((*benderscreatesub)),/**< creates a Benders' decomposition subproblem */
1214  SCIP_DECL_BENDERSPRESUBSOLVE((*benderspresubsolve)),/**< called prior to the subproblem solving loop */
1215  SCIP_DECL_BENDERSSOLVESUBCONVEX((*benderssolvesubconvex)),/**< the solving method for convex Benders' decomposition subproblems */
1216  SCIP_DECL_BENDERSSOLVESUB((*benderssolvesub)),/**< the solving method for the Benders' decomposition subproblems */
1217  SCIP_DECL_BENDERSPOSTSOLVE((*benderspostsolve)),/**< called after the subproblems are solved. */
1218  SCIP_DECL_BENDERSFREESUB((*bendersfreesub)),/**< the freeing method for the Benders' decomposition subproblems */
1219  SCIP_BENDERSDATA* bendersdata /**< Benders' decomposition data */
1220  )
1221 {
1222  assert(benders != NULL);
1223  assert(name != NULL);
1224  assert(desc != NULL);
1225 
1226  SCIP_CALL_FINALLY( doBendersCreate(benders, set, messagehdlr, blkmem, name, desc, priority, cutlp, cutpseudo,
1227  cutrelax, shareauxvars, benderscopy, bendersfree, bendersinit, bendersexit, bendersinitpre, bendersexitpre,
1228  bendersinitsol, bendersexitsol, bendersgetvar, benderscreatesub, benderspresubsolve, benderssolvesubconvex,
1229  benderssolvesub, benderspostsolve, bendersfreesub, bendersdata), (void) SCIPbendersFree(benders, set) );
1230 
1231  return SCIP_OKAY;
1232 }
1233 
1234 
1235 /** releases the variables that have been captured in the hashmap */
1236 static
1238  SCIP* scip, /**< the SCIP data structure */
1239  SCIP_BENDERS* benders /**< Benders' decomposition */
1240  )
1241 {
1242  int nentries;
1243  int i;
1244 
1245  assert(scip != NULL);
1246  assert(benders != NULL);
1247 
1248  assert(benders->mastervarsmap != NULL);
1249 
1250  nentries = SCIPhashmapGetNEntries(benders->mastervarsmap);
1251 
1252  for( i = 0; i < nentries; ++i )
1253  {
1254  SCIP_HASHMAPENTRY* entry;
1255  entry = SCIPhashmapGetEntry(benders->mastervarsmap, i);
1256 
1257  if( entry != NULL )
1258  {
1259  SCIP_VAR* var;
1260  var = (SCIP_VAR*) SCIPhashmapEntryGetImage(entry);
1261 
1262  SCIP_CALL( SCIPreleaseVar(scip, &var) );
1263  }
1264  }
1265 
1266  return SCIP_OKAY;
1267 }
1268 
1269 
1270 /** calls destructor and frees memory of Benders' decomposition */
1272  SCIP_BENDERS** benders, /**< pointer to Benders' decomposition data structure */
1273  SCIP_SET* set /**< global SCIP settings */
1274  )
1275 {
1276  int i;
1277 
1278  assert(benders != NULL);
1279  assert(*benders != NULL);
1280  assert(!(*benders)->initialized);
1281  assert(set != NULL);
1282 
1283  /* call destructor of Benders' decomposition */
1284  if( (*benders)->bendersfree != NULL )
1285  {
1286  SCIP_CALL( (*benders)->bendersfree(set->scip, *benders) );
1287  }
1288 
1289  /* if the Benders' decomposition is a copy and a varmap has been passed to SCIP_BENDERS, then the variable map
1290  * between the source and the target SCIP needs to be freed.
1291  */
1292  if( (*benders)->iscopy && (*benders)->mastervarsmap != NULL )
1293  {
1294  SCIP_CALL( releaseVarMappingHashmapVars((*benders)->sourcescip, (*benders)) );
1295  SCIPhashmapFree(&(*benders)->mastervarsmap);
1296  }
1297 
1298  /* freeing the Benders' cuts */
1299  for( i = 0; i < (*benders)->nbenderscuts; i++ )
1300  {
1301  SCIP_CALL( SCIPbenderscutFree(&((*benders)->benderscuts[i]), set) );
1302  }
1303  BMSfreeMemoryArrayNull(&(*benders)->benderscuts);
1304 
1305  SCIPclockFree(&(*benders)->bendersclock);
1306  SCIPclockFree(&(*benders)->setuptime);
1307  BMSfreeMemoryArray(&(*benders)->name);
1308  BMSfreeMemoryArray(&(*benders)->desc);
1309  BMSfreeMemory(benders);
1310 
1311  return SCIP_OKAY;
1312 }
1313 
1314 /* adds a slack variable to the given constraint */
1315 static
1317  SCIP* scip, /**< the SCIP data structure */
1318  SCIP_BENDERS* benders, /**< Benders' decomposition */
1319  SCIP_CONS* cons, /**< constraint to which the slack variable(s) is added to */
1320  SCIP_CONSHDLR** linearconshdlrs, /**< an array storing the linear constraint handlers */
1321  SCIP_CONSHDLR* nlconshdlr, /**< pointer to the nonlinear constraint handler */
1322  int nlinearconshdlrs /**< the number of linear constraint handlers */
1323  )
1324 {
1325  SCIP_CONSHDLR* conshdlr;
1326  SCIP_VAR* var;
1327  SCIP_Real rhs;
1328  SCIP_Real lhs;
1329  SCIP_Real objcoef;
1330  int i;
1331  SCIP_Bool linearcons;
1332  SCIP_Bool success;
1333  char name[SCIP_MAXSTRLEN];
1334 
1335  conshdlr = SCIPconsGetHdlr(cons);
1336 
1337  /* assume that the constraint is not linear, then we check whether it is linear */
1338  linearcons = FALSE;
1339 
1340  /* checking whether the constraint is a linear constraint. If so, we add a coefficient to the constraint */
1341  for( i = 0; i < nlinearconshdlrs; ++i )
1342  {
1343  if( conshdlr == linearconshdlrs[i] )
1344  {
1345  linearcons = TRUE;
1346  break;
1347  }
1348  }
1349 
1350  if( !linearcons && conshdlr != nlconshdlr )
1351  {
1352  SCIPwarningMessage(scip, "The subproblem includes constraint <%s>. "
1353  "This is not supported and the slack variable will not be added to the constraint. Feasibility cuts may be invalid.\n",
1354  SCIPconshdlrGetName(conshdlr));
1355  }
1356 
1357  if( linearcons )
1358  {
1359  rhs = SCIPconsGetRhs(scip, cons, &success);
1360  assert(success);
1361  lhs = SCIPconsGetLhs(scip, cons, &success);
1362  assert(success);
1363  }
1364  else
1365  {
1366  rhs = SCIPgetRhsNonlinear(cons);
1367  lhs = SCIPgetLhsNonlinear(cons);
1368  }
1369 
1370  /* getting the objective coefficient for the slack variables */
1371  objcoef = benders->slackvarcoef;
1372 
1373  /* if the right hand side is finite, then we need to add a slack variable with a negative coefficient */
1374  if( !SCIPisInfinity(scip, rhs) )
1375  {
1376  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_%s_neg", SLACKVAR_NAME, SCIPconsGetName(cons) );
1377 
1378  SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, SCIPinfinity(scip), objcoef, SCIP_VARTYPE_CONTINUOUS) );
1379 
1380  /* adding the slack variable to the subproblem */
1381  SCIP_CALL( SCIPaddVar(scip, var) );
1382 
1383  /* adds the slack variable to the constraint */
1384  if( linearcons )
1385  {
1386  SCIP_CALL( SCIPconsAddCoef(scip, cons, var, -1.0) );
1387  }
1388  else
1389  {
1390  SCIP_CALL( SCIPaddLinearVarNonlinear(scip, cons, var, -1.0) );
1391  }
1392 
1393  /* releasing the variable */
1394  SCIP_CALL( SCIPreleaseVar(scip, &var) );
1395  }
1396 
1397  /* if the left hand side if finite, then we need to add a slack variable with a positive coefficient */
1398  if( !SCIPisInfinity(scip, -lhs) )
1399  {
1400  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_%s_pos", SLACKVAR_NAME, SCIPconsGetName(cons) );
1401 
1402  SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, SCIPinfinity(scip), objcoef, SCIP_VARTYPE_CONTINUOUS) );
1403 
1404  /* adding the slack variable to the subproblem */
1405  SCIP_CALL( SCIPaddVar(scip, var) );
1406 
1407  /* adds the slack variable to the constraint */
1408  if( linearcons )
1409  {
1410  SCIP_CALL( SCIPconsAddCoef(scip, cons, var, 1.0) );
1411  }
1412  else
1413  {
1414  SCIP_CALL( SCIPaddLinearVarNonlinear(scip, cons, var, 1.0) );
1415  }
1416 
1417  /* releasing the variable */
1418  SCIP_CALL( SCIPreleaseVar(scip, &var) );
1419  }
1420 
1421  return SCIP_OKAY;
1422 }
1423 
1424 /** adds the slack variables to each of the constraints for the generation of feasibility cuts for the given non-linear
1425  * subproblem
1426  */
1427 static
1429  SCIP_BENDERS* benders, /**< Benders' decomposition */
1430  SCIP_SET* set, /**< global SCIP settings */
1431  int probnumber /**< the subproblem number */
1432  )
1433 {
1434  SCIP* subproblem;
1435  SCIP_CONSHDLR* linearconshdlrs[NLINEARCONSHDLRS];
1436  SCIP_CONSHDLR* nlconshdlr;
1437  SCIP_CONS* cons;
1438  int i;
1439 
1440  assert(benders != NULL);
1441  assert(set != NULL);
1442  assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
1443 
1444  subproblem = SCIPbendersSubproblem(benders, probnumber);
1445 
1446  /* get pointers to linear constraints handlers, so can avoid string comparisons */
1447  linearconshdlrs[0] = SCIPfindConshdlr(subproblem, "knapsack");
1448  linearconshdlrs[1] = SCIPfindConshdlr(subproblem, "linear");
1449  linearconshdlrs[2] = SCIPfindConshdlr(subproblem, "logicor");
1450  linearconshdlrs[3] = SCIPfindConshdlr(subproblem, "setppc");
1451  linearconshdlrs[4] = SCIPfindConshdlr(subproblem, "varbound");
1452 
1453  nlconshdlr = SCIPfindConshdlr(subproblem, "nonlinear");
1454 
1455  for( i = 0; i < SCIPgetNOrigConss(subproblem); ++i )
1456  {
1457  cons = SCIPgetOrigConss(subproblem)[i];
1458 
1459  /* adding the slack variables to the constraint */
1460  SCIP_CALL( addSlackVars(subproblem, benders, cons, linearconshdlrs, nlconshdlr, NLINEARCONSHDLRS) );
1461  }
1462 
1463  return SCIP_OKAY;
1464 }
1465 
1466 /** initialises a MIP subproblem by putting the problem into SCIP_STAGE_SOLVING. This is achieved by calling SCIPsolve
1467  * and then interrupting the solve in a node focus event handler.
1468  * The LP subproblem is also initialised using this method; however, a different event handler is added. This event
1469  * handler will put the LP subproblem into probing mode.
1470  * The MIP solving function is called to initialise the subproblem because this function calls SCIPsolve with the
1471  * appropriate parameter settings for Benders' decomposition.
1472  */
1473 static
1475  SCIP_BENDERS* benders, /**< Benders' decomposition */
1476  SCIP_SET* set, /**< global SCIP settings */
1477  int probnumber, /**< the subproblem number */
1478  SCIP_Bool* success /**< was the initialisation process successful */
1479  )
1480 {
1481  SCIP* subproblem;
1482  SCIP_STATUS solvestatus;
1483  SCIP_Bool cutoff;
1484 
1485  assert(benders != NULL);
1486  assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
1487  assert(success != NULL);
1488 
1489  (*success) = FALSE;
1490 
1491  subproblem = SCIPbendersSubproblem(benders, probnumber);
1492  assert(subproblem != NULL);
1493 
1494  /* Getting the problem into the right SCIP stage for solving */
1495  SCIP_CALL( SCIPbendersSolveSubproblemCIP(set->scip, benders, probnumber, &solvestatus, FALSE) );
1496 
1497  /* Constructing the LP that can be solved in later iterations */
1498  if( solvestatus != SCIP_STATUS_BESTSOLLIMIT && solvestatus != SCIP_STATUS_TIMELIMIT
1499  && solvestatus != SCIP_STATUS_MEMLIMIT )
1500  {
1501  assert(SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING);
1502 
1503  SCIP_CALL( SCIPconstructLP(subproblem, &cutoff) );
1504  (*success) = TRUE;
1505  }
1506 
1507  return SCIP_OKAY;
1508 }
1509 
1510 
1511 /** initialises an LP subproblem by putting the problem into probing mode. The probing mode is invoked in a node focus
1512  * event handler. This event handler is added just prior to calling the initialise subproblem function.
1513  */
1514 static
1516  SCIP_BENDERS* benders, /**< Benders' decomposition */
1517  SCIP_SET* set, /**< global SCIP settings */
1518  int probnumber /**< the subproblem number */
1519  )
1520 {
1521  SCIP* subproblem;
1522  SCIP_EVENTHDLR* eventhdlr;
1523  SCIP_EVENTHDLRDATA* eventhdlrdata;
1524  SCIP_Bool success;
1525 
1526  assert(benders != NULL);
1527  assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
1528 
1529  subproblem = SCIPbendersSubproblem(benders, probnumber);
1530  assert(subproblem != NULL);
1531 
1532  /* include event handler into SCIP */
1533  SCIP_CALL( SCIPallocBlockMemory(subproblem, &eventhdlrdata) );
1534 
1535  SCIP_CALL( initEventhandlerData(subproblem, eventhdlrdata) );
1536 
1538  eventExecBendersNodefocus, eventhdlrdata) );
1539  SCIP_CALL( SCIPsetEventhdlrInitsol(subproblem, eventhdlr, eventInitsolBendersNodefocus) );
1540  SCIP_CALL( SCIPsetEventhdlrExitsol(subproblem, eventhdlr, eventExitsolBendersNodefocus) );
1541  SCIP_CALL( SCIPsetEventhdlrExit(subproblem, eventhdlr, eventExitBendersNodefocus) );
1542  SCIP_CALL( SCIPsetEventhdlrFree(subproblem, eventhdlr, eventFreeBendersNodefocus) );
1543  assert(eventhdlr != NULL);
1544 
1545  /* calling an initial solve to put the problem into probing mode */
1546  SCIP_CALL( initialiseSubproblem(benders, set, probnumber, &success) );
1547 
1548  return SCIP_OKAY; /*lint !e438*/
1549 }
1550 
1551 /** checks whether the convex relaxation of the subproblem is sufficient to solve the original problem to optimality
1552  *
1553  * We check whether we can conclude that the CIP is actually an LP or a convex NLP.
1554  * To do this, we check that all variables are of continuous type and that every constraint is either handled by known
1555  * linear constraint handler (knapsack, linear, logicor, setppc, varbound) or the nonlinear constraint handler.
1556  * In the latter case, we also check whether the nonlinear constraint is convex.
1557  * Further, nonlinear constraints are only considered if an NLP solver interface is available, i.e., and NLP could
1558  * be solved.
1559  * If constraints are present that cannot be identified as linear or convex nonlinear, then we assume that the
1560  * problem is not convex, thus solving its LP or NLP relaxation will not be sufficient.
1561  */
1562 static
1564  SCIP_BENDERS* benders, /**< Benders' decomposition */
1565  SCIP_SET* set, /**< global SCIP settings */
1566  int probnumber /**< the subproblem number, or -1 for the master problem */
1567  )
1568 {
1569  SCIP* subproblem;
1570  SCIP_CONSHDLR* conshdlr;
1571  SCIP_CONS* cons;
1572  SCIP_HASHMAP* assumevarfixed;
1573  SCIP_VAR** vars;
1574  int nvars;
1575  int nbinvars;
1576  int nintvars;
1577  int nimplintvars;
1578  int i;
1579  int j;
1580  SCIP_Bool convexcons;
1581  SCIP_Bool discretevar;
1582  SCIP_Bool isnonlinear;
1583  SCIP_CONSHDLR* linearconshdlrs[NLINEARCONSHDLRS];
1584  SCIP_CONSHDLR* nlconshdlr = NULL;
1585 
1586  assert(benders != NULL);
1587  assert(set != NULL);
1588  assert(probnumber >= -1 && probnumber < SCIPbendersGetNSubproblems(benders));
1589 
1590  assumevarfixed = NULL;
1591  if( probnumber >= 0 )
1592  subproblem = SCIPbendersSubproblem(benders, probnumber);
1593  else
1594  subproblem = set->scip;
1595 
1596  assert(subproblem != NULL);
1597 
1598  convexcons = FALSE;
1599  discretevar = FALSE;
1600  isnonlinear = FALSE;
1601 
1602  /* getting the number of integer and binary variables to determine the problem type */
1603  SCIP_CALL( SCIPgetVarsData(subproblem, &vars, &nvars, &nbinvars, &nintvars, &nimplintvars, NULL) );
1604 
1605  /* if there are any binary, integer or implied integer variables, then the subproblems is marked as non-convex */
1606  if( nbinvars != 0 || nintvars != 0 || nimplintvars != 0 )
1607  {
1608  discretevar = TRUE;
1609  }
1610 
1611  /* get pointers to linear constraints handlers, so can avoid string comparisons */
1612  linearconshdlrs[0] = SCIPfindConshdlr(subproblem, "knapsack");
1613  linearconshdlrs[1] = SCIPfindConshdlr(subproblem, "linear");
1614  linearconshdlrs[2] = SCIPfindConshdlr(subproblem, "logicor");
1615  linearconshdlrs[3] = SCIPfindConshdlr(subproblem, "setppc");
1616  linearconshdlrs[4] = SCIPfindConshdlr(subproblem, "varbound");
1617 
1618  /* Get pointer to the nonlinear constraint handler if we also have an NLP solver to solve NLPs.
1619  * If there is no NLP solver, but there are (convex) nonlinear constraints, then the LP relaxation of subproblems
1620  * will (currently) not be sufficient to solve subproblems to optimality. Thus, we also take the presence of convex
1621  * nonlinear constraints as signal for having to solve the CIP eventually, thus, by abuse of notation,
1622  * return not-convex here. In summary, we do not need to have a special look onto non-linear constraints
1623  * if no NLP solver is present, and can treat them as any other constraint that is not of linear type.
1624  */
1625  if( SCIPgetNNlpis(subproblem) > 0 )
1626  {
1627  nlconshdlr = SCIPfindConshdlr(subproblem, "nonlinear");
1628  }
1629 
1630  /* if the nonlinear constraint handler exists, then we create a hashmap of variables that can be assumed to be fixed.
1631  * These variables correspond to the copies of the master variables in the subproblem
1632  */
1633  if( probnumber >= 0 && nlconshdlr != NULL )
1634  {
1635  SCIP_VAR* mappedvar;
1636 
1637  SCIP_CALL( SCIPhashmapCreate(&assumevarfixed, SCIPblkmem(set->scip), SCIPgetNVars(subproblem)) );
1638 
1639  /* finding the subproblem variables that correspond to master variables */
1640  for( i = 0; i < nvars; i++ )
1641  {
1642  /* getting the corresponding master problem variable for the given variable */
1643  SCIP_CALL( SCIPbendersGetVar(benders, set, vars[i], &mappedvar, -1) );
1644 
1645  /* if the mapped variable is not NULL, then it must be stored as a possible fixed variable */
1646  if( mappedvar != NULL )
1647  {
1648  SCIP_CALL( SCIPhashmapInsert(assumevarfixed, vars[i], vars[i]) );
1649  }
1650  }
1651  }
1652 
1653  for( i = 0; i < SCIPgetNOrigConss(subproblem); ++i )
1654  {
1655  cons = SCIPgetOrigConss(subproblem)[i];
1656  conshdlr = SCIPconsGetHdlr(cons);
1657 
1658  for( j = 0; j < NLINEARCONSHDLRS; ++j )
1659  if( conshdlr == linearconshdlrs[j] )
1660  break;
1661 
1662  /* if linear constraint, then we are good */
1663  if( j < NLINEARCONSHDLRS )
1664  {
1665 #ifdef SCIP_MOREDEBUG
1666  SCIPdebugMsg(subproblem, "subproblem <%s>: constraint <%s> is linear\n", SCIPgetProbName(subproblem), SCIPconsGetName(cons));
1667 #endif
1668  continue;
1669  }
1670 
1671  /* if cons_nonlinear (and nlconshdlr != NULL), then check whether convex */
1672  if( conshdlr == nlconshdlr )
1673  {
1674  SCIP_Bool isconvex;
1675  SCIP_EXPRCURV curv;
1676  SCIP_Bool havelhs;
1677  SCIP_Bool haverhs;
1678 
1679  isnonlinear = TRUE;
1680 
1681  havelhs = !SCIPisInfinity(subproblem, -SCIPgetLhsNonlinear(cons));
1682  haverhs = !SCIPisInfinity(subproblem, SCIPgetRhsNonlinear(cons));
1683  if( havelhs && haverhs )
1684  {
1685  isconvex = FALSE;
1686  }
1687  else
1688  {
1689  /* look at curvature stored in cons, though at this stage this will be unknown a.a. */
1690  curv = SCIPgetCurvatureNonlinear(cons);
1691  isconvex = ((!havelhs || (curv & SCIP_EXPRCURV_CONCAVE) == SCIP_EXPRCURV_CONCAVE)) &&
1692  ((!haverhs || (curv & SCIP_EXPRCURV_CONVEX) == SCIP_EXPRCURV_CONVEX));
1693 
1694  if( !isconvex )
1695  {
1696  /* if not found convex, compute curvature via nlhdlr_convex and decide again */
1697 
1698  /* make sure activities are up to date. SCIPhasExprCurvature currently assumes that this is already the case */
1699  SCIP_CALL( SCIPevalExprActivity(subproblem, SCIPgetExprNonlinear(cons)) );
1700 
1701  SCIP_CALL( SCIPhasExprCurvature(subproblem, SCIPgetExprNonlinear(cons), havelhs ? SCIP_EXPRCURV_CONCAVE : SCIP_EXPRCURV_CONVEX, &isconvex, assumevarfixed) );
1702  }
1703  }
1704 
1705  if( isconvex )
1706  {
1707 #ifdef SCIP_MOREDEBUG
1708  SCIPdebugMsg(subproblem, "subproblem <%s>: nonlinear constraint <%s> is convex\n", SCIPgetProbName(subproblem), SCIPconsGetName(cons));
1709 #endif
1710  continue;
1711  }
1712  else
1713  {
1714 #ifdef SCIP_MOREDEBUG
1715  SCIPdebugMsg(subproblem, "subproblem <%s>: nonlinear constraint <%s> not convex\n", SCIPgetProbName(subproblem), SCIPconsGetName(cons));
1716 #endif
1717  goto TERMINATE;
1718  }
1719  }
1720 
1721 #ifdef SCIP_MOREDEBUG
1722  SCIPdebugMsg(subproblem, "subproblem <%s>: potentially nonconvex constraint <%s>\n", SCIPgetProbName(subproblem), SCIPconsGetName(cons));
1723 #endif
1724  goto TERMINATE;
1725  }
1726 
1727  /* if we made it until here, then all constraints are known and convex */
1728  convexcons = TRUE;
1729 
1730 TERMINATE:
1731  /* setting the flag for the convexity of the subproblem. If convexity doesn't need to be checked, then it is assumed
1732  * that the subproblems are convex. However, if there are discrete variables, then the problem must be set as
1733  * non-convex. The discrete master variables will be changed to continuous, but this will happen at the first call to
1734  * SCIPbendersSetupSubproblem
1735  */
1736  if( probnumber >= 0 )
1737  {
1738  convexcons = convexcons || !benders->checkconsconvexity;
1739 
1740  if( convexcons && !discretevar )
1742  else if( convexcons && discretevar )
1744  else if( !convexcons && !discretevar )
1746  else if( !convexcons && discretevar )
1748  else
1749  SCIPABORT();
1750  }
1751 
1752  /* setting the non-linear subproblem flag */
1753  if( probnumber >= 0 )
1754  SCIPbendersSetSubproblemIsNonlinear(benders, probnumber, isnonlinear);
1755  else
1756  SCIPbendersSetMasterIsNonlinear(benders, isnonlinear);
1757 
1758  if( probnumber >= 0 )
1759  {
1760  SCIPsetDebugMsg(set, "subproblem <%s> has been found to be of type %d\n", SCIPgetProbName(subproblem),
1761  SCIPbendersGetSubproblemType(benders, probnumber));
1762  }
1763 
1764  /* releasing the fixed variable hashmap */
1765  if( assumevarfixed != NULL )
1766  SCIPhashmapFree(&assumevarfixed);
1767 
1768  return SCIP_OKAY;
1769 }
1770 
1771 /** creates the subproblems and registers it with the Benders' decomposition struct */
1772 static
1774  SCIP_BENDERS* benders, /**< Benders' decomposition */
1775  SCIP_SET* set /**< global SCIP settings */
1776  )
1777 {
1778  SCIP* subproblem;
1779  SCIP_EVENTHDLR* eventhdlr;
1780  SCIP_VAR* mastervar;
1781  SCIP_VAR** vars;
1782  int nvars;
1783  int nsubproblems;
1784  int i;
1785  int j;
1786 
1787  assert(benders != NULL);
1788  assert(set != NULL);
1789 
1790  /* if the subproblems have already been created, then they will not be created again. This is the case if the
1791  * transformed problem has been freed and then retransformed. The subproblems should only be created when the problem
1792  * is first transformed. */
1793  if( benders->subprobscreated )
1794  return SCIP_OKAY;
1795 
1796  nsubproblems = SCIPbendersGetNSubproblems(benders);
1797 
1798  /* creating all subproblems */
1799  for( i = 0; i < nsubproblems; i++ )
1800  {
1801  /* calling the create subproblem call back method */
1802  SCIP_CALL( benders->benderscreatesub(set->scip, benders, i) );
1803 
1804  subproblem = SCIPbendersSubproblem(benders, i);
1805 
1806  /* the subproblem SCIP instance could be set to NULL. This is because user defined subproblem solving methods
1807  * could be used that don't solve a SCIP instance. Thus, the following setup of the subproblem SCIP instance is
1808  * not required.
1809  *
1810  * NOTE: since the subproblems are supplied as NULL pointers, the internal convexity check can not be performed.
1811  * The user needs to explicitly specify the subproblem type.
1812  */
1813  if( subproblem != NULL )
1814  {
1815  /* setting global limits for the subproblems. This overwrites the limits set by the user */
1816  SCIP_CALL( SCIPsetIntParam(subproblem, "limits/maxorigsol", 0) );
1817 
1818  /* getting the number of integer and binary variables to determine the problem type */
1819  SCIP_CALL( SCIPgetVarsData(subproblem, &vars, &nvars, NULL, NULL, NULL, NULL) );
1820 
1821  /* The objective function coefficients of the master problem are set to zero. This is necessary for the Benders'
1822  * decomposition algorithm, since the cut methods and the objective function check assumes that the objective
1823  * coefficients of the master problem variables are zero.
1824  *
1825  * This only occurs if the Benders' decomposition is not a copy. It is assumed that the correct objective
1826  * coefficients are given during the first subproblem creation.
1827  *
1828  * If the subproblems were copied, then the master variables will be checked to ensure that they have a zero
1829  * objective value.
1830  */
1831  if( !benders->iscopy || benders->threadsafe )
1832  {
1833  SCIP_Bool objchanged = FALSE;
1834 
1835  assert(SCIPgetStage(subproblem) == SCIP_STAGE_PROBLEM);
1836  for( j = 0; j < nvars; j++ )
1837  {
1838  /* retrieving the master problem variable */
1839  SCIP_CALL( SCIPbendersGetVar(benders, set, vars[j], &mastervar, -1) );
1840 
1841  /* if mastervar is not NULL, then the subproblem variable has a corresponding master problem variable */
1842  if( mastervar != NULL && !SCIPisZero(subproblem, SCIPvarGetObj(vars[j])) )
1843  {
1844  SCIPverbMessage(subproblem, SCIP_VERBLEVEL_FULL, NULL, "Benders' decomposition: Changing the objective "
1845  "coefficient of copy of master problem variable <%s> in subproblem %d to zero.\n",
1846  SCIPvarGetName(mastervar), i);
1847  /* changing the subproblem variable objective coefficient to zero */
1848  SCIP_CALL( SCIPchgVarObj(subproblem, vars[j], 0.0) );
1849 
1850  objchanged = TRUE;
1851  }
1852  }
1853 
1854  if( objchanged )
1855  {
1856  SCIPverbMessage(subproblem, SCIP_VERBLEVEL_HIGH, NULL, "Benders' decomposition: Objective coefficients of "
1857  "copied of master problem variables has been changed to zero.\n");
1858  }
1859  }
1860 
1861  /* changing all of the master problem variable to continuous. */
1862  SCIP_CALL( SCIPbendersChgMastervarsToCont(benders, set, i) );
1863 
1864  /* checking the convexity of the subproblem. The convexity of the subproblem indicates whether the convex
1865  * relaxation is a valid relaxation for the problem
1866  */
1867  SCIP_CALL( checkSubproblemConvexity(benders, set, i) );
1868 
1869  /* if the problem is convex and has nonlinear constraints, then slack variables must be added to each of the
1870  * constraints
1871  */
1872  if( benders->execfeasphase ||
1874  && SCIPbendersSubproblemIsNonlinear(benders, i)) )
1875  {
1876  /* the slack variables are only added to the subproblem once. If the initialisation methods are called from a
1877  * copy, then the slack variables are not re-added. Alternatively, if the copy must be threadsafe, then the
1878  * subproblems are created from scratch again, so the slack variables need to be added.
1879  */
1880  if( !benders->iscopy || benders->threadsafe )
1881  {
1882  SCIP_CALL( addSlackVarsToConstraints(benders, set, i) );
1883  }
1884 
1885  /* setting the flag to indicate that slack variables have been added to the subproblem constraints. This is only
1886  * set if the slack variables have been added at the request of the user.
1887  */
1888  if( benders->execfeasphase )
1889  benders->feasibilityphase = TRUE;
1890  }
1891 
1892  /* after checking the subproblem for convexity, if the subproblem has convex constraints and continuous variables,
1893  * then the problem is entered into probing mode. Otherwise, it is initialised as a CIP
1894  */
1896  {
1897  /* if the user has not implemented a solve subproblem callback, then the subproblem solves are performed
1898  * internally. To be more efficient the subproblem is put into probing mode. */
1899  if( benders->benderssolvesubconvex == NULL && benders->benderssolvesub == NULL
1900  && SCIPgetStage(subproblem) <= SCIP_STAGE_PROBLEM )
1901  {
1902  SCIP_CALL( initialiseLPSubproblem(benders, set, i) );
1903  }
1904  }
1905  else
1906  {
1907  SCIP_EVENTHDLRDATA* eventhdlrdata_mipnodefocus;
1908  SCIP_EVENTHDLRDATA* eventhdlrdata_upperbound;
1909 
1910  /* because the subproblems could be reused in the copy, the event handler is not created again. If the
1911  * threadsafe is TRUE, then it is assumed that the subproblems are not reused.
1912  * NOTE: This currently works with the benders_default implementation. It may not be very general. */
1913  if( benders->benderssolvesubconvex == NULL && benders->benderssolvesub == NULL
1914  && (!benders->iscopy || benders->threadsafe) )
1915  {
1916  SCIP_CALL( SCIPallocBlockMemory(subproblem, &eventhdlrdata_mipnodefocus) );
1917  SCIP_CALL( SCIPallocBlockMemory(subproblem, &eventhdlrdata_upperbound) );
1918 
1919  SCIP_CALL( initEventhandlerData(subproblem, eventhdlrdata_mipnodefocus) );
1920  SCIP_CALL( initEventhandlerData(subproblem, eventhdlrdata_upperbound) );
1921 
1922  /* include the first LP solved event handler into the subproblem */
1924  MIPNODEFOCUS_EVENTHDLR_DESC, eventExecBendersMipnodefocus, eventhdlrdata_mipnodefocus) );
1925  SCIP_CALL( SCIPsetEventhdlrInitsol(subproblem, eventhdlr, eventInitsolBendersMipnodefocus) );
1926  SCIP_CALL( SCIPsetEventhdlrExitsol(subproblem, eventhdlr, eventExitsolBendersMipnodefocus) );
1927  SCIP_CALL( SCIPsetEventhdlrExit(subproblem, eventhdlr, eventExitBendersMipnodefocus) );
1928  SCIP_CALL( SCIPsetEventhdlrFree(subproblem, eventhdlr, eventFreeBendersMipnodefocus) );
1929  assert(eventhdlr != NULL);
1930 
1931  /* include the upper bound interrupt event handler into the subproblem */
1933  UPPERBOUND_EVENTHDLR_DESC, eventExecBendersUpperbound, eventhdlrdata_upperbound) );
1934  SCIP_CALL( SCIPsetEventhdlrInitsol(subproblem, eventhdlr, eventInitsolBendersUpperbound) );
1935  SCIP_CALL( SCIPsetEventhdlrExitsol(subproblem, eventhdlr, eventExitsolBendersUpperbound) );
1936  SCIP_CALL( SCIPsetEventhdlrExit(subproblem, eventhdlr, eventExitBendersUpperbound) );
1937  SCIP_CALL( SCIPsetEventhdlrFree(subproblem, eventhdlr, eventFreeBendersUpperbound) );
1938  assert(eventhdlr != NULL);
1939  }
1940  }
1941  }
1942  else
1943  {
1944  /* a user must specify the subproblem type if they are not supplying a SCIP instance. */
1946  {
1947  SCIPerrorMessage("If the subproblem is set to NULL, then the subproblem type must be specified.\n");
1948  SCIPerrorMessage("In the subproblem creation callback, call SCIPbendersSetSubproblemType with the appropriate problem type.\n");
1949 
1950  return SCIP_ERROR;
1951  }
1952  }
1953  }
1954 
1955  /* checking the convexity of the master problem. This information is useful for the cut generation methods, such as
1956  * non-good and integer cuts
1957  */
1958  SCIP_CALL( checkSubproblemConvexity(benders, set, -1) );
1959 
1960  benders->subprobscreated = TRUE;
1961 
1962  return SCIP_OKAY;
1963 }
1964 
1965 
1966 /** initializes Benders' decomposition */
1968  SCIP_BENDERS* benders, /**< Benders' decomposition */
1969  SCIP_SET* set /**< global SCIP settings */
1970  )
1971 {
1972  int i;
1973 
1974  assert(benders != NULL);
1975  assert(set != NULL);
1976 
1977  if( benders->initialized )
1978  {
1979  SCIPerrorMessage("Benders' decomposition <%s> already initialized\n", benders->name);
1980  return SCIP_INVALIDCALL;
1981  }
1982 
1983  if( set->misc_resetstat )
1984  {
1985  SCIPclockReset(benders->setuptime);
1986  SCIPclockReset(benders->bendersclock);
1987 
1988  benders->ncalls = 0;
1989  benders->ncutsfound = 0;
1990  benders->ntransferred = 0;
1991  }
1992 
1993  /* start timing */
1994  SCIPclockStart(benders->setuptime, set);
1995 
1996  if( benders->bendersinit != NULL )
1997  {
1998  SCIP_CALL( benders->bendersinit(set->scip, benders) );
1999  }
2000 
2001  benders->initialized = TRUE;
2002 
2003  /* if the Benders' decomposition is a copy, then the auxiliary variables already exist. So they are registered with
2004  * the Benders' decomposition struct during the init stage. If the Benders' decomposition is not a copy, then the
2005  * auxiliary variables need to be created, which occurs in the initpre stage
2006  */
2007  if( benders->iscopy )
2008  {
2009  /* the copied auxiliary variables must be assigned to the target Benders' decomposition */
2010  SCIP_CALL( assignAuxiliaryVariables(set->scip, benders) );
2011  }
2012 
2013  /* creates the subproblems and sets up the probing mode for LP subproblems. This function calls the benderscreatesub
2014  * callback. */
2015  SCIP_CALL( createSubproblems(benders, set) );
2016 
2017  /* storing the solution tolerance set by the SCIP parameters */
2018  SCIP_CALL( SCIPsetGetRealParam(set, "benders/solutiontol", &benders->solutiontol) );
2019 
2020  /* allocating memory for the stored constraints array */
2021  if( benders->storedcutssize == 0 )
2022  {
2024  benders->storedcutssize = BENDERS_ARRAYSIZE;
2025  benders->nstoredcuts = 0;
2026  }
2027 
2028  /* initialising the Benders' cuts */
2029  SCIPbendersSortBenderscuts(benders);
2030  for( i = 0; i < benders->nbenderscuts; i++ )
2031  {
2032  SCIP_CALL( SCIPbenderscutInit(benders->benderscuts[i], set) );
2033  }
2034 
2035  /* stop timing */
2036  SCIPclockStop(benders->setuptime, set);
2037 
2038  return SCIP_OKAY;
2039 }
2040 
2041 
2042 /** Transfers Benders' cuts that were generated while solving a sub-SCIP to the original SCIP instance. This involves
2043  * creating a constraint/cut that is equivalent to the generated cut in the sub-SCIP. This new constraint/cut is then
2044  * added to the original SCIP instance.
2045  */
2046 static
2048  SCIP* sourcescip, /**< the source SCIP from when the Benders' decomposition was copied */
2049  SCIP_BENDERS* benders, /**< the Benders' decomposition structure of the sub SCIP */
2050  SCIP_VAR** vars, /**< the variables from the source constraint */
2051  SCIP_Real* vals, /**< the coefficients of the variables in the source constriant */
2052  SCIP_Real lhs, /**< the LHS of the source constraint */
2053  SCIP_Real rhs, /**< the RHS of the source constraint */
2054  int nvars /**< the number of variables in the source constraint */
2055  )
2056 {
2057  SCIP_BENDERS* sourcebenders; /* the Benders' decomposition of the source SCIP */
2058  SCIP_CONSHDLR* consbenders; /* a helper variable for the Benders' decomposition constraint handler */
2059  SCIP_CONS* transfercons = NULL; /* the constraint that is generated to transfer the constraints/cuts */
2060  SCIP_ROW* transfercut = NULL; /* the cut that is generated to transfer the constraints/cuts */
2061  SCIP_VAR* sourcevar; /* the source variable that will be added to the transferred cut */
2062  SCIP_VAR* origvar;
2063  SCIP_Real scalar;
2064  SCIP_Real constant;
2065  char cutname[SCIP_MAXSTRLEN]; /* the name of the transferred cut */
2066  int i;
2067  SCIP_Bool fail;
2068 
2069  assert(sourcescip != NULL);
2070  assert(benders != NULL);
2071  assert(vars != NULL);
2072  assert(vals != NULL);
2073 
2074  /* retrieving the source Benders' decomposition structure */
2075  sourcebenders = SCIPfindBenders(sourcescip, SCIPbendersGetName(benders));
2076 
2077  /* retrieving the Benders' decomposition constraint handler */
2078  consbenders = SCIPfindConshdlr(sourcescip, "benders");
2079 
2080  /* setting the name of the transferred cut */
2081  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "transferredcut_%d",
2082  SCIPbendersGetNTransferredCuts(sourcebenders) );
2083 
2084  /* TODO: It could be more efficient to pass an updated vars array with the vals array to the
2085  * SCIPcreateConsBasicLinear/SCIPcreateEmptyRowConshdlr. This should be implemented to improve the performance of the
2086  * Large Neighbourhood Benders Search.
2087  */
2088 
2089  /* creating an empty row/constraint for the transferred cut */
2090  if( sourcebenders->cutsasconss )
2091  {
2092  SCIP_CALL( SCIPcreateConsBasicLinear(sourcescip, &transfercons, cutname, 0, NULL, NULL, lhs, rhs) );
2093  SCIP_CALL( SCIPsetConsRemovable(sourcescip, transfercons, TRUE) );
2094  }
2095  else
2096  {
2097  SCIP_CALL( SCIPcreateEmptyRowConshdlr(sourcescip, &transfercut, consbenders, cutname, lhs, rhs, FALSE,
2098  FALSE, TRUE) );
2099  }
2100 
2101  fail = FALSE;
2102  for( i = 0; i < nvars; i++ )
2103  {
2104  /* getting the original variable for the transformed variable */
2105  origvar = vars[i];
2106  scalar = 1.0;
2107  constant = 0.0;
2108  SCIP_CALL( SCIPvarGetOrigvarSum(&origvar, &scalar, &constant) );
2109 
2110  /* getting the source var from the hash map */
2111  sourcevar = (SCIP_VAR*) SCIPhashmapGetImage(benders->mastervarsmap, origvar);
2112 
2113  /* if the source variable is not found, then the mapping in incomplete. So the constraint can not be
2114  * transferred. */
2115  if( sourcevar == NULL )
2116  {
2117  fail = TRUE;
2118  break;
2119  }
2120 
2121  if( sourcebenders->cutsasconss )
2122  {
2123  assert( transfercons != NULL );
2124  SCIP_CALL( SCIPaddCoefLinear(sourcescip, transfercons, sourcevar, vals[i]) ); /*lint !e644*/
2125  }
2126  else
2127  {
2128  assert( transfercut != NULL );
2129  SCIP_CALL( SCIPaddVarToRow(sourcescip, transfercut, sourcevar, vals[i]) ); /*lint !e644*/
2130  }
2131  }
2132 
2133  /* if all of the source variables were found to generate the cut */
2134  if( !fail )
2135  {
2136  if( sourcebenders->cutsasconss )
2137  {
2138  SCIP_CALL( SCIPaddCons(sourcescip, transfercons) );
2139  }
2140  else
2141  {
2142  SCIP_CALL( SCIPaddPoolCut(sourcescip, transfercut) );
2143  }
2144 
2145  sourcebenders->ntransferred++;
2146  }
2147 
2148  /* release the row/constraint */
2149  if( sourcebenders->cutsasconss )
2150  {
2151  /* only release if the creation of the constraint failed. */
2152  SCIP_CALL( SCIPreleaseCons(sourcescip, &transfercons) );
2153  }
2154  else
2155  {
2156  SCIP_CALL( SCIPreleaseRow(sourcescip, &transfercut) );
2157  }
2158 
2159  return SCIP_OKAY;
2160 }
2161 
2162 
2163 /** transfers the cuts generated in a subscip to the source scip */
2164 static
2166  SCIP* sourcescip, /**< the source SCIP from when the Benders' decomposition was copied */
2167  SCIP* subscip, /**< the sub SCIP where the Benders' cuts were generated */
2168  SCIP_BENDERS* benders /**< the Benders' decomposition structure of the sub SCIP */
2169  )
2170 {
2171  SCIP_BENDERS* sourcebenders; /* the Benders' decomposition of the source SCIP */
2172  SCIP_VAR** vars; /* the variables of the added constraint/row */
2173  SCIP_Real* vals; /* the values of the added constraint/row */
2174  SCIP_Real lhs; /* the LHS of the added constraint/row */
2175  SCIP_Real rhs; /* the RHS of the added constraint/row */
2176  int naddedcuts;
2177  int nvars;
2178  int i;
2179 
2180  assert(subscip != NULL);
2181  assert(benders != NULL);
2182 
2183  /* retrieving the source Benders' decomposition structure */
2184  sourcebenders = SCIPfindBenders(sourcescip, SCIPbendersGetName(benders));
2185 
2186  /* exit if the cuts should not be transferred from the sub SCIP to the source SCIP. */
2187  if( !sourcebenders->transfercuts || benders->mastervarsmap == NULL )
2188  return SCIP_OKAY;
2189 
2190  /* retrieving the number of stored Benders' cuts */
2191  naddedcuts = SCIPbendersGetNStoredCuts(benders);
2192 
2193  /* looping over all added cuts to construct the cut for the source scip */
2194  for( i = 0; i < naddedcuts; i++ )
2195  {
2196  /* collecting the variable information from the constraint */
2197  SCIP_CALL( SCIPbendersGetStoredCutData(benders, i, &vars, &vals, &lhs, &rhs, &nvars) );
2198 
2199  if( nvars > 0 )
2200  {
2201  /* create and add the cut to be transferred from the sub SCIP to the source SCIP */
2202  SCIP_CALL( createAndAddTransferredCut(sourcescip, benders, vars, vals, lhs, rhs, nvars) );
2203  }
2204  }
2205 
2206  return SCIP_OKAY;
2207 }
2208 
2209 
2210 /** calls exit method of Benders' decomposition */
2212  SCIP_BENDERS* benders, /**< Benders' decomposition */
2213  SCIP_SET* set /**< global SCIP settings */
2214  )
2215 {
2216  int nsubproblems;
2217  int i;
2218 
2219  assert(benders != NULL);
2220  assert(set != NULL);
2221 
2222  if( !benders->initialized )
2223  {
2224  SCIPerrorMessage("Benders' decomposition <%s> not initialized\n", benders->name);
2225  return SCIP_INVALIDCALL;
2226  }
2227 
2228  /* start timing */
2229  SCIPclockStart(benders->setuptime, set);
2230 
2231  if( benders->bendersexit != NULL )
2232  {
2233  SCIP_CALL( benders->bendersexit(set->scip, benders) );
2234  }
2235 
2236  /* if the Benders' decomposition is a copy, then is a variable mapping was provided, then the generated cuts will
2237  * be transferred to the source scip
2238  */
2239  if( benders->iscopy && benders->mastervarsmap != NULL )
2240  {
2241  SCIP_CALL( transferBendersCuts(benders->sourcescip, set->scip, benders) );
2242  }
2243 
2244  /* releasing the stored constraints */
2245  for( i = benders->nstoredcuts - 1; i >= 0; i-- )
2246  {
2247  SCIPfreeBlockMemoryArray(set->scip, &benders->storedcuts[i]->vals, benders->storedcuts[i]->nvars);
2248  SCIPfreeBlockMemoryArray(set->scip, &benders->storedcuts[i]->vars, benders->storedcuts[i]->nvars);
2249  SCIPfreeBlockMemory(set->scip, &benders->storedcuts[i]); /*lint !e866*/
2250  }
2251 
2252  BMSfreeBlockMemoryArray(SCIPblkmem(set->scip), &benders->storedcuts, benders->storedcutssize);
2253  benders->storedcutssize = 0;
2254  benders->nstoredcuts = 0;
2255 
2256  /* releasing all of the auxiliary variables */
2257  nsubproblems = SCIPbendersGetNSubproblems(benders);
2258  for( i = 0; i < nsubproblems; i++ )
2259  {
2260  /* it is possible that the master problem is not solved. As such, the auxiliary variables will not be created. So
2261  * we don't need to release the variables
2262  */
2263  if( benders->auxiliaryvars[i] != NULL )
2264  {
2265  /* we need to remove the locks from the auxiliary variables. This will be called always for the highest priority
2266  * Benders' plugin and others if the auxiliary variables are not shared
2267  */
2268  if( !benders->iscopy && SCIPvarGetNLocksDown(benders->auxiliaryvars[i]) > 0 )
2269  SCIP_CALL( SCIPaddVarLocksType(set->scip, benders->auxiliaryvars[i], SCIP_LOCKTYPE_MODEL, -1, 0) );
2270 
2271  SCIP_CALL( SCIPreleaseVar(set->scip, &benders->auxiliaryvars[i]) );
2272  }
2273  }
2274 
2275  /* if a corepoint has been used for cut strengthening, then this needs to be freed */
2276  if( benders->corepoint != NULL )
2277  {
2278  SCIP_CALL( SCIPfreeSol(set->scip, &benders->corepoint) );
2279  }
2280 
2281  /* calling the exit method for the Benders' cuts */
2282  SCIPbendersSortBenderscuts(benders);
2283  for( i = 0; i < benders->nbenderscuts; i++ )
2284  {
2285  SCIP_CALL( SCIPbenderscutExit(benders->benderscuts[i], set) );
2286  }
2287 
2288  benders->initialized = FALSE;
2289 
2290  /* stop timing */
2291  SCIPclockStop(benders->setuptime, set);
2292 
2293  return SCIP_OKAY;
2294 }
2295 
2296 /** Checks whether a subproblem is independent. */
2297 static
2299  SCIP* scip, /**< the SCIP data structure */
2300  SCIP_BENDERS* benders /**< Benders' decomposition */
2301  )
2302 {
2303  SCIP_VAR** vars;
2304  int nvars;
2305  int nsubproblems;
2306  int i;
2307  int j;
2308 
2309  assert(scip != NULL);
2310  assert(benders != NULL);
2311 
2312  /* retrieving the master problem variables */
2313  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2314 
2315  nsubproblems = SCIPbendersGetNSubproblems(benders);
2316 
2317  /* looping over all subproblems to check whether there exists at least one master problem variable */
2318  for( i = 0; i < nsubproblems; i++ )
2319  {
2320  SCIP_Bool independent = FALSE;
2321 
2322  /* if there are user defined solving or freeing functions, then it is not possible to declare the independence of
2323  * the subproblems.
2324  */
2325  if( benders->benderssolvesubconvex == NULL && benders->benderssolvesub == NULL
2326  && benders->bendersfreesub == NULL )
2327  {
2328  independent = TRUE;
2329 
2330  for( j = 0; j < nvars; j++ )
2331  {
2332  SCIP_VAR* subprobvar;
2333 
2334  /* getting the subproblem problem variable corresponding to the master problem variable */
2335  SCIP_CALL( SCIPgetBendersSubproblemVar(scip, benders, vars[j], &subprobvar, i) );
2336 
2337  /* if the subporblem variable is not NULL, then the subproblem depends on the master problem */
2338  if( subprobvar != NULL )
2339  {
2340  independent = FALSE;
2341  break;
2342  }
2343  }
2344 
2345  /* setting the independent flag */
2346  SCIPbendersSetSubproblemIsIndependent(benders, i, independent);
2347  }
2348  }
2349 
2350  return SCIP_OKAY;
2351 }
2352 
2353 /** informs the Benders' decomposition that the presolving process is being started */
2355  SCIP_BENDERS* benders, /**< Benders' decomposition */
2356  SCIP_SET* set, /**< global SCIP settings */
2357  SCIP_STAT* stat /**< dynamic problem statistics */
2358  )
2359 {
2360  assert(benders != NULL);
2361  assert(set != NULL);
2362  assert(stat != NULL);
2363 
2364  /* if the Benders' decomposition is the original, then the auxiliary variables need to be created. If the Benders'
2365  * decomposition is a copy, then the auxiliary variables already exist. The assignment of the auxiliary variables
2366  * occurs in bendersInit
2367  */
2368  if( !benders->iscopy )
2369  {
2370  /* check the subproblem independence. This check is only performed if the user has not implemented a solve
2371  * subproblem function.
2372  */
2373  if( benders->benderssolvesubconvex == NULL && benders->benderssolvesub == NULL )
2374  SCIP_CALL( checkSubproblemIndependence(set->scip, benders) );
2375 
2376  /* adding the auxiliary variables to the master problem */
2377  SCIP_CALL( addAuxiliaryVariablesToMaster(set->scip, benders) );
2378  }
2379 
2380  /* call presolving initialization method of Benders' decomposition */
2381  if( benders->bendersinitpre != NULL )
2382  {
2383  /* start timing */
2384  SCIPclockStart(benders->setuptime, set);
2385 
2386  SCIP_CALL( benders->bendersinitpre(set->scip, benders) );
2387 
2388  /* stop timing */
2389  SCIPclockStop(benders->setuptime, set);
2390  }
2391 
2392  return SCIP_OKAY;
2393 }
2394 
2395 
2396 /** informs the Benders' decomposition that the presolving process has completed */
2398  SCIP_BENDERS* benders, /**< Benders' decomposition */
2399  SCIP_SET* set, /**< global SCIP settings */
2400  SCIP_STAT* stat /**< dynamic problem statistics */
2401  )
2402 {
2403  assert(benders != NULL);
2404  assert(set != NULL);
2405  assert(stat != NULL);
2406 
2407  /* call presolving deinitialization method of Benders' decomposition */
2408  if( benders->bendersexitpre != NULL )
2409  {
2410  /* start timing */
2411  SCIPclockStart(benders->setuptime, set);
2412 
2413  SCIP_CALL( benders->bendersexitpre(set->scip, benders) );
2414 
2415  /* stop timing */
2416  SCIPclockStop(benders->setuptime, set);
2417  }
2418 
2419  return SCIP_OKAY;
2420 }
2421 
2422 /** informs Benders' decomposition that the branch and bound process is being started */
2424  SCIP_BENDERS* benders, /**< Benders' decomposition */
2425  SCIP_SET* set /**< global SCIP settings */
2426  )
2427 {
2428  int i;
2429 
2430  assert(benders != NULL);
2431  assert(set != NULL);
2432 
2433  /* call solving process initialization method of Benders' decomposition */
2434  if( benders->bendersinitsol != NULL )
2435  {
2436  /* start timing */
2437  SCIPclockStart(benders->setuptime, set);
2438 
2439  SCIP_CALL( benders->bendersinitsol(set->scip, benders) );
2440 
2441  /* stop timing */
2442  SCIPclockStop(benders->setuptime, set);
2443  }
2444 
2445  /* calling the initsol method for the Benders' cuts */
2446  SCIPbendersSortBenderscuts(benders);
2447  for( i = 0; i < benders->nbenderscuts; i++ )
2448  {
2449  SCIP_CALL( SCIPbenderscutInitsol(benders->benderscuts[i], set) );
2450  }
2451 
2452  return SCIP_OKAY;
2453 }
2454 
2455 /** informs Benders' decomposition that the branch and bound process data is being freed */
2457  SCIP_BENDERS* benders, /**< Benders' decomposition */
2458  SCIP_SET* set /**< global SCIP settings */
2459  )
2460 {
2461  int nsubproblems;
2462  int i;
2463 
2464  assert(benders != NULL);
2465  assert(set != NULL);
2466 
2467  nsubproblems = SCIPbendersGetNSubproblems(benders);
2468  /* freeing all subproblems that are independent, this is because they have not bee freed during the subproblem
2469  * solving loop.
2470  */
2471  for( i = 0; i < nsubproblems; i++ )
2472  {
2473  if( SCIPbendersSubproblemIsIndependent(benders, i) )
2474  {
2475  /* disabling the independence of the subproblem so that it can be freed */
2477 
2478  /* freeing the independent subproblem */
2479  SCIP_CALL( SCIPbendersFreeSubproblem(benders, set, i) );
2480  }
2481  }
2482 
2483  /* call solving process deinitialization method of Benders' decomposition */
2484  if( benders->bendersexitsol != NULL )
2485  {
2486  /* start timing */
2487  SCIPclockStart(benders->setuptime, set);
2488 
2489  SCIP_CALL( benders->bendersexitsol(set->scip, benders) );
2490 
2491  /* stop timing */
2492  SCIPclockStop(benders->setuptime, set);
2493  }
2494 
2495  /* sorting the Benders' decomposition cuts in order of priority. Only a single cut is generated for each subproblem
2496  * per solving iteration. This is particularly important in the case of the optimality and feasibility cuts. Since
2497  * these work on two different solutions to the subproblem, it is not necessary to generate both cuts. So, once the
2498  * feasibility cut is generated, then no other cuts will be generated.
2499  */
2500  SCIPbendersSortBenderscuts(benders);
2501 
2502  /* calling the exitsol method for the Benders' cuts */
2503  for( i = 0; i < benders->nbenderscuts; i++ )
2504  {
2505  SCIP_CALL( SCIPbenderscutExitsol(benders->benderscuts[i], set) );
2506  }
2507 
2508  return SCIP_OKAY;
2509 }
2510 
2511 /** activates Benders' decomposition such that it is called in LP solving loop */
2513  SCIP_BENDERS* benders, /**< the Benders' decomposition structure */
2514  SCIP_SET* set, /**< global SCIP settings */
2515  int nsubproblems /**< the number subproblems used in this decomposition */
2516  )
2517 {
2518  SCIP_EVENTHDLR* eventhdlr;
2519  SCIP_EVENTHDLRDATA* eventhdlrdata;
2520  int i;
2521 
2522  assert(benders != NULL);
2523  assert(set != NULL);
2524  assert(set->stage == SCIP_STAGE_INIT || set->stage == SCIP_STAGE_PROBLEM);
2525 
2526  if( !benders->active )
2527  {
2528  benders->active = TRUE;
2529  set->nactivebenders++;
2530  set->benderssorted = FALSE;
2531 
2532  benders->nsubproblems = nsubproblems;
2533  benders->nactivesubprobs = nsubproblems;
2534  benders->prevlowerbound = -SCIPsetInfinity(set);
2535  benders->strengthenround = FALSE;
2536 
2537  /* allocating memory for the subproblems arrays */
2538  SCIP_ALLOC( BMSallocMemoryArray(&benders->subproblems, benders->nsubproblems) );
2539  SCIP_ALLOC( BMSallocMemoryArray(&benders->auxiliaryvars, benders->nsubproblems) );
2540  SCIP_ALLOC( BMSallocMemoryArray(&benders->solvestat, benders->nsubproblems) );
2541  SCIP_ALLOC( BMSallocMemoryArray(&benders->subprobobjval, benders->nsubproblems) );
2544  SCIP_ALLOC( BMSallocMemoryArray(&benders->subprobtype, benders->nsubproblems) );
2547  SCIP_ALLOC( BMSallocMemoryArray(&benders->subprobsetup, benders->nsubproblems) );
2548  SCIP_ALLOC( BMSallocMemoryArray(&benders->indepsubprob, benders->nsubproblems) );
2551 
2552  /* creating the priority queue for the subproblem solving status */
2553  SCIP_CALL( SCIPpqueueCreate(&benders->subprobqueue, benders->nsubproblems, 1.1,
2554  benders->benderssubcomp == NULL ? benderssubcompdefault : benders->benderssubcomp, NULL) );
2555 
2556  for( i = 0; i < benders->nsubproblems; i++ )
2557  {
2558  SCIP_SUBPROBLEMSOLVESTAT* solvestat;
2559 
2560  benders->subproblems[i] = NULL;
2561  benders->auxiliaryvars[i] = NULL;
2562  benders->subprobobjval[i] = SCIPsetInfinity(set);
2563  benders->bestsubprobobjval[i] = SCIPsetInfinity(set);
2564  benders->subproblowerbound[i] = -SCIPsetInfinity(set);
2566  benders->subprobisconvex[i] = FALSE;
2567  benders->subprobisnonlinear[i] = FALSE;
2568  benders->subprobsetup[i] = FALSE;
2569  benders->indepsubprob[i] = FALSE;
2570  benders->subprobenabled[i] = TRUE;
2571  benders->mastervarscont[i] = FALSE;
2572 
2573  /* initialising the subproblem solving status */
2574  SCIP_ALLOC( BMSallocMemory(&solvestat) );
2575  solvestat->idx = i;
2576  solvestat->ncalls = 0;
2577  solvestat->avgiter = 0;
2578  benders->solvestat[i] = solvestat;
2579 
2580  /* inserting the initial elements into the priority queue */
2581  SCIP_CALL( SCIPpqueueInsert(benders->subprobqueue, benders->solvestat[i]) );
2582  }
2583 
2584  /* adding an eventhandler for updating the lower bound when the root node is solved. */
2585  eventhdlrdata = (SCIP_EVENTHDLRDATA*)benders;
2586 
2587  /* include event handler into SCIP */
2589  eventExecBendersNodesolved, eventhdlrdata) );
2590  SCIP_CALL( SCIPsetEventhdlrInitsol(set->scip, eventhdlr, eventInitsolBendersNodesolved) );
2591  assert(eventhdlr != NULL);
2592  }
2593 
2594  return SCIP_OKAY;
2595 }
2596 
2597 /** deactivates Benders' decomposition such that it is no longer called in LP solving loop */
2599  SCIP_BENDERS* benders, /**< the Benders' decomposition structure */
2600  SCIP_SET* set /**< global SCIP settings */
2601  )
2602 {
2603  int i;
2604 
2605  assert(benders != NULL);
2606  assert(set != NULL);
2607  assert(set->stage == SCIP_STAGE_INIT || set->stage == SCIP_STAGE_PROBLEM);
2608 
2609  if( benders->active )
2610  {
2611  int nsubproblems;
2612 
2613  nsubproblems = SCIPbendersGetNSubproblems(benders);
2614 
2615 #ifndef NDEBUG
2616  /* checking whether the auxiliary variables and subproblems are all NULL */
2617  for( i = 0; i < nsubproblems; i++ )
2618  assert(benders->auxiliaryvars[i] == NULL);
2619 #endif
2620 
2621  /* if the subproblems were created by the Benders' decomposition core, then they need to be freed */
2622  if( benders->freesubprobs )
2623  {
2624  for( i = SCIPbendersGetNSubproblems(benders) - 1; i >= 0; i-- )
2625  {
2626  SCIP* subproblem = SCIPbendersSubproblem(benders, i);
2627  SCIP_CALL( SCIPfree(&subproblem) );
2628  }
2629  }
2630 
2631  benders->active = FALSE;
2632  set->nactivebenders--;
2633  set->benderssorted = FALSE;
2634 
2635  /* freeing the priority queue memory */
2636  SCIPpqueueFree(&benders->subprobqueue);
2637 
2638  for( i = nsubproblems - 1; i >= 0; i-- )
2639  BMSfreeMemory(&benders->solvestat[i]);
2640 
2641  /* freeing the memory allocated during the activation of the Benders' decomposition */
2644  BMSfreeMemoryArray(&benders->indepsubprob);
2645  BMSfreeMemoryArray(&benders->subprobsetup);
2648  BMSfreeMemoryArray(&benders->subprobtype);
2651  BMSfreeMemoryArray(&benders->subprobobjval);
2652  BMSfreeMemoryArray(&benders->auxiliaryvars);
2653  BMSfreeMemoryArray(&benders->solvestat);
2654  BMSfreeMemoryArray(&benders->subproblems);
2655  }
2656 
2657  return SCIP_OKAY;
2658 }
2659 
2660 /** returns whether the given Benders' decomposition is in use in the current problem */
2662  SCIP_BENDERS* benders /**< the Benders' decomposition structure */
2663  )
2664 {
2665  assert(benders != NULL);
2666 
2667  return benders->active;
2668 }
2669 
2670 /** updates the lower bound for all auxiliary variables. This is called if the first LP enforced is unbounded. */
2671 static
2673  SCIP_BENDERS* benders, /**< Benders' decomposition */
2674  SCIP_SET* set, /**< global SCIP settings */
2675  SCIP_RESULT* result /**< the result from updating the auxiliary variable lower bound */
2676  )
2677 {
2678  int nsubproblems;
2679  int i;
2680 
2681  assert(benders != NULL);
2682  assert(set != NULL);
2683 
2684  (*result) = SCIP_DIDNOTRUN;
2685 
2686  nsubproblems = SCIPbendersGetNSubproblems(benders);
2687 
2688  for( i = 0; i < nsubproblems; i++ )
2689  {
2690  SCIP_VAR* auxiliaryvar;
2691  SCIP_Real lowerbound;
2692  SCIP_Bool infeasible;
2693 
2694  infeasible = FALSE;
2695 
2696  /* computing the lower bound of the subproblem by solving it without any variable fixings */
2697  SCIP_CALL( SCIPbendersComputeSubproblemLowerbound(benders, set, i, &lowerbound, &infeasible) );
2698 
2699  /* if the subproblem is infeasible, then the original problem is infeasible */
2700  if( infeasible )
2701  {
2702  (*result) = SCIP_INFEASIBLE;
2703  break;
2704  }
2705 
2706  /* retrieving the auxiliary variable */
2707  auxiliaryvar = SCIPbendersGetAuxiliaryVar(benders, i);
2708 
2709  /* only update the lower bound if it is greater than the current lower bound */
2710  if( SCIPsetIsGT(set, lowerbound, SCIPvarGetLbGlobal(auxiliaryvar)) )
2711  {
2712  SCIPsetDebugMsg(set, "Tightened lower bound of <%s> to %g\n", SCIPvarGetName(auxiliaryvar), lowerbound);
2713  /* updating the lower bound of the auxiliary variable */
2714  SCIP_CALL( SCIPchgVarLb(set->scip, auxiliaryvar, lowerbound) );
2715  (*result) = SCIP_REDUCEDDOM;
2716  }
2717 
2718  /* stores the lower bound for the subproblem */
2719  SCIPbendersUpdateSubproblemLowerbound(benders, i, lowerbound);
2720  }
2721 
2722  return SCIP_OKAY;
2723 }
2724 
2725 /** sets the core point used for cut strengthening. If the strenghtenintpoint is set to 'i', then the core point is
2726  * reinitialised each time the incumbent is updated
2727  */
2728 static
2730  SCIP* scip, /**< the SCIP data structure */
2731  SCIP_BENDERS* benders /**< Benders' decomposition */
2732  )
2733 {
2734  SCIP_SOL* bestsol;
2735 
2736  assert(scip != NULL);
2737  assert(benders != NULL);
2738 
2739  /* if the core point is not NULL and the interior point is not reinitialised, then nothing is done */
2740  if( benders->corepoint != NULL && benders->strengthenintpoint != 'i' )
2741  return SCIP_OKAY;
2742 
2743  bestsol = SCIPgetBestSol(scip);
2744 
2745  /* if the core point should be updated, then this only happens if the incumbent solution has been updated */
2746  if( benders->strengthenintpoint == 'i' && benders->initcorepoint == bestsol )
2747  return SCIP_OKAY;
2748 
2749  /* if a corepoint has been used for cut strengthening, then this needs to be freed */
2750  if( benders->corepoint != NULL )
2751  {
2752  SCIP_CALL( SCIPfreeSol(scip, &benders->corepoint) );
2753  }
2754 
2755  switch( benders->strengthenintpoint )
2756  {
2757  SCIP_VAR** vars;
2758  SCIP_Real timelimit;
2759  int nvars;
2760  int i;
2761 
2762  case 'l':
2763  SCIP_CALL( SCIPcreateLPSol(scip, &benders->corepoint, NULL) );
2764  SCIP_CALL( SCIPunlinkSol(scip, benders->corepoint) );
2765  break;
2766  case 'f':
2767  case 'i':
2768  SCIP_CALL( SCIPcreateSolCopy(scip, &benders->corepoint, bestsol) );
2769  SCIP_CALL( SCIPunlinkSol(scip, benders->corepoint) );
2770  benders->initcorepoint = bestsol;
2771  break;
2772  case 'r':
2773  /* prepare time limit */
2774  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
2775  if ( ! SCIPisInfinity(scip, timelimit) )
2776  timelimit -= SCIPgetSolvingTime(scip);
2777 
2778  /* if there is time remaining, then compute the relative interior point. Otherwise, return the LP solution */
2779  if ( timelimit > 0.0 )
2780  {
2781  SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, 0, "Computing relative interior point (time limit: %g, iter limit: %d) ...\n", timelimit, INT_MAX);
2782  SCIP_CALL( SCIPcomputeLPRelIntPoint(scip, TRUE, FALSE, timelimit, INT_MAX, &benders->corepoint) );
2783  }
2784  else
2785  {
2786  SCIP_CALL( SCIPcreateLPSol(scip, &benders->corepoint, NULL) );
2787  SCIP_CALL( SCIPunlinkSol(scip, benders->corepoint) );
2788  }
2789  break;
2790  case 'z':
2791  SCIP_CALL( SCIPcreateSol(scip, &benders->corepoint, NULL) );
2792  break;
2793  case 'o':
2794  SCIP_CALL( SCIPcreateSol(scip, &benders->corepoint, NULL) );
2795 
2796  /* getting the variable data so that the */
2797  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2798 
2799  /* setting all variable values to 1.0 */
2800  for( i = 0; i < nvars; i++ )
2801  {
2802  SCIP_CALL( SCIPsetSolVal(scip, benders->corepoint, vars[i], 1.0) );
2803  }
2804  break;
2805  default:
2806  SCIP_CALL( SCIPcreateLPSol(scip, &benders->corepoint, NULL) );
2807  SCIP_CALL( SCIPunlinkSol(scip, benders->corepoint) );
2808  }
2809 
2810  return SCIP_OKAY;
2811 }
2812 
2813 /** performs cut strengthening by using an interior solution to generate cuts */
2814 static
2816  SCIP_BENDERS* benders, /**< Benders' decomposition */
2817  SCIP_SET* set, /**< global SCIP settings */
2818  SCIP_SOL* sol, /**< primal CIP solution */
2819  SCIP_BENDERSENFOTYPE type, /**< the type of solution being enforced */
2820  SCIP_Bool checkint, /**< are the subproblems called during a check/enforce of integer sols? */
2821  SCIP_Bool perturbsol, /**< should the solution be perturbed to escape infeasibility? */
2822  SCIP_Bool* auxviol, /**< set to TRUE only if the solution is feasible but the aux vars are violated */
2823  SCIP_Bool* infeasible, /**< is the master problem infeasible with respect to the Benders' cuts? */
2824  SCIP_Bool* skipsolve, /**< should the main solve be skipped as a result of this strengthening? */
2825  SCIP_RESULT* result /**< result of the pricing process */
2826  )
2827 {
2828  SCIP_SOL* sepapoint;
2829  SCIP_VAR** vars;
2830  int prevcutsfound;
2831  int nvars;
2832  int i;
2833 
2834  assert(benders != NULL);
2835  assert(set != NULL);
2836 
2837  (*result) = SCIP_DIDNOTRUN;
2838  (*skipsolve) = FALSE;
2839 
2840  /* the cut stabilisation is only performed when enforcing LP solutions. The solution is not NULL if the stabilisation
2841  * is currently being performed. It is important to avoid recursion
2842  */
2843  if( type != SCIP_BENDERSENFOTYPE_LP || sol != NULL )
2844  return SCIP_OKAY;
2845 
2846  /* checking if a change to the lower bound has occurred */
2847  if( SCIPsetIsGT(set, SCIPgetLowerbound(set->scip), benders->prevlowerbound)
2848  || SCIPgetCurrentNode(set->scip) != benders->prevnode )
2849  {
2850  benders->prevnode = SCIPgetCurrentNode(set->scip);
2851  benders->prevlowerbound = SCIPgetLowerbound(set->scip);
2852  benders->noimprovecount = 0;
2853  }
2854  else
2855  benders->noimprovecount++;
2856 
2857  /* if the number of iterations without improvement exceeds 3*noimprovelimit, then the no stabilisation is performed
2858  */
2859  if( benders->noimprovecount > 3*benders->noimprovelimit )
2860  return SCIP_OKAY;
2861 
2862  /* if there is no incumbent solution, then it is not possible to create the core point and hence the strengthening
2863  * can not be performed
2864  */
2865  if( SCIPgetBestSol(set->scip) == NULL )
2866  return SCIP_OKAY;
2867 
2868  /* if no LP iterations have been performed since the last call of the cut strenghtening, then the strengthening is
2869  * aborted
2870  */
2871  if( benders->prevnlpiter == SCIPgetNLPIterations(set->scip) )
2872  return SCIP_OKAY;
2873 
2874  benders->prevnlpiter = SCIPgetNLPIterations(set->scip);
2875 
2876  /* if the separation point solution is NULL, then we create the solution using the current LP relaxation. */
2877  SCIP_CALL( setAndUpdateCorePoint(set->scip, benders) );
2878 
2879  /* creating the separation point
2880  * TODO: This could be a little to memory heavy, it may be better just to create the separation point once and then
2881  * update it each time.
2882  */
2883  SCIP_CALL( SCIPcreateLPSol(set->scip, &sepapoint, NULL) );
2884  SCIP_CALL( SCIPunlinkSol(set->scip, sepapoint) );
2885 
2886  SCIP_CALL( SCIPgetVarsData(set->scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2887  assert(vars != NULL);
2888 
2889  /* creating a solution that is a convex combination of the LP solution and the separation point */
2890  for( i = 0; i < nvars; i++ )
2891  {
2892  SCIP_VAR* subvar;
2893  SCIP_Real corepointval;
2894  SCIP_Real lpsolval;
2895  SCIP_Real newsolval;
2896  int j;
2897 
2898  corepointval = SCIPgetSolVal(set->scip, benders->corepoint, vars[i]);
2899  lpsolval = SCIPgetSolVal(set->scip, sol, vars[i]);
2900  newsolval = lpsolval;
2901 
2902  /* checking whether the master variable is mapped to any subproblem variables */
2903  subvar = NULL;
2904  j = 0;
2905  while( subvar == NULL && j < SCIPgetBendersNSubproblems(set->scip, benders) )
2906  {
2907  SCIP_CALL( SCIPgetBendersSubproblemVar(set->scip, benders, vars[i], &subvar, j) );
2908  j++;
2909  }
2910 
2911  /* if the variable is a linking variable and it is not fixed, then a convex combination with the corepoint is
2912  * computed.
2913  */
2914  if( subvar != NULL && SCIPvarGetStatus(vars[i]) != SCIP_VARSTATUS_FIXED )
2915  {
2916  /* if the number of iterations without improvement exceeds noimprovelimit, then no convex combination is
2917  * created
2918  */
2919  if( !perturbsol && benders->noimprovecount <= benders->noimprovelimit )
2920  {
2921  newsolval = lpsolval*benders->convexmult + corepointval*(1 - benders->convexmult);
2922 
2923  /* updating the core point */
2924  SCIP_CALL( SCIPsetSolVal(set->scip, benders->corepoint, vars[i], newsolval) );
2925  }
2926 
2927  /* if the number of iterations without improvement is less than 2*noimprovelimit, then perturbation is
2928  * performed
2929  * TODO: This should be a random vector!!!!
2930  */
2931  if( perturbsol || benders->noimprovecount <= 2*benders->noimprovelimit )
2932  newsolval += benders->perturbeps;
2933  }
2934 
2935  /* updating the separation point */
2936  SCIP_CALL( SCIPsetSolVal(set->scip, sepapoint, vars[i], newsolval) );
2937  }
2938 
2939  /* storing the number of cuts found */
2940  prevcutsfound = SCIPbendersGetNCutsFound(benders);
2941 
2942  SCIPsetDebugMsg(set, "solving Benders' decomposition subproblems with stabilised point.\n");
2943 
2944  /* calling the subproblem solving method to generate cuts from the separation solution */
2945  SCIP_CALL( SCIPsolveBendersSubproblems(set->scip, benders, sepapoint, result, infeasible, auxviol, type, checkint) );
2946 
2947  SCIPsetDebugMsg(set, "solved Benders' decomposition subproblems with stabilised point. noimprovecount %d result %d\n",
2948  benders->noimprovecount, (*result));
2949 
2950  /* if constraints were added, then the main Benders' solving loop is skipped. */
2951  if( !(*infeasible) && ((*result) == SCIP_CONSADDED || (*result) == SCIP_SEPARATED) )
2952  (*skipsolve) = TRUE;
2953 
2954  /* capturing cut strengthening statistics */
2955  benders->nstrengthencalls++;
2956  benders->nstrengthencuts += (SCIPbendersGetNCutsFound(benders) - prevcutsfound);
2957 
2958  /* if no cuts were added, then the strengthening round is marked as failed */
2959  if( SCIPbendersGetNCutsFound(benders) == prevcutsfound )
2960  benders->nstrengthenfails++;
2961 
2962  /* freeing the sepapoint solution */
2963  SCIP_CALL( SCIPfreeSol(set->scip, &sepapoint) );
2964 
2965  return SCIP_OKAY;
2966 }
2967 
2968 
2969 /** Returns whether only the convex relaxations will be checked in this solve loop
2970  * when Benders' is used in the LNS heuristics, only the convex relaxations of the master/subproblems are checked,
2971  * i.e. no integer cuts are generated. In this case, then Benders' decomposition is performed under the assumption
2972  * that all subproblems are convex relaxations.
2973  */
2975  SCIP_BENDERS* benders, /**< Benders' decomposition */
2976  SCIP_Bool subscipsoff /**< flag indicating whether plugins using sub-SCIPs are deactivated */
2977  )
2978 {
2979  return benders->iscopy && benders->lnscheck && subscipsoff;
2980 }
2981 
2982 /** returns the number of subproblems that will be checked in this iteration */
2983 static
2985  SCIP_BENDERS* benders, /**< Benders' decomposition */
2986  SCIP_SET* set, /**< global SCIP settings */
2987  SCIP_BENDERSENFOTYPE type /**< the type of solution being enforced */
2988  )
2989 {
2990  if( benders->ncalls == 0 || type == SCIP_BENDERSENFOTYPE_CHECK
2992  return SCIPbendersGetNSubproblems(benders);
2993  else
2994  return (int) SCIPsetCeil(set, (SCIP_Real) SCIPbendersGetNSubproblems(benders)*benders->subprobfrac);
2995 }
2996 
2997 /** returns whether the solving of the given subproblem needs to be executed */
2998 static
3000  SCIP_BENDERS* benders, /**< Benders' decomposition */
3001  int probnumber /**< the subproblem index */
3002  )
3003 {
3004  return (!SCIPbendersSubproblemIsIndependent(benders, probnumber)
3005  && SCIPbendersSubproblemIsEnabled(benders, probnumber));
3006 }
3007 
3008 /** creates an ordered list of subproblem indices to be solved */
3009 static
3011  SCIP_BENDERS* benders, /**< Benders' decomposition */
3012  SCIP_SET* set, /**< global SCIP settings */
3013  SCIP_BENDERSENFOTYPE type, /**< the type of solution being enforced */
3014  int** solveidx, /**< a list of subproblem indices to the solved in the current iteration */
3015  int* nsolveidx /**< the number of subproblem indices in the list */
3016  )
3017 {
3018  int nsubproblems;
3019  int numtocheck;
3020  int subproblemcount;
3021 
3022  assert(benders != NULL);
3023  assert(set != NULL);
3024  assert((*solveidx) != NULL);
3025  assert(nsolveidx != NULL);
3026  assert(SCIPpqueueNElems(benders->subprobqueue) <= SCIPbendersGetNSubproblems(benders));
3027 
3028  nsubproblems = SCIPbendersGetNSubproblems(benders);
3029 
3030  /* it is possible to only solve a subset of subproblems. This is given by a parameter. */
3031  numtocheck = numSubproblemsToCheck(benders, set, type);
3032 
3033  (*nsolveidx) = 0;
3034 
3035  subproblemcount = 0;
3036  while( subproblemcount < nsubproblems && subproblemcount < numtocheck )
3037  {
3038  SCIP_SUBPROBLEMSOLVESTAT* solvestat;
3039 
3040  solvestat = (SCIP_SUBPROBLEMSOLVESTAT*)SCIPpqueueRemove(benders->subprobqueue);
3041  (*solveidx)[(*nsolveidx)] = solvestat->idx;
3042  (*nsolveidx)++;
3043 
3044  subproblemcount++;
3045  }
3046 }
3047 
3048 /** updates the subproblem solving statistics and inserts the indices into the queue */
3049 static
3051  SCIP_BENDERS* benders, /**< Benders' decomposition */
3052  int* solveidx, /**< the list of indices of subproblems that were solved */
3053  int nsolveidx, /**< the number of subproblem indices */
3054  SCIP_Bool updatestat /**< should the statistics be updated */
3055  )
3056 {
3057  int i;
3058 
3059  assert(benders != NULL);
3060  assert(solveidx != NULL);
3061 
3062  for( i = 0; i < nsolveidx; i++ )
3063  {
3064  SCIP* subproblem;
3065  SCIP_SUBPROBLEMSOLVESTAT* solvestat;
3066 
3067  subproblem = SCIPbendersSubproblem(benders, solveidx[i]);
3068  solvestat = benders->solvestat[solveidx[i]];
3069  assert(solvestat->idx == solveidx[i]);
3070 
3071  /* updating the solving statistics */
3072  if( updatestat )
3073  {
3074  if( subproblem == NULL )
3075  solvestat->avgiter = 1;
3076  else
3077  solvestat->avgiter = (SCIP_Real)(solvestat->avgiter*solvestat->ncalls + SCIPgetNLPIterations(subproblem))
3078  /(SCIP_Real)(solvestat->ncalls + 1);
3079  solvestat->ncalls++;
3080  }
3081 
3082  /* inserting the solving statistics into the priority queue */
3083  SCIP_CALL( SCIPpqueueInsert(benders->subprobqueue, solvestat) );
3084  }
3085 
3086  assert(SCIPpqueueNElems(benders->subprobqueue) == SCIPbendersGetNSubproblems(benders));
3087 
3088  return SCIP_OKAY;
3089 }
3090 
3091 /** Solves each of the Benders' decomposition subproblems for the given solution. All, or a fraction, of subproblems are
3092  * solved before the Benders' decomposition cuts are generated.
3093  * Since a convex relaxation of the subproblem could be solved to generate cuts, a parameter nverified is used to
3094  * identified the number of subproblems that have been solved in their "original" form. For example, if the subproblem
3095  * is a MIP, then if the LP is solved to generate cuts, this does not constitute a verification. The verification is
3096  * only performed when the MIP is solved.
3097  */
3098 static
3100  SCIP_BENDERS* benders, /**< Benders' decomposition */
3101  SCIP_SET* set, /**< global SCIP settings */
3102  SCIP_SOL* sol, /**< primal CIP solution */
3103  SCIP_BENDERSENFOTYPE type, /**< the type of solution being enforced */
3104  SCIP_BENDERSSOLVELOOP solveloop, /**< the current solve loop */
3105  SCIP_Bool checkint, /**< are the subproblems called during a check/enforce of integer sols? */
3106  int* nverified, /**< the number of subproblems verified in the current loop */
3107  int* solveidx, /**< the indices of subproblems to be solved in this loop */
3108  int nsolveidx, /**< the number of subproblems to be solved in this loop */
3109  SCIP_Bool** subprobsolved, /**< an array indicating the subproblems that were solved in this loop. */
3110  SCIP_BENDERSSUBSTATUS** substatus, /**< array to store the status of the subsystem */
3111  SCIP_Bool* infeasible, /**< is the master problem infeasible with respect to the Benders' cuts? */
3112  SCIP_Bool* optimal, /**< is the current solution optimal? */
3113  SCIP_Bool* stopped /**< was the solving process stopped? */
3114  )
3115 {
3116  SCIP_Bool onlyconvexcheck;
3117 #ifdef _OPENMP
3118  int numthreads;
3119  int maxnthreads;
3120 #endif
3121  int i;
3122  int j;
3123 
3124  /* local variables for parallelisation of the solving loop */
3125  int locnverified = *nverified;
3126  SCIP_Bool locinfeasible = *infeasible;
3127  SCIP_Bool locoptimal = *optimal;
3128  SCIP_Bool locstopped = *stopped;
3129 
3130  SCIP_RETCODE retcode = SCIP_OKAY;
3131 
3132  assert(benders != NULL);
3133  assert(set != NULL);
3134 
3135  /* getting the number of threads to use when solving the subproblems. This will be either be
3136  * min(numthreads, maxnthreads).
3137  * NOTE: This may not be correct. The Benders' decomposition parallelisation should not take all minimum threads if
3138  * they are specified. The number of threads should be specified with the Benders' decomposition parameters.
3139  */
3140 #ifdef _OPENMP
3141  SCIP_CALL( SCIPsetGetIntParam(set, "parallel/maxnthreads", &maxnthreads) );
3142  numthreads = MIN(benders->numthreads, maxnthreads);
3143 #endif
3144 
3145  /* in the case of an LNS check, only the convex relaxations of the subproblems will be solved. This is a performance
3146  * feature, since solving the convex relaxation is typically much faster than solving the corresponding CIP. While
3147  * the CIP is not solved during the LNS check, the solutions are still of higher quality than when Benders' is not
3148  * employed.
3149  */
3150  onlyconvexcheck = SCIPbendersOnlyCheckConvexRelax(benders, SCIPsetGetSubscipsOff(set));
3151 
3152  SCIPsetDebugMsg(set, "Performing the subproblem solving process. Number of subproblems to check %d\n", nsolveidx);
3153 
3154  SCIPsetDebugMsg(set, "Benders' decomposition - solve loop %d\n", solveloop);
3155 
3156  if( type == SCIP_BENDERSENFOTYPE_CHECK && sol == NULL )
3157  {
3158  /* TODO: Check whether this is absolutely necessary. I think that this if statment can be removed. */
3159  locinfeasible = TRUE;
3160  }
3161  else
3162  {
3163  /* solving each of the subproblems for Benders' decomposition */
3164  /* TODO: ensure that the each of the subproblems solve and update the parameters with the correct return values
3165  */
3166 #ifndef __INTEL_COMPILER
3167  #pragma omp parallel for num_threads(numthreads) private(i) reduction(&&:locoptimal) reduction(||:locinfeasible) reduction(+:locnverified) reduction(||:locstopped) reduction(min:retcode)
3168 #endif
3169  for( j = 0; j < nsolveidx; j++ )
3170  {
3171  SCIP_Bool subinfeas = FALSE;
3172  SCIP_Bool convexsub;
3173  SCIP_Bool solvesub = TRUE;
3174  SCIP_Bool solved;
3175 
3176  i = solveidx[j];
3178 
3179  /* the subproblem is initially flagged as not solved for this solving loop */
3180  (*subprobsolved)[i] = FALSE;
3181 
3182  /* setting the subsystem status to UNKNOWN at the start of each solve loop */
3183  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_UNKNOWN;
3184 
3185  /* for the second solving loop, if the problem is an LP, it is not solved again. If the problem is a MIP,
3186  * then the subproblem objective function value is set to infinity. However, if the subproblem is proven
3187  * infeasible from the LP, then the IP loop is not performed.
3188  * If the solve loop is SCIP_BENDERSSOLVELOOP_USERCIP, then nothing is done. It is assumed that the user will
3189  * correctly update the objective function within the user-defined solving function.
3190  */
3191  if( solveloop == SCIP_BENDERSSOLVELOOP_CIP )
3192  {
3193  if( convexsub || (*substatus)[i] == SCIP_BENDERSSUBSTATUS_INFEAS )
3194  solvesub = FALSE;
3195  else
3196  {
3197  SCIPbendersSetSubproblemObjval(benders, i, SCIPbendersSubproblem(benders, i) != NULL ?
3198  SCIPinfinity(SCIPbendersSubproblem(benders, i)) : SCIPsetInfinity(set));
3199  }
3200  }
3201 
3202  /* if the subproblem is independent, then it does not need to be solved. In this case, the nverified flag will
3203  * increase by one. When the subproblem is not independent, then it needs to be checked.
3204  */
3205  if( !subproblemIsActive(benders, i) )
3206  {
3207  /* NOTE: There is no need to update the optimal flag. This is because optimal is always TRUE until a
3208  * non-optimal subproblem is found.
3209  */
3210  /* if the auxiliary variable value is infinity, then the subproblem has not been solved yet. Currently the
3211  * subproblem statue is unknown. */
3212  if( SCIPsetIsInfinity(set, SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i))
3213  || SCIPsetIsInfinity(set, -SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i))
3215  {
3216  SCIPbendersSetSubproblemObjval(benders, i, SCIPbendersSubproblem(benders, i) != NULL ?
3217  SCIPinfinity(SCIPbendersSubproblem(benders, i)) : SCIPsetInfinity(set));
3218 
3219  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_UNKNOWN;
3220  locoptimal = FALSE;
3221 
3222  SCIPsetDebugMsg(set, "Benders' decomposition: subproblem %d is not active, but has not been solved."
3223  " setting status to UNKNOWN\n", i);
3224  }
3225  else
3226  {
3228  SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i)) < benders->solutiontol )
3229  {
3230  SCIPbendersSetSubproblemObjval(benders, i, SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i));
3231  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_OPTIMAL;
3232  }
3233  else
3234  {
3236  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_AUXVIOL;
3237  }
3238 
3239  SCIPsetDebugMsg(set, "Benders' decomposition: subproblem %d is not active, setting status to OPTIMAL\n", i);
3240  }
3241 
3242  (*subprobsolved)[i] = TRUE;
3243 
3244  /* the nverified counter is only increased in the convex solving loop */
3245  if( solveloop == SCIP_BENDERSSOLVELOOP_CONVEX || solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX )
3246  locnverified++;
3247  }
3248  else if( solvesub )
3249  {
3250  retcode = SCIPbendersExecSubproblemSolve(benders, set, sol, i, solveloop, FALSE, &solved, &subinfeas, type);
3251 
3252  /* the solution for the subproblem is only processed if the return code is SCIP_OKAY */
3253  if( retcode == SCIP_OKAY )
3254  {
3255 #ifdef SCIP_DEBUG
3256  if( type == SCIP_BENDERSENFOTYPE_LP )
3257  {
3258  SCIPsetDebugMsg(set, "Enfo LP: Subproblem %d Type %d (%f < %f)\n", i,
3259  SCIPbendersGetSubproblemType(benders, i), SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i),
3260  SCIPbendersGetSubproblemObjval(benders, i));
3261  }
3262 #endif
3263  (*subprobsolved)[i] = solved;
3264 
3265  locinfeasible = locinfeasible || subinfeas;
3266  if( subinfeas )
3267  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_INFEAS;
3268 
3269  /* if the subproblems are solved to check integer feasibility, then the optimality check must be performed.
3270  * This will only be performed if checkint is TRUE and the subproblem was solved. The subproblem may not be
3271  * solved if the user has defined a solving function
3272  */
3273  if( checkint && (*subprobsolved)[i] )
3274  {
3275  /* if the subproblem is feasible, then it is necessary to update the value of the auxiliary variable to the
3276  * objective function value of the subproblem.
3277  */
3278  if( !subinfeas )
3279  {
3280  SCIP_Bool subproboptimal;
3281 
3282  subproboptimal = SCIPbendersSubproblemIsOptimal(benders, set, sol, i);
3283 
3284  if( subproboptimal )
3285  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_OPTIMAL;
3286  else
3287  (*substatus)[i] = SCIP_BENDERSSUBSTATUS_AUXVIOL;
3288 
3289  /* It is only possible to determine the optimality of a solution within a given subproblem in four
3290  * different cases:
3291  * i) solveloop == SCIP_BENDERSSOLVELOOP_CONVEX or USERCONVEX and the subproblem is convex.
3292  * ii) solveloop == SCIP_BENDERSOLVELOOP_CONVEX and only the convex relaxations will be checked.
3293  * iii) solveloop == SCIP_BENDERSSOLVELOOP_USERCIP and the subproblem was solved, since the user has
3294  * defined a solve function, it is expected that the solving is correctly executed.
3295  * iv) solveloop == SCIP_BENDERSSOLVELOOP_CIP and the MIP for the subproblem has been solved.
3296  */
3297  if( convexsub || onlyconvexcheck
3298  || solveloop == SCIP_BENDERSSOLVELOOP_CIP
3299  || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP )
3300  locoptimal = locoptimal && subproboptimal;
3301 
3302 #ifdef SCIP_DEBUG
3303  if( convexsub || solveloop >= SCIP_BENDERSSOLVELOOP_CIP )
3304  {
3305  if( subproboptimal )
3306  {
3307  SCIPsetDebugMsg(set, "Subproblem %d is Optimal (%f >= %f)\n", i,
3308  SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i), SCIPbendersGetSubproblemObjval(benders, i));
3309  }
3310  else
3311  {
3312  SCIPsetDebugMsg(set, "Subproblem %d is NOT Optimal (%f < %f)\n", i,
3313  SCIPbendersGetAuxiliaryVarVal(benders, set, sol, i), SCIPbendersGetSubproblemObjval(benders, i));
3314  }
3315  }
3316 #endif
3317 
3318  /* the nverified variable is only incremented when the original form of the subproblem has been solved.
3319  * What is meant by "original" is that the LP relaxation of CIPs are solved to generate valid cuts. So
3320  * if the subproblem is defined as a CIP, then it is only classified as checked if the CIP is solved.
3321  * There are three cases where the "original" form is solved are:
3322  * i) solveloop == SCIP_BENDERSSOLVELOOP_CONVEX or USERCONVEX and the subproblem is an LP
3323  * - the original form has been solved.
3324  * ii) solveloop == SCIP_BENDERSSOLVELOOP_CIP or USERCIP and the CIP for the subproblem has been
3325  * solved.
3326  * iii) or, only a convex check is performed.
3327  */
3328  if( ((solveloop == SCIP_BENDERSSOLVELOOP_CONVEX || solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX)
3329  && convexsub)
3330  || ((solveloop == SCIP_BENDERSSOLVELOOP_CIP || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP)
3331  && !convexsub)
3332  || onlyconvexcheck )
3333  locnverified++;
3334  }
3335  }
3336  }
3337  }
3338 
3339  /* checking whether the limits have been exceeded in the master problem */
3340  locstopped = SCIPisStopped(set->scip);
3341  }
3342  }
3343 
3344  /* setting the input parameters to the local variables */
3345  SCIPsetDebugMsg(set, "Local variable values: nverified %d infeasible %u optimal %u stopped %u\n", locnverified,
3346  locinfeasible, locoptimal, locstopped);
3347  *nverified = locnverified;
3348  *infeasible = locinfeasible;
3349  *optimal = locoptimal;
3350  *stopped = locstopped;
3351 
3352  return retcode;
3353 }
3354 
3355 /** Calls the Benders' decompsition cuts for the given solve loop. There are four cases:
3356  * i) solveloop == SCIP_BENDERSSOLVELOOP_CONVEX - only the LP Benders' cuts are called
3357  * ii) solveloop == SCIP_BENDERSSOLVELOOP_CIP - only the CIP Benders' cuts are called
3358  * iii) solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX - only the LP Benders' cuts are called
3359  * iv) solveloop == SCIP_BENDERSSOLVELOOP_USERCIP - only the CIP Benders' cuts are called
3360  *
3361  * The priority of the results are: SCIP_CONSADDED (SCIP_SEPARATED), SCIP_DIDNOTFIND, SCIP_FEASIBLE, SCIP_DIDNOTRUN. In
3362  * this function, there are four levels of results that need to be assessed. These are:
3363  * i) The result from the individual cut for the subproblem
3364  * ii) The overall result for the subproblem from all cuts
3365  * iii) the overall result for the solve loop from all cuts
3366  * iv) the over all result from all solve loops.
3367  * In each level, the priority of results must be adhered to.
3368  */
3369 static
3371  SCIP_BENDERS* benders, /**< Benders' decomposition */
3372  SCIP_SET* set, /**< global SCIP settings */
3373  SCIP_SOL* sol, /**< primal CIP solution */
3374  SCIP_RESULT* result, /**< result of the pricing process */
3375  SCIP_BENDERSENFOTYPE type, /**< the type of solution being enforced */
3376  SCIP_BENDERSSOLVELOOP solveloop, /**< the current solve loop */
3377  SCIP_Bool checkint, /**< are the subproblems called during a check/enforce of integer sols? */
3378  SCIP_Bool* subprobsolved, /**< an array indicating the subproblems that were solved in this loop. */
3379  SCIP_BENDERSSUBSTATUS* substatus, /**< array to store the status of the subsystem */
3380  int* solveidx, /**< the indices of subproblems to be solved in this loop */
3381  int nsolveidx, /**< the number of subproblems to be solved in this loop */
3382  int** mergecands, /**< the subproblems that are merge candidates */
3383  int* npriomergecands, /**< the number of priority merge candidates. */
3384  int* nmergecands, /**< the number of merge candidates. */
3385  int* nsolveloops /**< the number of solve loops, is updated w.r.t added cuts */
3386  )
3387 {
3388  SCIP_BENDERSCUT** benderscuts;
3389  SCIP_RESULT solveloopresult;
3390  int nbenderscuts;
3391  SCIP_Longint addedcuts = 0;
3392  int i;
3393  int j;
3394  int k;
3395  SCIP_Bool onlyconvexcheck;
3396 
3397  assert(benders != NULL);
3398  assert(set != NULL);
3399 
3400  /* getting the Benders' decomposition cuts */
3401  benderscuts = SCIPbendersGetBenderscuts(benders);
3402  nbenderscuts = SCIPbendersGetNBenderscuts(benders);
3403 
3404  solveloopresult = SCIP_DIDNOTRUN;
3405 
3406  /* in the case of an LNS check, only the convex relaxations of the subproblems will be solved. This is a performance
3407  * feature, since solving the convex relaxation is typically much faster than solving the corresponding CIP. While
3408  * the CIP is not solved during the LNS check, the solutions are still of higher quality than when Benders' is not
3409  * employed.
3410  */
3411  onlyconvexcheck = SCIPbendersOnlyCheckConvexRelax(benders, SCIPsetGetSubscipsOff(set));
3412 
3413  /* It is only possible to add cuts to the problem if it has not already been solved */
3416  && (benders->cutcheck || type != SCIP_BENDERSENFOTYPE_CHECK) )
3417  {
3418  /* This is done in two loops. The first is by subproblem and the second is by cut type. */
3419  for( k = 0; k < nsolveidx; k++ )
3420  {
3421  SCIP_RESULT subprobresult;
3422  SCIP_Bool convexsub;
3423 
3424  i = solveidx[k];
3425 
3427 
3428  /* cuts can only be generated if the subproblem is not independent and if it has been solved. Additionally, the
3429  * status of the subproblem solving must not be INFEASIBLE while in a cut strengthening round.
3430  * The subproblem solved flag is important for the user-defined subproblem solving methods
3431  */
3432  if( subproblemIsActive(benders, i) && subprobsolved[i]
3433  && !(substatus[i] == SCIP_BENDERSSUBSTATUS_INFEAS && benders->strengthenround) )
3434  {
3435  subprobresult = SCIP_DIDNOTRUN;
3436  for( j = 0; j < nbenderscuts; j++ )
3437  {
3438  SCIP_RESULT cutresult;
3439  SCIP_Longint prevaddedcuts;
3440 
3441  assert(benderscuts[j] != NULL);
3442 
3443  prevaddedcuts = SCIPbenderscutGetNFound(benderscuts[j]);
3444  cutresult = SCIP_DIDNOTRUN;
3445 
3446  /* the result is updated only if a Benders' cut is generated or one was not found. However, if a cut has
3447  * been found in a previous iteration, then the result is returned as SCIP_CONSADDED or SCIP_SEPARATED.
3448  * This result is permitted because if a constraint was added, the solution that caused the error in the cut
3449  * generation will be cutoff from the master problem.
3450  */
3451  if( (SCIPbenderscutIsLPCut(benderscuts[j]) && (solveloop == SCIP_BENDERSSOLVELOOP_CONVEX
3452  || solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX))
3453  || (!SCIPbenderscutIsLPCut(benderscuts[j]) && ((solveloop == SCIP_BENDERSSOLVELOOP_CIP && !convexsub)
3454  || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP)) )
3455  SCIP_CALL( SCIPbenderscutExec(benderscuts[j], set, benders, sol, i, type, &cutresult) );
3456 
3457  addedcuts += (SCIPbenderscutGetNFound(benderscuts[j]) - prevaddedcuts);
3458 
3459  /* the result is updated only if a Benders' cut is generated */
3460  if( cutresult == SCIP_CONSADDED || cutresult == SCIP_SEPARATED )
3461  {
3462  subprobresult = cutresult;
3463 
3464  benders->ncutsfound++;
3465 
3466  /* at most a single cut is generated for each subproblem */
3467  break;
3468  }
3469  else
3470  {
3471  /* checking from lowest priority result */
3472  if( subprobresult == SCIP_DIDNOTRUN )
3473  subprobresult = cutresult;
3474  else if( subprobresult == SCIP_FEASIBLE && cutresult == SCIP_DIDNOTFIND )
3475  subprobresult = cutresult;
3476  /* if the subprobresult is SCIP_DIDNOTFIND, then it can't be updated. */
3477  }
3478  }
3479 
3480  /* the highest priority for the results is CONSADDED and SEPARATED. The solveloopresult will always be
3481  * updated if the subprobresult is either of these.
3482  */
3483  if( subprobresult == SCIP_CONSADDED || subprobresult == SCIP_SEPARATED )
3484  {
3485  solveloopresult = subprobresult;
3486  }
3487  else if( subprobresult == SCIP_FEASIBLE )
3488  {
3489  /* updating the solve loop result based upon the priority */
3490  if( solveloopresult == SCIP_DIDNOTRUN )
3491  solveloopresult = subprobresult;
3492  }
3493  else if( subprobresult == SCIP_DIDNOTFIND )
3494  {
3495  /* updating the solve loop result based upon the priority */
3496  if( solveloopresult == SCIP_DIDNOTRUN || solveloopresult == SCIP_FEASIBLE )
3497  solveloopresult = subprobresult;
3498 
3499  /* since a cut was not found, then merging could be useful to avoid this in subsequent iterations. The
3500  * candidate is labelled as a non-priority merge candidate
3501  */
3502  if( substatus[i] != SCIP_BENDERSSUBSTATUS_OPTIMAL )
3503  {
3504  (*mergecands)[(*nmergecands)] = i;
3505  (*nmergecands)++;
3506  }
3507  }
3508  else if( subprobresult == SCIP_DIDNOTRUN )
3509  {
3510  /* if the subproblem is infeasible and no cut generation methods were run, then the infeasibility will
3511  * never be resolved. As such, the subproblem will be merged into the master problem. If the subproblem
3512  * was not infeasible, then it is added as a possible merge candidate
3513  */
3514  if( substatus[i] == SCIP_BENDERSSUBSTATUS_INFEAS )
3515  {
3516  (*mergecands)[(*nmergecands)] = (*mergecands)[(*npriomergecands)];
3517  (*mergecands)[(*npriomergecands)] = i;
3518  (*npriomergecands)++;
3519  (*nmergecands)++;
3520  }
3521  else if( substatus[i] != SCIP_BENDERSSUBSTATUS_OPTIMAL )
3522  {
3523  (*mergecands)[(*nmergecands)] = i;
3524  (*nmergecands)++;
3525  }
3526  }
3527  }
3528  }
3529  }
3530 
3531  /* updating the overall result based upon the priorities */
3532  if( solveloopresult == SCIP_CONSADDED || solveloopresult == SCIP_SEPARATED )
3533  {
3534  (*result) = solveloopresult;
3535  }
3536  else if( solveloopresult == SCIP_FEASIBLE )
3537  {
3538  /* updating the solve loop result based upon the priority */
3539  if( (*result) == SCIP_DIDNOTRUN )
3540  (*result) = solveloopresult;
3541  }
3542  else if( solveloopresult == SCIP_DIDNOTFIND )
3543  {
3544  /* updating the solve loop result based upon the priority */
3545  if( (*result) == SCIP_DIDNOTRUN || (*result) == SCIP_FEASIBLE )
3546  (*result) = solveloopresult;
3547  }
3548 
3549  /* if no cuts were added, then the number of solve loops is increased */
3550  if( addedcuts == 0 && SCIPbendersGetNConvexSubproblems(benders) < SCIPbendersGetNSubproblems(benders)
3551  && checkint && !onlyconvexcheck )
3552  (*nsolveloops) = 2;
3553 
3554  return SCIP_OKAY;
3555 }
3556 
3557 /** Solves the subproblem using the current master problem solution.
3558  *
3559  * The checkint flag indicates whether integer feasibility can be assumed. If it is not assumed, i.e. checkint ==
3560  * FALSE, then only the convex relaxations of the subproblems are solved. If integer feasibility is assumed, i.e.
3561  * checkint == TRUE, then the convex relaxations and the full CIP are solved to generate Benders' cuts and check
3562  * solution feasibility.
3563  *
3564  * TODO: consider allowing the possibility to pass solution information back from the subproblems instead of the scip
3565  * instance. This would allow the use of different solvers for the subproblems, more importantly allowing the use of an
3566  * LP solver for LP subproblems.
3567  */
3569  SCIP_BENDERS* benders, /**< Benders' decomposition */
3570  SCIP_SET* set, /**< global SCIP settings */
3571  SCIP_SOL* sol, /**< primal CIP solution */
3572  SCIP_RESULT* result, /**< result of the pricing process */
3573  SCIP_Bool* infeasible, /**< is the master problem infeasible with respect to the Benders' cuts? */
3574  SCIP_Bool* auxviol, /**< set to TRUE only if the solution is feasible but the aux vars are violated */
3575  SCIP_BENDERSENFOTYPE type, /**< the type of solution being enforced */
3576  SCIP_Bool checkint /**< should the integer solution be checked by the subproblems */
3577  )
3578 {
3579  int nsubproblems;
3580  int subproblemcount;
3581  int nsolveloops;
3582  int nverified;
3583  int nsolved;
3584  int* mergecands;
3585  int npriomergecands;
3586  int nmergecands;
3587  int* solveidx;
3588  int* executedidx;
3589  int nsolveidx;
3590  int nexecutedidx;
3591  int nfree;
3592  SCIP_Bool* subprobsolved;
3593  SCIP_BENDERSSUBSTATUS* substatus;
3594  SCIP_Bool optimal;
3595  SCIP_Bool allverified;
3596  SCIP_Bool success;
3597  SCIP_Bool stopped;
3598  int i;
3599  int l;
3600 
3601  success = TRUE;
3602  stopped = FALSE;
3603 
3604  SCIPsetDebugMsg(set, "Starting Benders' decomposition subproblem solving. type %d checkint %u\n", type, checkint);
3605 
3606 #ifdef SCIP_MOREDEBUG
3607  SCIP_CALL( SCIPprintSol(set->scip, sol, NULL, FALSE) );
3608 #endif
3609 
3610  /* start timing */
3611  SCIPclockStart(benders->bendersclock, set);
3612 
3613  nsubproblems = SCIPbendersGetNSubproblems(benders);
3614 
3615  (*auxviol) = FALSE;
3616  (*infeasible) = FALSE;
3617 
3618  /* It is assumed that the problem is optimal, until a subproblem is found not to be optimal. However, not all
3619  * subproblems could be checked in each iteration. As such, it is not possible to state that the problem is optimal
3620  * if not all subproblems are checked. Situations where this may occur is when a subproblem is a MIP and only the LP
3621  * is solved. Also, in a distributed computation, then it may be advantageous to only solve some subproblems before
3622  * resolving the master problem. As such, for a problem to be optimal, then (optimal && allverified) == TRUE
3623  */
3624  optimal = TRUE;
3625  nverified = 0;
3626  nsolved = 0;
3627 
3628  assert(benders != NULL);
3629  assert(result != NULL);
3630  assert(infeasible != NULL);
3631  assert(auxviol != NULL);
3632 
3633  /* if the Benders' decomposition is called from a sub-SCIP and the sub-SCIPs have been deactivated, then it is
3634  * assumed that this is an LNS heuristic. As such, the check is not performed and the solution is assumed to be
3635  * feasible
3636  */
3637  if( benders->iscopy && set->subscipsoff
3638  && (!benders->lnscheck
3639  || (benders->lnsmaxdepth > -1 && SCIPgetDepth(benders->sourcescip) >= benders->lnsmaxdepth)
3640  || (benders->lnsmaxcalls > -1 && SCIPbendersGetNCalls(benders) >= benders->lnsmaxcalls)
3641  || (type != SCIP_BENDERSENFOTYPE_CHECK && SCIPgetDepth(set->scip) == 0 && benders->lnsmaxcallsroot > -1
3642  && SCIPbendersGetNCalls(benders) >= benders->lnsmaxcallsroot)) )
3643  {
3644  (*result) = SCIP_DIDNOTRUN;
3645  return SCIP_OKAY;
3646  }
3647 
3648  /* it is not necessary to check all primal solutions by solving the Benders' decomposition subproblems.
3649  * Only the improving solutions are checked to improve efficiency of the algorithm.
3650  * If the solution is non-improving, the result FEASIBLE is returned. While this may be incorrect w.r.t to the
3651  * Benders' subproblems, this solution will never be the optimal solution. A non-improving solution may be used
3652  * within LNS primal heuristics. If this occurs, the improving solution, if found, will be checked by the solving
3653  * the Benders' decomposition subproblems.
3654  * TODO: Add a parameter to control this behaviour.
3655  */
3656  if( checkint && SCIPsetIsLE(set, SCIPgetPrimalbound(set->scip)*(int)SCIPgetObjsense(set->scip),
3657  SCIPgetSolOrigObj(set->scip, sol)*(int)SCIPgetObjsense(set->scip)) )
3658  {
3659  (*result) = SCIP_DIDNOTRUN;
3660  return SCIP_OKAY;
3661  }
3662 
3663  /* if the enforcement type is SCIP_BENDERSENFOTYPE_LP and the LP is currently unbounded. This could mean that there
3664  * is no lower bound on the auxiliary variables. In this case, we try to update the lower bound for the auxiliary
3665  * variables.
3666  */
3668  && benders->updateauxvarbound )
3669  {
3670  SCIP_CALL( updateAuxiliaryVarLowerbound(benders, set, result) );
3671 
3672  /* the auxiliary variable bound will only be updated once. */
3673  benders->updateauxvarbound = FALSE;
3674  }
3675 
3676  /* sets the stored objective function values of the subproblems to infinity */
3677  resetSubproblemObjectiveValue(benders, set);
3678 
3679  *result = SCIP_DIDNOTRUN;
3680 
3681  if( benders->benderspresubsolve != NULL && !benders->strengthenround )
3682  {
3683  SCIP_Bool skipsolve;
3684 
3685  skipsolve = FALSE;
3686  SCIP_CALL( benders->benderspresubsolve(set->scip, benders, sol, type, checkint, infeasible, auxviol, &skipsolve,
3687  result) );
3688 
3689  /* evaluate result */
3690  if( (*result) != SCIP_DIDNOTRUN
3691  && (*result) != SCIP_FEASIBLE
3692  && (*result) != SCIP_INFEASIBLE
3693  && (*result) != SCIP_CONSADDED
3694  && (*result) != SCIP_SEPARATED )
3695  {
3696  SCIPerrorMessage("the user-defined pre subproblem solving method for the Benders' decomposition <%s> returned "
3697  "invalid result <%d>\n", benders->name, *result);
3698  return SCIP_INVALIDRESULT;
3699  }
3700 
3701  /* if the solve must be skipped, then the solving loop is exited and the user defined result is returned */
3702  if( skipsolve )
3703  {
3704  SCIPsetDebugMsg(set, "skipping the subproblem solving for Benders' decomposition <%s>. "
3705  "returning result <%d>\n", benders->name, *result);
3706  return SCIP_OKAY;
3707  }
3708  }
3709 
3710  /* the cut strengthening is performed before the regular subproblem solve is called. To avoid recursion, the flag
3711  * strengthenround is set to TRUE when the cut strengthening is performed. The cut strengthening is not performed as
3712  * part of the large neighbourhood Benders' search.
3713  *
3714  * NOTE: cut strengthening is only applied for fractional solutions and integer solutions if there are no CIP
3715  * subproblems.
3716  */
3717  if( benders->strengthenenabled && !benders->strengthenround && !benders->iscopy
3718  && (!checkint || SCIPbendersGetNConvexSubproblems(benders) == SCIPbendersGetNSubproblems(benders)) )
3719  {
3720  SCIP_Bool skipsolve;
3721 
3722  benders->strengthenround = TRUE;
3723  /* if the user has not requested the solve to be skipped, then the cut strengthening is performed */
3724  SCIP_CALL( performInteriorSolCutStrengthening(benders, set, sol, type, checkint, FALSE, infeasible, auxviol,
3725  &skipsolve, result) );
3726  benders->strengthenround = FALSE;
3727 
3728  /* if the solve must be skipped, then the solving loop is exited and the user defined result is returned */
3729  if( skipsolve )
3730  {
3731  SCIPsetDebugMsg(set, "skipping the subproblem solving because cut strengthening found a cut "
3732  "for Benders' decomposition <%s>. Returning result <%d>\n", benders->name, *result);
3733  return SCIP_OKAY;
3734  }
3735 
3736  /* the result flag need to be reset to DIDNOTRUN for the main subproblem solve */
3737  (*result) = SCIP_DIDNOTRUN;
3738  }
3739 
3740  /* allocating memory for the infeasible subproblem array */
3741  SCIP_CALL( SCIPallocClearBlockMemoryArray(set->scip, &subprobsolved, nsubproblems) );
3742  SCIP_CALL( SCIPallocClearBlockMemoryArray(set->scip, &substatus, nsubproblems) );
3743  SCIP_CALL( SCIPallocClearBlockMemoryArray(set->scip, &mergecands, nsubproblems) );
3744  npriomergecands = 0;
3745  nmergecands = 0;
3746 
3747  /* allocating the memory for the subproblem solving and cut generation indices */
3748  SCIP_CALL( SCIPallocClearBlockMemoryArray(set->scip, &solveidx, nsubproblems) );
3749  SCIP_CALL( SCIPallocClearBlockMemoryArray(set->scip, &executedidx, nsubproblems) );
3750  nsolveidx = 0;
3751  nexecutedidx = 0;
3752 
3753  /* only a subset of the subproblems are initially solved. Both solving loops are executed for the subproblems to
3754  * check whether any cuts are generated. If a cut is generated, then no further subproblems are solved. If a cut is
3755  * not generated, then an additional set of subproblems are solved.
3756  */
3757  while( nsolved < nsubproblems )
3758  {
3759  /* getting the indices for the subproblems that will be solved */
3760  createSolveSubproblemIndexList(benders, set, type, &solveidx, &nsolveidx);
3761 
3762  /* by default the number of solve loops is 1. This is the case if all subproblems are LP or the user has defined a
3763  * benderssolvesub callback. If there is a subproblem that is not an LP, then 2 solve loops are performed. The first
3764  * loop is the LP solving loop, the second solves the subproblem to integer optimality.
3765  */
3766  nsolveloops = 1;
3767 
3768  for( l = 0; l < nsolveloops; l++ )
3769  {
3770  SCIP_BENDERSSOLVELOOP solveloop; /* identifies what problem type is solve in this solve loop */
3771 
3772  /* if either benderssolvesubconvex or benderssolvesub are implemented, then the user callbacks are invoked */
3773  if( benders->benderssolvesubconvex != NULL || benders->benderssolvesub != NULL )
3774  {
3775  if( l == 0 )
3777  else
3778  solveloop = SCIP_BENDERSSOLVELOOP_USERCIP;
3779  }
3780  else
3781  solveloop = (SCIP_BENDERSSOLVELOOP) l;
3782 
3783  /* solving the subproblems for this round of enforcement/checking. */
3784  SCIP_CALL( solveBendersSubproblems(benders, set, sol, type, solveloop, checkint, &nverified,
3785  solveidx, nsolveidx, &subprobsolved, &substatus, infeasible, &optimal, &stopped) );
3786 
3787  /* if the solving has been stopped, then the subproblem solving and cut generation must terminate */
3788  if( stopped )
3789  break;
3790 
3791  /* Generating cuts for the subproblems. Cuts are only generated when the solution is from primal heuristics,
3792  * relaxations or the LP
3793  */
3794  if( type != SCIP_BENDERSENFOTYPE_PSEUDO )
3795  {
3796  SCIP_CALL( generateBendersCuts(benders, set, sol, result, type, solveloop, checkint, subprobsolved,
3797  substatus, solveidx, nsolveidx, &mergecands, &npriomergecands, &nmergecands, &nsolveloops) );
3798  }
3799  else
3800  {
3801  /* The first solving loop solves the convex subproblems and the convex relaxations of the CIP subproblems. The
3802  * second solving loop solves the CIP subproblems. The second solving loop is only called if the integer
3803  * feasibility is being checked and if the convex subproblems and convex relaxations are not infeasible.
3804  */
3805  if( !(*infeasible) && checkint && !SCIPbendersOnlyCheckConvexRelax(benders, SCIPsetGetSubscipsOff(set))
3807  nsolveloops = 2;
3808  }
3809  }
3810 
3811  nsolved += nsolveidx;
3812 
3813  /* storing the indices of the subproblems for which the solving loop was executed */
3814  for( i = 0; i < nsolveidx; i++ )
3815  executedidx[nexecutedidx++] = solveidx[i];
3816 
3817  /* if the result is CONSADDED or SEPARATED, then a cut is generated and no further subproblem processing is
3818  * required
3819  */
3820  if( (*result) == SCIP_CONSADDED || (*result) == SCIP_SEPARATED )
3821  break;
3822  }
3823 
3824  /* inserting the subproblems into the priority queue for the next solve call */
3825  SCIP_CALL( updateSubproblemStatQueue(benders, executedidx, nexecutedidx, TRUE) );
3826 
3827  if( stopped )
3828  goto TERMINATE;
3829 
3830  allverified = (nverified == nsubproblems);
3831 
3832  SCIPsetDebugMsg(set, "End Benders' decomposition subproblem solve. result %d infeasible %u auxviol %u nverified %d\n",
3833  *result, *infeasible, *auxviol, nverified);
3834 
3835 #ifdef SCIP_DEBUG
3836  if( (*result) == SCIP_CONSADDED )
3837  {
3838  SCIPsetDebugMsg(set, "Benders' decomposition: Cut added\n");
3839  }
3840 #endif
3841 
3842  /* if the number of checked pseudo solutions exceeds a set limit, then all subproblems are passed as merge
3843  * candidates. Currently, merging subproblems into the master problem is the only method for resolving numerical
3844  * troubles.
3845  *
3846  * We are only interested in the pseudo solutions that have been checked completely for integrality. This is
3847  * identified by checkint == TRUE. This means that the Benders' decomposition constraint is one of the last
3848  * constraint handlers that must resolve the infeasibility. If the Benders' decomposition framework can't resolve the
3849  * infeasibility, then this will result in an error.
3850  */
3851  if( type == SCIP_BENDERSENFOTYPE_PSEUDO && checkint )
3852  {
3853  benders->npseudosols++;
3854 
3855  if( benders->npseudosols > BENDERS_MAXPSEUDOSOLS )
3856  {
3857  /* if a priority merge candidate already exists, then no other merge candidates need to be added.*/
3858  if( npriomergecands == 0 )
3859  {
3860  /* all subproblems are added to the merge candidate list. The first active subproblem is added as a
3861  * priority merge candidate
3862  */
3863  nmergecands = 0;
3864  npriomergecands = 1;
3865  for( i = 0; i < nsubproblems; i++ )
3866  {
3867  /* only active subproblems are added to the merge candidate list */
3868  if( subproblemIsActive(benders, i) )
3869  {
3870  mergecands[nmergecands] = i;
3871  nmergecands++;
3872  }
3873  }
3874 
3875  SCIPverbMessage(set->scip, SCIP_VERBLEVEL_HIGH, NULL, " The number of checked pseudo solutions exceeds the "
3876  "limit of %d. All active subproblems are merge candidates, with subproblem %d a priority candidate.\n",
3877  BENDERS_MAXPSEUDOSOLS, mergecands[0]);
3878  }
3879  }
3880  }
3881  else
3882  benders->npseudosols = 0;
3883 
3884  /* if the result is SCIP_DIDNOTFIND, then there was a error in generating cuts in all subproblems that are not
3885  * optimal. This result does not cutoff any solution, so the Benders' decomposition algorithm will fail.
3886  *
3887  * It could happen that the cut strengthening approach causes an error the cut generation. In this case, an error
3888  * should not be thrown. So, this check will be skipped when in a strengthening round.
3889  * TODO: Work out a way to ensure Benders' decomposition does not terminate due to a SCIP_DIDNOTFIND result.
3890  */
3891  if( (*result) == SCIP_DIDNOTFIND && !benders->strengthenround )
3892  {
3893  if( type == SCIP_BENDERSENFOTYPE_PSEUDO )
3894  (*result) = SCIP_SOLVELP;
3895  else
3896  (*result) = SCIP_INFEASIBLE;
3897 
3898  SCIPerrorMessage("An error was found when generating cuts for non-optimal subproblems of Benders' "
3899  "decomposition <%s>. Consider merging the infeasible subproblems into the master problem.\n", SCIPbendersGetName(benders));
3900 
3901  /* since no other cuts are generated, then this error will result in a crash. It is possible to avoid the error,
3902  * by merging the affected subproblem into the master problem.
3903  *
3904  * NOTE: If the error occurs while checking solutions, i.e. SCIP_BENDERSENFOTYPE_CHECK, then it is valid to set
3905  * the result to SCIP_INFEASIBLE and the success flag to TRUE
3906  */
3907  if( type != SCIP_BENDERSENFOTYPE_CHECK )
3908  success = FALSE;
3909 
3910  goto POSTSOLVE;
3911  }
3912 
3913  if( type == SCIP_BENDERSENFOTYPE_PSEUDO )
3914  {
3915  if( (*infeasible) || !allverified )
3916  (*result) = SCIP_SOLVELP;
3917  else
3918  {
3919  (*result) = SCIP_FEASIBLE;
3920 
3921  /* if the subproblems are not infeasible, but they are also not optimal. This means that there is a violation
3922  * in the auxiliary variable values. In this case, a feasible result is returned with the auxviol flag set to
3923  * TRUE.
3924  */
3925  (*auxviol) = !optimal;
3926  }
3927  }
3928  else if( checkint && (type == SCIP_BENDERSENFOTYPE_CHECK
3929  || ((*result) != SCIP_CONSADDED && (*result) != SCIP_SEPARATED)) )
3930  {
3931  /* if the subproblems are being solved as part of conscheck, then the results flag must be returned after the solving
3932  * has completed.
3933  */
3934  if( (*infeasible) || !allverified )
3935  (*result) = SCIP_INFEASIBLE;
3936  else
3937  {
3938  (*result) = SCIP_FEASIBLE;
3939 
3940  /* if the subproblems are not infeasible, but they are also not optimal. This means that there is a violation
3941  * in the auxiliary variable values. In this case, a feasible result is returned with the auxviol flag set to
3942  * TRUE.
3943  */
3944  (*auxviol) = !optimal;
3945  }
3946  }
3947 
3948 POSTSOLVE:
3949  /* calling the post-solve call back for the Benders' decomposition algorithm. This allows the user to work directly
3950  * with the solved subproblems and the master problem */
3951  if( benders->benderspostsolve != NULL )
3952  {
3953  SCIP_Bool merged;
3954 
3955  merged = FALSE;
3956 
3957  SCIP_CALL( benders->benderspostsolve(set->scip, benders, sol, type, mergecands, npriomergecands, nmergecands,
3958  checkint, (*infeasible), &merged) );
3959 
3960  if( merged )
3961  {
3962  (*result) = SCIP_CONSADDED;
3963 
3964  /* since subproblems have been merged, then constraints have been added. This could resolve the unresolved
3965  * infeasibility, so the error has been corrected.
3966  */
3967  success = TRUE;
3968  }
3969  else if( !success )
3970  {
3971  SCIPerrorMessage("An error occurred during Benders' decomposition cut generations and no merging had been "
3972  "performed. It is not possible to continue solving the problem by Benders' decomposition\n");
3973  }
3974  }
3975 
3976 TERMINATE:
3977  /* if the solving process has stopped, then all subproblems need to be freed */
3978  if( stopped )
3979  nfree = nsubproblems;
3980  else
3981  nfree = nexecutedidx;
3982 
3983  /* freeing the subproblems after the cuts are generated */
3984  subproblemcount = 0;
3985  while( subproblemcount < nfree )
3986  {
3987  int subidx;
3988 
3989  if( stopped )
3990  subidx = subproblemcount;
3991  else
3992  subidx = executedidx[subproblemcount];
3993 
3994  SCIP_CALL( SCIPbendersFreeSubproblem(benders, set, subidx) );
3995 
3996  subproblemcount++;
3997  }
3998 
3999 #ifndef NDEBUG
4000  for( i = 0; i < nsubproblems; i++ )
4001  assert(SCIPbendersSubproblem(benders, i) == NULL
4003  || !SCIPinProbing(SCIPbendersSubproblem(benders, i))
4004  || !subproblemIsActive(benders, i));
4005 #endif
4006 
4007  /* increment the number of calls to the Benders' decomposition subproblem solve */
4008  benders->ncalls++;
4009 
4010  SCIPsetDebugMsg(set, "End Benders' decomposition execution method. result %d infeasible %u auxviol %u\n", *result,
4011  *infeasible, *auxviol);
4012 
4013  /* end timing */
4014  SCIPclockStop(benders->bendersclock, set);
4015 
4016  /* freeing memory */
4017  SCIPfreeBlockMemoryArray(set->scip, &executedidx, nsubproblems);
4018  SCIPfreeBlockMemoryArray(set->scip, &solveidx, nsubproblems);
4019  SCIPfreeBlockMemoryArray(set->scip, &mergecands, nsubproblems);
4020  SCIPfreeBlockMemoryArray(set->scip, &substatus, nsubproblems);
4021  SCIPfreeBlockMemoryArray(set->scip, &subprobsolved, nsubproblems);
4022 
4023  /* if there was an error in generating cuts and merging was not performed, then the solution is perturbed in an
4024  * attempt to generate a cut and correct the infeasibility
4025  */
4026  if( !success && !stopped )
4027  {
4028  SCIP_Bool skipsolve;
4029  SCIP_RESULT perturbresult;
4030 
4031  skipsolve = FALSE;
4032 
4033  benders->strengthenround = TRUE;
4034  /* if the user has not requested the solve to be skipped, then the cut strengthening is performed */
4035  SCIP_CALL( performInteriorSolCutStrengthening(benders, set, sol, type, checkint, TRUE, infeasible, auxviol,
4036  &skipsolve, &perturbresult) );
4037  benders->strengthenround = FALSE;
4038 
4039  if( perturbresult == SCIP_CONSADDED || perturbresult == SCIP_SEPARATED )
4040  (*result) = perturbresult;
4041 
4042  success = skipsolve;
4043  }
4044 
4045  /* if the Benders' decomposition subproblem check stopped, then we don't have a valid result. In this case, the
4046  * safest thing to do is report INFEASIBLE.
4047  */
4048  if( stopped )
4049  (*result) = SCIP_INFEASIBLE;
4050 
4051  /* if the subproblem verification identifies the solution as feasible, then a check whether slack variables have been
4052  * used is necessary. If any slack variables are non-zero, then the solution is reverified after the objective
4053  * coefficient for the slack variables is increased.
4054  */
4055  if( (*result) == SCIP_FEASIBLE )
4056  {
4057  SCIP_Bool activeslack;
4058 
4059  SCIP_CALL( SCIPbendersSolSlackVarsActive(benders, &activeslack) );
4060  SCIPsetDebugMsg(set, "Type: %d Active slack: %u Feasibility Phase: %u\n", type, activeslack,
4061  benders->feasibilityphase);
4062  if( activeslack )
4063  {
4064  if( type == SCIP_BENDERSENFOTYPE_CHECK )
4065  (*result) = SCIP_INFEASIBLE;
4066  else
4067  {
4068  /* increasing the value of the slack variable by a factor of 10 */
4069  benders->slackvarcoef *= 10.0;
4070 
4071  if( benders->slackvarcoef <= benders->maxslackvarcoef )
4072  {
4073  SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(set->scip), set->disp_verblevel, SCIP_VERBLEVEL_HIGH, "Increasing the slack variable coefficient to %g.\n", benders->slackvarcoef);
4074  }
4075  else
4076  {
4077  SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(set->scip), set->disp_verblevel, SCIP_VERBLEVEL_HIGH, "Fixing the slack variables to zero.\n");
4078  }
4079 
4080  /* resolving the subproblems with an increased slack variable */
4081  SCIP_CALL( SCIPsolveBendersSubproblems(set->scip, benders, sol, result, infeasible, auxviol, type, checkint) );
4082  }
4083  }
4084  else if( benders->feasibilityphase )
4085  {
4086  if( type != SCIP_BENDERSENFOTYPE_CHECK )
4087  {
4088  /* disabling the feasibility phase */
4089  benders->feasibilityphase = FALSE;
4090 
4091  /* resolving the subproblems with the slack variables set to zero */
4092  SCIP_CALL( SCIPsolveBendersSubproblems(set->scip, benders, sol, result, infeasible, auxviol, type, checkint) );
4093  }
4094  }
4095  }
4096 
4097  if( !success )
4098  return SCIP_ERROR;
4099  else
4100  return SCIP_OKAY;
4101 }
4102 
4103 /** solves the user-defined subproblem solving function */
4104 static
4106  SCIP_BENDERS* benders, /**< Benders' decomposition */
4107  SCIP_SET* set, /**< global SCIP settings */
4108  SCIP_SOL* sol, /**< primal CIP solution */
4109  int probnumber, /**< the subproblem number */
4110  SCIP_BENDERSSOLVELOOP solveloop, /**< the solve loop iteration. The first iter is for LP, the second for IP */
4111  SCIP_Bool* infeasible, /**< returns whether the current subproblem is infeasible */
4112  SCIP_Real* objective, /**< the objective function value of the subproblem */
4113  SCIP_RESULT* result /**< the result from solving the subproblem */
4114  )
4115 {
4116  assert(benders != NULL);
4117  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
4118  assert(benders->benderssolvesubconvex != NULL || benders->benderssolvesub != NULL);
4119 
4120  assert(solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP);
4121 
4122  (*objective) = -SCIPsetInfinity(set);
4123 
4124  /* calls the user defined subproblem solving method. Only the convex relaxations are solved during the Large
4125  * Neighbourhood Benders' Search. */
4126  if( solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX )
4127  {
4128  if( benders->benderssolvesubconvex != NULL )
4129  {
4130  SCIP_CALL( benders->benderssolvesubconvex(set->scip, benders, sol, probnumber,
4131  SCIPbendersOnlyCheckConvexRelax(benders, SCIPsetGetSubscipsOff(set)), objective, result) );
4132  }
4133  else
4134  (*result) = SCIP_DIDNOTRUN;
4135  }
4136  else if( solveloop == SCIP_BENDERSSOLVELOOP_USERCIP )
4137  {
4138  if( benders->benderssolvesub != NULL )
4139  {
4140  SCIP_CALL( benders->benderssolvesub(set->scip, benders, sol, probnumber, objective, result) );
4141  }
4142  else
4143  (*result) = SCIP_DIDNOTRUN;
4144  }
4145 
4146  /* evaluate result */
4147  if( (*result) != SCIP_DIDNOTRUN
4148  && (*result) != SCIP_FEASIBLE
4149  && (*result) != SCIP_INFEASIBLE
4150  && (*result) != SCIP_UNBOUNDED )
4151  {
4152  SCIPerrorMessage("the user-defined solving method for the Benders' decomposition <%s> returned invalid result <%d>\n",
4153  benders->name, *result);
4154  return SCIP_INVALIDRESULT;
4155  }
4156 
4157  if( (*result) == SCIP_INFEASIBLE )
4158  (*infeasible) = TRUE;
4159 
4160  if( (*result) == SCIP_FEASIBLE
4161  && (SCIPsetIsInfinity(set, -(*objective)) || SCIPsetIsInfinity(set, (*objective))) )
4162  {
4163  SCIPerrorMessage("the user-defined solving method for the Benders' decomposition <%s> returned objective value %g\n",
4164  benders->name, (*objective));
4165  return SCIP_ERROR;
4166  }
4167 
4168  /* if the result is SCIP_DIDNOTFIND, then an error is returned and SCIP will terminate. */
4169  if( (*result) == SCIP_DIDNOTFIND )
4170  return SCIP_ERROR;
4171  else
4172  return SCIP_OKAY;
4173 }
4174 
4175 /** executes the subproblem solving process */
4177  SCIP_BENDERS* benders, /**< Benders' decomposition */
4178  SCIP_SET* set, /**< global SCIP settings */
4179  SCIP_SOL* sol, /**< primal CIP solution */
4180  int probnumber, /**< the subproblem number */
4181  SCIP_BENDERSSOLVELOOP solveloop, /**< the solve loop iteration. The first iter is for LP, the second for IP */
4182  SCIP_Bool enhancement, /**< is the solve performed as part of and enhancement? */
4183  SCIP_Bool* solved, /**< flag to indicate whether the subproblem was solved */
4184  SCIP_Bool* infeasible, /**< returns whether the current subproblem is infeasible */
4185  SCIP_BENDERSENFOTYPE type /**< the enforcement type calling this function */
4186  )
4187 { /*lint --e{715}*/
4188  SCIP* subproblem;
4189  SCIP_RESULT result;
4190  SCIP_Real objective;
4191  SCIP_STATUS solvestatus = SCIP_STATUS_UNKNOWN;
4192 
4193  assert(benders != NULL);
4194  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
4195 
4196  SCIPsetDebugMsg(set, "Benders' decomposition: solving subproblem %d\n", probnumber);
4197 
4198  result = SCIP_DIDNOTRUN;
4199  objective = SCIPsetInfinity(set);
4200 
4201  subproblem = SCIPbendersSubproblem(benders, probnumber);
4202 
4203  if( subproblem == NULL && (benders->benderssolvesubconvex == NULL || benders->benderssolvesub == NULL) )
4204  {
4205  SCIPerrorMessage("The subproblem %d is set to NULL, but both bendersSolvesubconvex%s and bendersSolvesub%s "
4206  "are not defined.\n", probnumber, benders->name, benders->name);
4207  SCIPABORT();
4208  return SCIP_ERROR;
4209  }
4210 
4211  /* initially setting the solved flag to FALSE */
4212  (*solved) = FALSE;
4213 
4214  /* if the subproblem solve callback is implemented, then that is used instead of the default setup */
4215  if( solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP )
4216  {
4217  /* calls the user defined subproblem solving method. Only the convex relaxations are solved during the Large
4218  * Neighbourhood Benders' Search. */
4219  SCIP_CALL( executeUserDefinedSolvesub(benders, set, sol, probnumber, solveloop, infeasible, &objective, &result) );
4220 
4221  /* if the result is DIDNOTRUN, then the subproblem was not solved */
4222  (*solved) = (result != SCIP_DIDNOTRUN);
4223  }
4224  else if( subproblem != NULL )
4225  {
4226  /* setting up the subproblem */
4227  if( solveloop == SCIP_BENDERSSOLVELOOP_CONVEX )
4228  {
4229  SCIP_CALL( SCIPbendersSetupSubproblem(benders, set, sol, probnumber, type) );
4230 
4231  /* if the limits of the master problem were hit during the setup process, then the subproblem will not have
4232  * been setup. In this case, the solving function must be exited.
4233  */
4234  if( !SCIPbendersSubproblemIsSetup(benders, probnumber) )
4235  {
4236  SCIPbendersSetSubproblemObjval(benders, probnumber, SCIPsetInfinity(set));
4237  (*solved) = FALSE;
4238  return SCIP_OKAY;
4239  }
4240  }
4241  else
4242  {
4243  SCIP_CALL( updateEventhdlrUpperbound(benders, probnumber, SCIPbendersGetAuxiliaryVarVal(benders, set, sol, probnumber)) );
4244  }
4245 
4246  /* solving the subproblem
4247  * the LP of the subproblem is solved in the first solveloop.
4248  * In the second solve loop, the MIP problem is solved */
4249  if( solveloop == SCIP_BENDERSSOLVELOOP_CONVEX
4251  {
4252  SCIP_CALL( SCIPbendersSolveSubproblemLP(set->scip, benders, probnumber, &solvestatus, &objective) );
4253 
4254  /* if the (N)LP was solved without error, then the subproblem is labelled as solved */
4255  if( solvestatus == SCIP_STATUS_OPTIMAL || solvestatus == SCIP_STATUS_INFEASIBLE )
4256  (*solved) = TRUE;
4257 
4258  if( solvestatus == SCIP_STATUS_INFEASIBLE )
4259  (*infeasible) = TRUE;
4260  }
4261  else
4262  {
4263  SCIP_SOL* bestsol;
4264 
4265  SCIP_CALL( SCIPbendersSolveSubproblemCIP(set->scip, benders, probnumber, &solvestatus, FALSE) );
4266 
4267  if( solvestatus == SCIP_STATUS_INFEASIBLE )
4268  (*infeasible) = TRUE;
4269 
4270  /* if the generic subproblem solving methods are used, then the CIP subproblems are always solved. */
4271  (*solved) = TRUE;
4272 
4273  bestsol = SCIPgetBestSol(subproblem);
4274  if( bestsol != NULL )
4275  objective = SCIPgetSolOrigObj(subproblem, bestsol)*(int)SCIPgetObjsense(set->scip);
4276  else
4277  objective = SCIPsetInfinity(set);
4278  }
4279  }
4280  else
4281  {
4282  SCIPABORT();
4283  }
4284 
4285  if( !enhancement )
4286  {
4287  /* The following handles the cases when the subproblem is OPTIMAL, INFEASIBLE and UNBOUNDED.
4288  * If a subproblem is unbounded, then the auxiliary variables are set to -infinity and the unbounded flag is
4289  * returned as TRUE. No cut will be generated, but the result will be set to SCIP_FEASIBLE.
4290  */
4291  if( solveloop == SCIP_BENDERSSOLVELOOP_CONVEX || solveloop == SCIP_BENDERSSOLVELOOP_CIP )
4292  {
4293  /* TODO: Consider whether other solutions status should be handled */
4294  if( solvestatus == SCIP_STATUS_OPTIMAL )
4295  SCIPbendersSetSubproblemObjval(benders, probnumber, objective);
4296  else if( solvestatus == SCIP_STATUS_INFEASIBLE )
4297  SCIPbendersSetSubproblemObjval(benders, probnumber, SCIPsetInfinity(set));
4298  else if( solvestatus == SCIP_STATUS_USERINTERRUPT || solvestatus == SCIP_STATUS_BESTSOLLIMIT )
4299  SCIPbendersSetSubproblemObjval(benders, probnumber, objective);
4300  else if( solvestatus == SCIP_STATUS_MEMLIMIT || solvestatus == SCIP_STATUS_TIMELIMIT
4301  || solvestatus == SCIP_STATUS_UNKNOWN )
4302  {
4303  SCIPverbMessage(set->scip, SCIP_VERBLEVEL_FULL, NULL, " Benders' decomposition: Error solving "
4304  "subproblem %d. No cut will be generated for this subproblem.\n", probnumber);
4305  SCIPbendersSetSubproblemObjval(benders, probnumber, SCIPsetInfinity(set));
4306  }
4307  else if( solvestatus == SCIP_STATUS_UNBOUNDED )
4308  {
4309  SCIPerrorMessage("The Benders' decomposition subproblem %d is unbounded. This should not happen.\n",
4310  probnumber);
4311  SCIPABORT();
4312  }
4313  else
4314  {
4315  SCIPerrorMessage("Invalid status returned from solving Benders' decomposition subproblem %d. Solution status: %d\n",
4316  probnumber, solvestatus);
4317  SCIPABORT();
4318  }
4319  }
4320  else
4321  {
4322  assert(solveloop == SCIP_BENDERSSOLVELOOP_USERCONVEX || solveloop == SCIP_BENDERSSOLVELOOP_USERCIP);
4323  if( result == SCIP_FEASIBLE )
4324  SCIPbendersSetSubproblemObjval(benders, probnumber, objective);
4325  else if( result == SCIP_INFEASIBLE )
4326  SCIPbendersSetSubproblemObjval(benders, probnumber, SCIPsetInfinity(set));
4327  else if( result == SCIP_UNBOUNDED )
4328  {
4329  SCIPerrorMessage("The Benders' decomposition subproblem %d is unbounded. This should not happen.\n",
4330  probnumber);
4331  SCIPABORT();
4332  }
4333  else if( result != SCIP_DIDNOTRUN )
4334  {
4335  SCIPerrorMessage("Invalid result <%d> from user-defined subproblem solving method. This should not happen.\n",
4336  result);
4337  }
4338  }
4339  }
4340 
4341  return SCIP_OKAY;
4342 }
4343 
4344 /** sets up the subproblem using the solution to the master problem */
4346  SCIP_BENDERS* benders, /**< Benders' decomposition */
4347  SCIP_SET* set, /**< global SCIP settings */
4348  SCIP_SOL* sol, /**< primal CIP solution */
4349  int probnumber, /**< the subproblem number */
4350  SCIP_BENDERSENFOTYPE type /**< the enforcement type calling this function */
4351  )
4352 {
4353  SCIP* subproblem;
4354  SCIP_VAR** vars;
4355  SCIP_VAR* mastervar;
4356  SCIP_Real solval;
4357  int nvars;
4358  int i;
4359 
4360  assert(benders != NULL);
4361  assert(set != NULL);
4362  assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
4363 
4364  subproblem = SCIPbendersSubproblem(benders, probnumber);
4365 
4366  /* the subproblem setup can only be performed if the subproblem is not NULL */
4367  if( subproblem == NULL )
4368  {
4369  SCIPerrorMessage("The subproblem %d is NULL. Thus, the subproblem setup must be performed manually in either "
4370  "bendersSolvesubconvex%s or bendersSolvesub%s.\n", probnumber, benders->name, benders->name);
4371  return SCIP_ERROR;
4372  }
4373  assert(subproblem != NULL);
4374 
4375  /* changing all of the master problem variable to continuous. */
4376  SCIP_CALL( SCIPbendersChgMastervarsToCont(benders, set, probnumber) );
4377 
4378  /* if the Benders' decomposition subproblem is convex and has continuous variables, then probing mode
4379  * must be started.
4380  * If the subproblem contains non-convex constraints or discrete variables, then the problem must be initialised,
4381  * and then put into SCIP_STAGE_SOLVING to be able to change the variable bounds. The probing mode is entered once
4382  * the variable bounds are set.
4383  * In the latter case, the transformed problem is freed after each subproblem solve round. */
4384  if( SCIPbendersGetSubproblemType(benders, probnumber) == SCIP_BENDERSSUBTYPE_CONVEXCONT )
4385  {
4386  SCIP_CALL( SCIPstartProbing(subproblem) );
4387  }
4388  else
4389  {
4390  SCIP_Bool success;
4391 
4392  SCIP_CALL( initialiseSubproblem(benders, set, probnumber, &success) );
4393 
4394  if( !success )
4395  {
4396  /* set the flag to indicate that the subproblems have been set up */
4397  SCIPbendersSetSubproblemIsSetup(benders, probnumber, FALSE);
4398 
4399  return SCIP_OKAY;
4400  }
4401  }
4402 
4403  vars = SCIPgetVars(subproblem);
4404  nvars = SCIPgetNVars(subproblem);
4405 
4406  /* looping over all variables in the subproblem to find those corresponding to the master problem variables. */
4407  /* TODO: It should be possible to store the pointers to the master variables to speed up the subproblem setup */
4408  for( i = 0; i < nvars; i++ )
4409  {
4410  SCIP_CALL( SCIPbendersGetVar(benders, set, vars[i], &mastervar, -1) );
4411 
4412  if( mastervar != NULL )
4413  {
4414  /* It is possible due to numerics that the solution value exceeds the upper or lower bounds. When this
4415  * happens, it causes an error in the LP solver as a result of inconsistent bounds. So the following statements
4416  * are used to ensure that the bounds are not exceeded when applying the fixings for the Benders'
4417  * decomposition subproblems
4418  */
4419  solval = SCIPgetSolVal(set->scip, sol, mastervar);
4420  if( !SCIPisLT(set->scip, solval, SCIPvarGetUbLocal(vars[i])) )
4421  solval = SCIPvarGetUbLocal(vars[i]);
4422  else if( !SCIPisGT(set->scip, solval, SCIPvarGetLbLocal(vars[i])) )
4423  solval = SCIPvarGetLbLocal(vars[i]);
4424 
4425  /* fixing the variable in the subproblem */
4426  if( !SCIPisEQ(subproblem, SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])) )
4427  {
4428  if( SCIPisGT(subproblem, solval, SCIPvarGetLbLocal(vars[i])) )
4429  {
4430  SCIP_CALL( SCIPchgVarLb(subproblem, vars[i], solval) );
4431  }
4432  if( SCIPisLT(subproblem, solval, SCIPvarGetUbLocal(vars[i])) )
4433  {
4434  SCIP_CALL( SCIPchgVarUb(subproblem, vars[i], solval) );
4435  }
4436  }
4437 
4438  assert(SCIPisEQ(subproblem, SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])));
4439  }
4440  else if( strstr(SCIPvarGetName(vars[i]), SLACKVAR_NAME) != NULL )
4441  {
4442  /* if the slack variables have been added to help improve feasibility, then they remain unfixed with a large
4443  * objective coefficient. Once the root node has been solved to optimality, then the slack variables are
4444  * fixed to zero.
4445  */
4446  if( benders->feasibilityphase && SCIPgetDepth(set->scip) == 0 && type != SCIP_BENDERSENFOTYPE_CHECK )
4447  {
4448  /* The coefficient update or variable fixing can only be performed if the subproblem is in probing mode.
4449  * If the slack var coef gets very large, then we fix the slack variable to 0 instead.
4450  */
4451  if( SCIPinProbing(subproblem) )
4452  {
4453  if( benders->slackvarcoef <= benders->maxslackvarcoef )
4454  {
4455  SCIP_CALL( SCIPchgVarObjProbing(subproblem, vars[i], benders->slackvarcoef) );
4456  }
4457  else
4458  {
4459  SCIP_CALL( SCIPchgVarUbProbing(subproblem, vars[i], 0.0) );
4460  }
4461  }
4462  }
4463  else
4464  {
4465  /* if the subproblem is non-linear and convex, then slack variables have been added to the subproblem. These
4466  * need to be fixed to zero when first solving the subproblem. However, if the slack variables have been added
4467  * by setting the execfeasphase runtime parameter, then they must not get fixed to zero
4468  */
4469  assert( !SCIPisEQ(subproblem, SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])) );
4470  assert( SCIPisZero(subproblem, SCIPvarGetLbLocal(vars[i])) );
4471 
4472  if( SCIPisLT(subproblem, 0.0, SCIPvarGetUbLocal(vars[i])) )
4473  {
4474  SCIP_CALL( SCIPchgVarUb(subproblem, vars[i], 0.0) );
4475  }
4476  }
4477  }
4478  }
4479 
4480  /* if the subproblem contain non-convex constraints or discrete variables, then the probing mode is entered after
4481  * setting up the subproblem
4482  */
4483  if( SCIPbendersGetSubproblemType(benders, probnumber) != SCIP_BENDERSSUBTYPE_CONVEXCONT )
4484  {
4485  SCIP_CALL( SCIPstartProbing(subproblem) );
4486  }
4487 
4488  /* set the flag to indicate that the subproblems have been set up */
4489  SCIPbendersSetSubproblemIsSetup(benders, probnumber, TRUE);
4490 
4491  return SCIP_OKAY;
4492 }
4493 
4494 /** Solve a Benders' decomposition subproblems. This will either call the user defined method or the generic solving
4495  * methods. If the generic method is called, then the subproblem must be set up before calling this method. */
4497  SCIP_BENDERS* benders, /**< Benders' decomposition */
4498  SCIP_SET* set, /**< global SCIP settings */
4499  SCIP_SOL* sol, /**< primal CIP solution, can be NULL */
4500  int probnumber, /**< the subproblem number */
4501  SCIP_Bool* infeasible, /**< returns whether the current subproblem is infeasible */
4502  SCIP_Bool solvecip, /**< directly solve the CIP subproblem */
4503  SCIP_Real* objective /**< the objective function value of the subproblem, can be NULL */
4504  )
4505 {
4506  assert(benders != NULL);
4507  assert(set != NULL);
4508  assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
4509 
4510  (*infeasible) = FALSE;
4511 
4512  /* the subproblem must be set up before this function is called. */
4513  if( SCIPbendersSubproblem(benders, probnumber) != NULL && !SCIPbendersSubproblemIsSetup(benders, probnumber)
4514  && !SCIPbendersSubproblemIsIndependent(benders, probnumber) )
4515  {
4516  SCIPerrorMessage("Benders' decomposition subproblem %d must be set up before calling SCIPbendersSolveSubproblem(). Call SCIPsetupSubproblem() first.\n", probnumber);
4517  return SCIP_ERROR;
4518  }
4519 
4520  /* if the subproblem solve callback is implemented, then that is used instead of the default setup */
4521  if( benders->benderssolvesubconvex != NULL || benders->benderssolvesub != NULL)
4522  {
4523  SCIP_BENDERSSOLVELOOP solveloop;
4524  SCIP_RESULT result;
4525  SCIP_Real subobj;
4526 
4527  if( solvecip )
4528  solveloop = SCIP_BENDERSSOLVELOOP_USERCIP;
4529  else
4531 
4532  SCIP_CALL( executeUserDefinedSolvesub(benders, set, sol, probnumber, solveloop, infeasible, &subobj, &result) );
4533 
4534  if( objective != NULL )
4535  (*objective) = subobj;
4536  }
4537  else
4538  {
4539  SCIP* subproblem;
4540 
4541  subproblem = SCIPbendersSubproblem(benders, probnumber);
4542  assert(subproblem != NULL);
4543 
4544  /* solving the subproblem */
4545  if( solvecip && SCIPbendersGetSubproblemType(benders, probnumber) != SCIP_BENDERSSUBTYPE_CONVEXCONT )
4546  {
4547  SCIP_STATUS solvestatus;
4548 
4549  SCIP_CALL( SCIPbendersSolveSubproblemCIP(set->scip, benders, probnumber, &solvestatus, solvecip) );
4550 
4551  if( solvestatus == SCIP_STATUS_INFEASIBLE )
4552  (*infeasible) = TRUE;
4553  if( objective != NULL )
4554  (*objective) = SCIPgetSolOrigObj(subproblem, SCIPgetBestSol(subproblem))*(int)SCIPgetObjsense(subproblem);
4555  }
4556  else
4557  {
4558  SCIP_Bool success;
4559 
4560  /* if the subproblem has convex constraints and continuous variables, then it should have been initialised and
4561  * in SCIP_STAGE_SOLVING. In this case, the subproblem only needs to be put into probing mode.
4562  */
4563  if( SCIPbendersGetSubproblemType(benders, probnumber) == SCIP_BENDERSSUBTYPE_CONVEXCONT )
4564  {
4565  /* if the subproblem is not in probing mode, then it must be put into that mode for the LP solve. */
4566  if( !SCIPinProbing(subproblem) )
4567  {
4568  SCIP_CALL( SCIPstartProbing(subproblem) );
4569  }
4570 
4571  success = TRUE;
4572  }
4573  else
4574  {
4575  SCIP_CALL( initialiseSubproblem(benders, set, probnumber, &success) );
4576  }
4577 
4578  /* if setting up the subproblem was successful */
4579  if( success )
4580  {
4581  SCIP_STATUS solvestatus;
4582  SCIP_Real lpobjective;
4583 
4584  SCIP_CALL( SCIPbendersSolveSubproblemLP(set->scip, benders, probnumber, &solvestatus, &lpobjective) );
4585 
4586  if( solvestatus == SCIP_STATUS_INFEASIBLE )
4587  (*infeasible) = TRUE;
4588  else if( objective != NULL )
4589  (*objective) = lpobjective;
4590  }
4591  else
4592  {
4593  if( objective != NULL )
4594  (*objective) = SCIPinfinity(subproblem);
4595  }
4596  }
4597  }
4598 
4599  return SCIP_OKAY;
4600 }
4601 
4602 /** copies the time and memory limit from the master problem to the subproblem */
4603 static
4605  SCIP* scip, /**< the SCIP data structure */
4606  SCIP* subproblem /**< the Benders' decomposition subproblem */
4607  )
4608 {
4609  SCIP_Real mastertimelimit;
4610  SCIP_Real subtimelimit;
4611  SCIP_Real maxsubtimelimit;
4612  SCIP_Real mastermemorylimit;
4613  SCIP_Real submemorylimit;
4614  SCIP_Real maxsubmemorylimit;
4615 
4616  assert(scip != NULL);
4617 
4618  /* setting the time limit for the Benders' decomposition subproblems. It is set to 102% of the remaining time. */
4619  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &mastertimelimit) );
4620  maxsubtimelimit = SCIPparamGetRealMax(SCIPgetParam(subproblem, "limits/time"));
4621  subtimelimit = (mastertimelimit - SCIPgetSolvingTime(scip)) * 1.02;
4622  subtimelimit = MIN(subtimelimit, maxsubtimelimit);
4623  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/time", MAX(0.0, subtimelimit)) );
4624 
4625  /* setting the memory limit for the Benders' decomposition subproblems. */
4626  SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &mastermemorylimit) );
4627  maxsubmemorylimit = SCIPparamGetRealMax(SCIPgetParam(subproblem, "limits/memory"));
4628  submemorylimit = mastermemorylimit - (SCIPgetMemUsed(scip) + SCIPgetMemExternEstim(scip))/1048576.0;
4629  submemorylimit = MIN(submemorylimit, maxsubmemorylimit);
4630  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/memory", MAX(0.0, submemorylimit)) );
4631 
4632  return SCIP_OKAY;
4633 }
4634 
4635 /** stores the original parameters from the subproblem */
4636 static
4638  SCIP* subproblem, /**< the SCIP data structure */
4639  SCIP_SUBPROBPARAMS* origparams /**< the original subproblem parameters */
4640  )
4641 {
4642  assert(subproblem != NULL);
4643  assert(origparams != NULL);
4644 
4645  SCIP_CALL( SCIPgetRealParam(subproblem, "limits/memory", &origparams->limits_memory) );
4646  SCIP_CALL( SCIPgetRealParam(subproblem, "limits/time", &origparams->limits_time) );
4647  SCIP_CALL( SCIPgetBoolParam(subproblem, "conflict/enable", &origparams->conflict_enable) );
4648  SCIP_CALL( SCIPgetIntParam(subproblem, "lp/disablecutoff", &origparams->lp_disablecutoff) );
4649  SCIP_CALL( SCIPgetIntParam(subproblem, "lp/scaling", &origparams->lp_scaling) );
4650  SCIP_CALL( SCIPgetCharParam(subproblem, "lp/initalgorithm", &origparams->lp_initalg) );
4651  SCIP_CALL( SCIPgetCharParam(subproblem, "lp/resolvealgorithm", &origparams->lp_resolvealg) );
4652  SCIP_CALL( SCIPgetBoolParam(subproblem, "lp/alwaysgetduals", &origparams->lp_alwaysgetduals) );
4653  SCIP_CALL( SCIPgetBoolParam(subproblem, "misc/scaleobj", &origparams->misc_scaleobj) );
4654  SCIP_CALL( SCIPgetBoolParam(subproblem, "misc/catchctrlc", &origparams->misc_catchctrlc) );
4655  SCIP_CALL( SCIPgetIntParam(subproblem, "propagating/maxrounds", &origparams->prop_maxrounds) );
4656  SCIP_CALL( SCIPgetIntParam(subproblem, "propagating/maxroundsroot", &origparams->prop_maxroundsroot) );
4657  SCIP_CALL( SCIPgetIntParam(subproblem, "constraints/linear/propfreq", &origparams->cons_linear_propfreq) );
4658 
4659  return SCIP_OKAY;
4660 }
4661 
4662 /** sets the parameters for the subproblem */
4663 static
4665  SCIP* scip, /**< the SCIP data structure */
4666  SCIP* subproblem /**< the subproblem SCIP instance */
4667  )
4668 {
4669  assert(scip != NULL);
4670  assert(subproblem != NULL);
4671 
4672  /* copying memory and time limits */
4673  SCIP_CALL( copyMemoryAndTimeLimits(scip, subproblem) );
4674 
4675  /* Do we have to disable presolving? If yes, we have to store all presolving parameters. */
4677 
4678  /* Disabling heuristics so that the problem is not trivially solved */
4680 
4681  /* store parameters that are changed for the generation of the subproblem cuts */
4682  SCIP_CALL( SCIPsetBoolParam(subproblem, "conflict/enable", FALSE) );
4683 
4684  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/disablecutoff", 1) );
4685  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/scaling", 0) );
4686 
4687  SCIP_CALL( SCIPsetCharParam(subproblem, "lp/initalgorithm", 'd') );
4688  SCIP_CALL( SCIPsetCharParam(subproblem, "lp/resolvealgorithm", 'd') );
4689 
4690  SCIP_CALL( SCIPsetBoolParam(subproblem, "lp/alwaysgetduals", TRUE) );
4691  SCIP_CALL( SCIPsetBoolParam(subproblem, "misc/scaleobj", FALSE) );
4692 
4693  /* do not abort subproblem on CTRL-C */
4694  SCIP_CALL( SCIPsetBoolParam(subproblem, "misc/catchctrlc", FALSE) );
4695 
4696 #ifndef SCIP_MOREDEBUG
4697  SCIP_CALL( SCIPsetIntParam(subproblem, "display/verblevel", (int)SCIP_VERBLEVEL_NONE) );
4698 #endif
4699 
4700  SCIP_CALL( SCIPsetIntParam(subproblem, "propagating/maxrounds", 0) );
4701  SCIP_CALL( SCIPsetIntParam(subproblem, "propagating/maxroundsroot", 0) );
4702 
4703  SCIP_CALL( SCIPsetIntParam(subproblem, "constraints/linear/propfreq", -1) );
4704 
4705  SCIP_CALL( SCIPsetIntParam(subproblem, "heuristics/alns/freq", -1) );
4706 
4707  SCIP_CALL( SCIPsetIntParam(subproblem, "separating/aggregation/freq", -1) );
4708  SCIP_CALL( SCIPsetIntParam(subproblem, "separating/gomory/freq", -1) );
4709 
4710  return SCIP_OKAY;
4711 }
4712 
4713 /** resets the original parameters from the subproblem */
4714 static
4716  SCIP* subproblem, /**< the SCIP data structure */
4717  SCIP_SUBPROBPARAMS* origparams /**< the original subproblem parameters */
4718  )
4719 {
4720  assert(subproblem != NULL);
4721  assert(origparams != NULL);
4722 
4723  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/memory", origparams->limits_memory) );
4724  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/time", origparams->limits_time) );
4725  SCIP_CALL( SCIPsetBoolParam(subproblem, "conflict/enable", origparams->conflict_enable) );
4726  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/disablecutoff", origparams->lp_disablecutoff) );
4727  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/scaling", origparams->lp_scaling) );
4728  SCIP_CALL( SCIPsetCharParam(subproblem, "lp/initalgorithm", origparams->lp_initalg) );
4729  SCIP_CALL( SCIPsetCharParam(subproblem, "lp/resolvealgorithm", origparams->lp_resolvealg) );
4730  SCIP_CALL( SCIPsetBoolParam(subproblem, "lp/alwaysgetduals", origparams->lp_alwaysgetduals) );
4731  SCIP_CALL( SCIPsetBoolParam(subproblem, "misc/scaleobj", origparams->misc_scaleobj) );
4732  SCIP_CALL( SCIPsetBoolParam(subproblem, "misc/catchctrlc", origparams->misc_catchctrlc) );
4733  SCIP_CALL( SCIPsetIntParam(subproblem, "propagating/maxrounds", origparams->prop_maxrounds) );
4734  SCIP_CALL( SCIPsetIntParam(subproblem, "propagating/maxroundsroot", origparams->prop_maxroundsroot) );
4735  SCIP_CALL( SCIPsetIntParam(subproblem, "constraints/linear/propfreq", origparams->cons_linear_propfreq) );
4736 
4737  return SCIP_OKAY;
4738 }
4739 
4740 /** solves the LP of the Benders' decomposition subproblem
4741  *
4742  * This requires that the subproblem is in probing mode.
4743  */
4745  SCIP* scip, /**< the SCIP data structure */
4746  SCIP_BENDERS* benders, /**< the Benders' decomposition data structure */
4747  int probnumber, /**< the subproblem number */
4748  SCIP_STATUS* solvestatus, /**< status of subproblem solve */
4749  SCIP_Real* objective /**< optimal value of subproblem, if solved to optimality */
4750  )
4751 {
4752  SCIP* subproblem;
4753  SCIP_SUBPROBPARAMS* origparams;
4754  SCIP_Bool solvenlp;
4755 
4756  assert(benders != NULL);
4757  assert(solvestatus != NULL);
4758  assert(objective != NULL);
4759  assert(SCIPbendersSubproblemIsSetup(benders, probnumber));
4760 
4761  /* TODO: This should be solved just as an LP, so as a MIP. There is too much overhead with the MIP.
4762  * Need to change status check for checking the LP. */
4763  subproblem = SCIPbendersSubproblem(benders, probnumber);
4764  assert(subproblem != NULL);
4765 
4766  /* only solve the NLP relaxation if the NLP has been constructed and there exists an NLPI. If it is not possible to
4767  * solve the NLP relaxation, then the LP relaxation is used to generate Benders' cuts
4768  */
4769  solvenlp = FALSE;
4770  if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) > 0
4771  && SCIPbendersGetSubproblemType(benders, probnumber) <= SCIP_BENDERSSUBTYPE_CONVEXDIS )
4772  solvenlp = TRUE;
4773 
4774  *objective = SCIPinfinity(subproblem);
4775 
4776  assert(SCIPisNLPConstructed(subproblem) || SCIPisLPConstructed(subproblem));
4777  assert(SCIPinProbing(subproblem));
4778 
4779  /* allocating memory for the parameter storage */
4780  SCIP_CALL( SCIPallocBlockMemory(subproblem, &origparams) );
4781 
4782  /* store the original parameters of the subproblem */
4783  SCIP_CALL( storeOrigSubproblemParams(subproblem, origparams) );
4784 
4785  /* setting the subproblem parameters */
4786  SCIP_CALL( setSubproblemParams(scip, subproblem) );
4787 
4788  if( solvenlp )
4789  {
4790  SCIP_NLPSOLSTAT nlpsolstat;
4791  SCIP_NLPTERMSTAT nlptermstat;
4792 #ifdef SCIP_MOREDEBUG
4793  SCIP_SOL* nlpsol;
4794 #endif
4795 
4796  SCIP_CALL( SCIPsolveNLP(subproblem) ); /*lint !e666*/
4797 
4798  nlpsolstat = SCIPgetNLPSolstat(subproblem);
4799  nlptermstat = SCIPgetNLPTermstat(subproblem);
4800  SCIPdebugMsg(scip, "NLP solstat %d termstat %d\n", nlpsolstat, nlptermstat);
4801 
4802  if( nlptermstat == SCIP_NLPTERMSTAT_OKAY && (nlpsolstat == SCIP_NLPSOLSTAT_LOCINFEASIBLE || nlpsolstat == SCIP_NLPSOLSTAT_GLOBINFEASIBLE) )
4803  {
4804  /* trust infeasible only if terminated "okay" */
4805  (*solvestatus) = SCIP_STATUS_INFEASIBLE;
4806  }
4807  else if( nlpsolstat == SCIP_NLPSOLSTAT_LOCOPT || nlpsolstat == SCIP_NLPSOLSTAT_GLOBOPT
4808  || nlpsolstat == SCIP_NLPSOLSTAT_FEASIBLE )
4809  {
4810 #ifdef SCIP_MOREDEBUG
4811  SCIP_CALL( SCIPcreateNLPSol(subproblem, &nlpsol, NULL) );
4812  SCIP_CALL( SCIPprintSol(subproblem, nlpsol, NULL, FALSE) );
4813  SCIP_CALL( SCIPfreeSol(subproblem, &nlpsol) );
4814 #endif
4815 
4816  (*solvestatus) = SCIP_STATUS_OPTIMAL;
4817  (*objective) = SCIPretransformObj(subproblem, SCIPgetNLPObjval(subproblem));
4818  }
4819  else if( nlpsolstat == SCIP_NLPSOLSTAT_UNBOUNDED )
4820  {
4821  (*solvestatus) = SCIP_STATUS_UNBOUNDED;
4822  SCIPerrorMessage("The NLP of Benders' decomposition subproblem %d is unbounded. This should not happen.\n",
4823  probnumber);
4824  SCIPABORT();
4825  }
4826  else if( nlptermstat == SCIP_NLPTERMSTAT_TIMELIMIT )
4827  {
4828  (*solvestatus) = SCIP_STATUS_TIMELIMIT;
4829  }
4830  else if( nlptermstat == SCIP_NLPTERMSTAT_INTERRUPT )
4831  {
4832  (*solvestatus) = SCIP_STATUS_USERINTERRUPT;
4833  }
4834  else
4835  {
4836  SCIPerrorMessage("Invalid solution status: %d. Termination status: %d. Solving the NLP relaxation of Benders' decomposition subproblem %d.\n",
4837  nlpsolstat, nlptermstat, probnumber);
4838  SCIPABORT();
4839  }
4840  }
4841  else
4842  {
4843  SCIP_Bool lperror;
4844  SCIP_Bool cutoff;
4845 
4846  SCIP_CALL( SCIPsolveProbingLP(subproblem, -1, &lperror, &cutoff) );
4847 
4848  switch( SCIPgetLPSolstat(subproblem) )
4849  {
4851  {
4852  (*solvestatus) = SCIP_STATUS_INFEASIBLE;
4853  break;
4854  }
4855 
4856  case SCIP_LPSOLSTAT_OPTIMAL :
4857  {
4858  (*solvestatus) = SCIP_STATUS_OPTIMAL;
4859  (*objective) = SCIPgetSolOrigObj(subproblem, NULL)*(int)SCIPgetObjsense(scip);
4860  break;
4861  }
4862 
4864  {
4865  (*solvestatus) = SCIP_STATUS_UNBOUNDED;
4866  SCIPerrorMessage("The LP of Benders' decomposition subproblem %d is unbounded. This should not happen.\n",
4867  probnumber);
4868  SCIPABORT();
4869  break;
4870  }
4871 
4872  case SCIP_LPSOLSTAT_ERROR :
4875  {
4876  if( SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_TIMELIMIT )
4877  (*solvestatus) = SCIP_STATUS_TIMELIMIT;
4878  else
4879  (*solvestatus) = SCIP_STATUS_UNKNOWN;
4880 
4881  SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, " Benders' decomposition: Error solving LP "
4882  "relaxation of subproblem %d. No cut will be generated for this subproblem.\n", probnumber);
4883  break;
4884  }
4885 
4888  default:
4889  {
4890  SCIPerrorMessage("Invalid status: %d. Solving the LP relaxation of Benders' decomposition subproblem %d.\n",
4891  SCIPgetLPSolstat(subproblem), probnumber);
4892  SCIPABORT();
4893  break;
4894  }
4895  }
4896  }
4897 
4898  /* resetting the subproblem parameters */
4899  SCIP_CALL( resetOrigSubproblemParams(subproblem, origparams) );
4900 
4901  /* freeing the parameter storage */
4902  SCIPfreeBlockMemory(subproblem, &origparams);
4903 
4904  return SCIP_OKAY;
4905 }
4906 
4907 /** solves the Benders' decomposition subproblem */
4909  SCIP* scip, /**< the SCIP data structure */
4910  SCIP_BENDERS* benders, /**< the Benders' decomposition data structure */
4911  int probnumber, /**< the subproblem number */
4912  SCIP_STATUS* solvestatus, /**< status of subproblem solve */
4913  SCIP_Bool solvecip /**< directly solve the CIP subproblem */
4914  )
4915 {
4916  SCIP* subproblem;
4917  SCIP_SUBPROBPARAMS* origparams;
4918 
4919  assert(benders != NULL);
4920  assert(solvestatus != NULL);
4921 
4922  subproblem = SCIPbendersSubproblem(benders, probnumber);
4923  assert(subproblem != NULL);
4924 
4925  /* allocating memory for the parameter storage */
4926  SCIP_CALL( SCIPallocBlockMemory(subproblem, &origparams) );
4927 
4928  /* store the original parameters of the subproblem */
4929  SCIP_CALL( storeOrigSubproblemParams(subproblem, origparams) );
4930 
4931  /* If the solve has been stopped for the subproblem, then we need to restart it to complete the solve. The subproblem
4932  * is stopped when it is a MIP so that LP cuts and IP cuts can be generated. */
4933  if( SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING )
4934  {
4935  /* the subproblem should be in probing mode. Otherwise, the event handler did not work correctly */
4936  assert( SCIPinProbing(subproblem) );
4937 
4938  /* the probing mode needs to be stopped so that the MIP can be solved */
4939  SCIP_CALL( SCIPendProbing(subproblem) );
4940 
4941  /* the problem was interrupted in the event handler, so SCIP needs to be informed that the problem is to be restarted */
4942  SCIP_CALL( SCIPrestartSolve(subproblem) );
4943  }
4944  else if( solvecip )
4945  {
4946  /* if the MIP will be solved directly, then the probing mode needs to be skipped.
4947  * This is achieved by setting the solvecip flag in the event handler data to TRUE
4948  */
4949  SCIP_EVENTHDLR* eventhdlr;
4950  SCIP_EVENTHDLRDATA* eventhdlrdata;
4951 
4952  eventhdlr = SCIPfindEventhdlr(subproblem, MIPNODEFOCUS_EVENTHDLR_NAME);
4953  eventhdlrdata = SCIPeventhdlrGetData(eventhdlr);
4954 
4955  eventhdlrdata->solvecip = TRUE;
4956  }
4957  else
4958  {
4959  /* if the problem is not in probing mode, then we need to solve the LP. That requires all methods that will
4960  * modify the structure of the problem need to be deactivated */
4961 
4962  /* setting the subproblem parameters */
4963  SCIP_CALL( setSubproblemParams(scip, subproblem) );
4964 
4965 #ifdef SCIP_EVENMOREDEBUG
4966  SCIP_CALL( SCIPsetBoolParam(subproblem, "display/lpinfo", TRUE) );
4967 #endif
4968  }
4969 
4970 #ifdef SCIP_MOREDEBUG
4971  SCIP_CALL( SCIPsetIntParam(subproblem, "display/verblevel", (int)SCIP_VERBLEVEL_FULL) );
4972  SCIP_CALL( SCIPsetIntParam(subproblem, "display/freq", 1) );
4973 #endif
4974 
4975  SCIP_CALL( SCIPsolve(subproblem) );
4976 
4977  *solvestatus = SCIPgetStatus(subproblem);
4978 
4979  if( *solvestatus != SCIP_STATUS_OPTIMAL && *solvestatus != SCIP_STATUS_UNBOUNDED
4980  && *solvestatus != SCIP_STATUS_INFEASIBLE && *solvestatus != SCIP_STATUS_USERINTERRUPT
4981  && *solvestatus != SCIP_STATUS_BESTSOLLIMIT && *solvestatus != SCIP_STATUS_TIMELIMIT
4982  && *solvestatus != SCIP_STATUS_MEMLIMIT )
4983  {
4984  SCIPerrorMessage("Invalid status: %d. Solving the CIP of Benders' decomposition subproblem %d.\n",
4985  *solvestatus, probnumber);
4986  SCIPABORT();
4987  }
4988 
4989  /* resetting the subproblem parameters */
4990  SCIP_CALL( resetOrigSubproblemParams(subproblem, origparams) );
4991 
4992  /* freeing the parameter storage */
4993  SCIPfreeBlockMemory(subproblem, &origparams);
4994 
4995  return SCIP_OKAY;
4996 }
4997 
4998 /** frees the subproblems */
5000  SCIP_BENDERS* benders, /**< Benders' decomposition */
5001  SCIP_SET* set, /**< global SCIP settings */
5002  int probnumber /**< the subproblem number */
5003  )
5004 {
5005  assert(benders != NULL);
5006  assert(benders->bendersfreesub != NULL
5007  || (benders->bendersfreesub == NULL && benders->benderssolvesubconvex == NULL && benders->benderssolvesub == NULL));
5008  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
5009 
5010  if( benders->bendersfreesub != NULL )
5011  {
5012  SCIP_CALL( benders->bendersfreesub(set->scip, benders, probnumber) );
5013  }
5014  else
5015  {
5016  /* the subproblem is only freed if it is not independent */
5017  if( subproblemIsActive(benders, probnumber) )
5018  {
5019  SCIP* subproblem = SCIPbendersSubproblem(benders, probnumber);
5020 
5021  if( SCIPbendersGetSubproblemType(benders, probnumber) == SCIP_BENDERSSUBTYPE_CONVEXCONT )
5022  {
5023  /* ending probing mode to reset the current node. The probing mode will be restarted at the next solve */
5024  if( SCIPinProbing(subproblem) )
5025  {
5026  SCIP_CALL( SCIPendProbing(subproblem) );
5027  }
5028  }
5029  else
5030  {
5031  /* if the subproblems were solved as part of an enforcement stage, then they will still be in probing mode. The
5032  * probing mode must first be finished and then the problem can be freed */
5033  if( SCIPgetStage(subproblem) >= SCIP_STAGE_TRANSFORMED && SCIPinProbing(subproblem) )
5034  {
5035  SCIP_CALL( SCIPendProbing(subproblem) );
5036  }
5037 
5038  SCIP_CALL( SCIPfreeTransform(subproblem) );
5039  }
5040  }
5041  }
5042 
5043  /* setting the setup flag for the subproblem to FALSE */
5044  SCIPbendersSetSubproblemIsSetup(benders, probnumber, FALSE);
5045  return SCIP_OKAY;
5046 }
5047 
5048 /** compares the subproblem objective value with the auxiliary variable value for optimality */
5050  SCIP_BENDERS* benders, /**< the benders' decomposition structure */
5051  SCIP_SET* set, /**< global SCIP settings */
5052  SCIP_SOL* sol, /**< primal CIP solution */
5053  int probnumber /**< the subproblem number */
5054  )
5055 {
5056  SCIP_Real auxiliaryvarval;
5057  SCIP_Bool optimal;
5058 
5059  assert(benders != NULL);
5060  assert(set != NULL);
5061  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
5062 
5063  optimal = FALSE;
5064 
5065  auxiliaryvarval = SCIPbendersGetAuxiliaryVarVal(benders, set, sol, probnumber);
5066 
5067  SCIPsetDebugMsg(set, "Subproblem %d - Auxiliary Variable: %g Subproblem Objective: %g Reldiff: %g Soltol: %g\n",
5068  probnumber, auxiliaryvarval, SCIPbendersGetSubproblemObjval(benders, probnumber),
5069  SCIPrelDiff(SCIPbendersGetSubproblemObjval(benders, probnumber), auxiliaryvarval), benders->solutiontol);
5070 
5071  if( SCIPrelDiff(SCIPbendersGetSubproblemObjval(benders, probnumber), auxiliaryvarval) < benders->solutiontol )
5072  optimal = TRUE;
5073 
5074  return optimal;
5075 }
5076 
5077 /** returns the value of the auxiliary variable value in a master problem solution */
5079  SCIP_BENDERS* benders, /**< the benders' decomposition structure */
5080  SCIP_SET* set, /**< global SCIP settings */
5081  SCIP_SOL* sol, /**< primal CIP solution */
5082  int probnumber /**< the subproblem number */
5083  )
5084 {
5085  SCIP_VAR* auxiliaryvar;
5086 
5087  assert(benders != NULL);
5088  assert(set != NULL);
5089 
5090  auxiliaryvar = SCIPbendersGetAuxiliaryVar(benders, probnumber);
5091  assert(auxiliaryvar != NULL);
5092 
5093  return SCIPgetSolVal(set->scip, sol, auxiliaryvar);
5094 }
5095 
5096 /** Solves an independent subproblem to identify its lower bound. The lower bound is then used to update the bound on
5097  * the auxiliary variable.
5098  */
5100  SCIP_BENDERS* benders, /**< Benders' decomposition */
5101  SCIP_SET* set, /**< global SCIP settings */
5102  int probnumber, /**< the subproblem to be evaluated */
5103  SCIP_Real* lowerbound, /**< the lowerbound for the subproblem */
5104  SCIP_Bool* infeasible /**< was the subproblem found to be infeasible? */
5105  )
5106 {
5107  SCIP* subproblem;
5108  SCIP_Real dualbound;
5109  SCIP_Real memorylimit;
5110  SCIP_Real timelimit;
5111  SCIP_Longint totalnodes;
5112  int disablecutoff;
5113  int verblevel;
5114  SCIP_Bool lperror;
5115  SCIP_Bool cutoff;
5116 
5117  assert(benders != NULL);
5118  assert(set != NULL);
5119 
5120  if( benders->benderssolvesub != NULL || benders->benderssolvesubconvex != NULL )
5121  {
5122  (*lowerbound) = SCIPvarGetLbGlobal(SCIPbendersGetAuxiliaryVar(benders, probnumber));
5123  (*infeasible) = FALSE;
5124 
5125  SCIPinfoMessage(set->scip, NULL, "Benders' decomposition: a bendersSolvesub or bendersSolvesubconvex has been "
5126  "implemented. SCIPbendersComputeSubproblemLowerbound can not be executed.\n");
5127  SCIPinfoMessage(set->scip, NULL, "Set the auxiliary variable lower bound by calling "
5128  "SCIPbendersUpdateSubproblemLowerbound in bendersCreatesub. The auxiliary variable %d will remain as %g\n",
5129  probnumber, (*lowerbound));
5130 
5131  return SCIP_OKAY;
5132  }
5133  else
5134  {
5135  SCIPverbMessage(set->scip, SCIP_VERBLEVEL_FULL, NULL, "Benders' decomposition: Computing a lower bound for"
5136  " subproblem %d\n", probnumber);
5137  }
5138 
5139  /* getting the subproblem to evaluate */
5140  subproblem = SCIPbendersSubproblem(benders, probnumber);
5141 
5142  (*lowerbound) = -SCIPinfinity(subproblem);
5143  (*infeasible) = FALSE;
5144 
5145  SCIP_CALL( SCIPgetIntParam(subproblem, "display/verblevel", &verblevel) );
5146  SCIP_CALL( SCIPsetIntParam(subproblem, "display/verblevel", (int)SCIP_VERBLEVEL_NONE) );
5147 #ifdef SCIP_MOREDEBUG
5148  SCIP_CALL( SCIPsetIntParam(subproblem, "display/verblevel", (int)SCIP_VERBLEVEL_HIGH) );
5149 #endif
5150 
5151  /* copying memory and time limits */
5152  SCIP_CALL( SCIPgetRealParam(subproblem, "limits/time", &timelimit) );
5153  SCIP_CALL( SCIPgetRealParam(subproblem, "limits/memory", &memorylimit) );
5154  SCIP_CALL( copyMemoryAndTimeLimits(set->scip, subproblem) );
5155 
5156  /* if the subproblem is independent, then the default SCIP settings are used. Otherwise, only the root node is solved
5157  * to compute a lower bound on the subproblem
5158  */
5159  SCIP_CALL( SCIPgetLongintParam(subproblem, "limits/totalnodes", &totalnodes) );
5160  SCIP_CALL( SCIPgetIntParam(subproblem, "lp/disablecutoff", &disablecutoff) );
5161  if( !SCIPbendersSubproblemIsIndependent(benders, probnumber) )
5162  {
5163  SCIP_CALL( SCIPsetLongintParam(subproblem, "limits/totalnodes", 1LL) );
5164  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/disablecutoff", 1) );
5165  }
5166 
5167  /* if the subproblem not independent and is convex, then the probing LP is solved. Otherwise, the MIP is solved */
5168  dualbound = -SCIPinfinity(subproblem);
5169  if( SCIPbendersGetSubproblemType(benders, probnumber) == SCIP_BENDERSSUBTYPE_CONVEXCONT )
5170  {
5171  SCIP_Bool solvenlp = FALSE;
5172 
5173  assert(SCIPisLPConstructed(subproblem) || SCIPisNLPConstructed(subproblem));
5174 
5175  if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) > 0
5176  && SCIPbendersGetSubproblemType(benders, probnumber) <= SCIP_BENDERSSUBTYPE_CONVEXDIS )
5177  solvenlp = TRUE;
5178 
5179  SCIP_CALL( SCIPstartProbing(subproblem) );
5180  if( solvenlp )
5181  {
5182  SCIP_NLPSOLSTAT nlpsolstat;
5183  SCIP_NLPTERMSTAT nlptermstat;
5184 
5185  SCIP_CALL( SCIPsolveNLP(subproblem) ); /*lint !e666*/
5186 
5187  nlpsolstat = SCIPgetNLPSolstat(subproblem);
5188  nlptermstat = SCIPgetNLPTermstat(subproblem);
5189  SCIPdebugMsg(set->scip, "NLP solstat %d termstat %d\n", nlpsolstat, nlptermstat);
5190 
5191  if( nlptermstat == SCIP_NLPTERMSTAT_OKAY && (nlpsolstat == SCIP_NLPSOLSTAT_LOCINFEASIBLE || nlpsolstat == SCIP_NLPSOLSTAT_GLOBINFEASIBLE) )
5192  {
5193  /* trust infeasible only if terminated "okay" */
5194  (*infeasible) = TRUE;
5195  }
5196  else if( nlpsolstat == SCIP_NLPSOLSTAT_LOCOPT || nlpsolstat == SCIP_NLPSOLSTAT_GLOBOPT
5197  || nlpsolstat == SCIP_NLPSOLSTAT_FEASIBLE )
5198  {
5199  dualbound = SCIPretransformObj(subproblem, SCIPgetNLPObjval(subproblem));
5200  }
5201  }
5202  else
5203  {
5204  SCIP_CALL( SCIPsolveProbingLP(subproblem, -1, &lperror, &cutoff) );
5205 
5206  if( SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_INFEASIBLE )
5207  (*infeasible) = TRUE;
5208  else if( SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_OPTIMAL )
5209  dualbound = SCIPgetSolOrigObj(subproblem, NULL)*(int)SCIPgetObjsense(set->scip);
5210  }
5211  }
5212  else
5213  {
5214  SCIP_EVENTHDLRDATA* eventhdlrdata;
5215 
5216  /* if the subproblem is not convex, then event handlers have been added to interrupt the solve. These must be
5217  * disabled
5218  */
5220  eventhdlrdata->solvecip = TRUE;
5221 
5222  SCIP_CALL( SCIPsolve(subproblem) );
5223 
5224  if( SCIPgetStatus(subproblem) == SCIP_STATUS_INFEASIBLE )
5225  (*infeasible) = TRUE;
5226  else
5227  dualbound = SCIPgetDualbound(subproblem);
5228  }
5229 
5230  /* getting the lower bound value */
5231  (*lowerbound) = dualbound;
5232 
5233  if( !SCIPbendersSubproblemIsIndependent(benders, probnumber) )
5234  {
5235  SCIP_CALL( SCIPsetLongintParam(subproblem, "limits/totalnodes", totalnodes) );
5236  SCIP_CALL( SCIPsetIntParam(subproblem, "lp/disablecutoff", disablecutoff) );
5237  }
5238  SCIP_CALL( SCIPsetIntParam(subproblem, "display/verblevel", verblevel) );
5239  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/memory", memorylimit) );
5240  SCIP_CALL( SCIPsetRealParam(subproblem, "limits/time", timelimit) );
5241 
5242  /* the subproblem must be freed so that it is reset for the subsequent Benders' decomposition solves. If the
5243  * subproblems are independent, they are not freed. SCIPfreeBendersSubproblem must still be called, but in this
5244  * function the independent subproblems are not freed. However, they will still be freed at the end of the
5245  * solving process for the master problem.
5246  */
5247  SCIP_CALL( SCIPbendersFreeSubproblem(benders, set, probnumber) );
5248 
5249  return SCIP_OKAY;
5250 }
5251 
5252 /** Merges a subproblem into the master problem. This process just adds a copy of the subproblem variables and
5253  * constraints to the master problem, but keeps the subproblem stored in the Benders' decomposition data structure. The reason for
5254  * keeping the subproblem available is for when it is queried for solutions after the problem is solved.
5255  *
5256  * Once the subproblem is merged into the master problem, then the subproblem is flagged as disabled. This means that
5257  * it will not be solved in the subsequent subproblem solving loops.
5258  *
5259  * The associated auxiliary variables are kept in the master problem. The objective function of the merged subproblem
5260  * is added as an underestimator constraint.
5261  */
5263  SCIP_BENDERS* benders, /**< Benders' decomposition */
5264  SCIP_SET* set, /**< global SCIP settings */
5265  SCIP_HASHMAP* varmap, /**< a hashmap to store the mapping of subproblem variables corresponding
5266  * to the newly created master variables, or NULL */
5267  SCIP_HASHMAP* consmap, /**< a hashmap to store the mapping of subproblem constraints to the
5268  * corresponding newly created constraints, or NULL */
5269  int probnumber /**< the number of the subproblem that will be merged into the master problem*/
5270  )
5271 {
5272  SCIP* subproblem;
5273  SCIP_HASHMAP* localvarmap;
5274  SCIP_HASHMAP* localconsmap;
5275  SCIP_VAR** vars;
5276  SCIP_VAR* auxiliaryvar;
5277  SCIP_CONS** conss;
5278  SCIP_CONS* objcons;
5279  int nvars;
5280  int nconss;
5281  int i;
5282  SCIP_Bool uselocalvarmap;
5283  SCIP_Bool uselocalconsmap;
5284  char varname[SCIP_MAXSTRLEN];
5285  char consname[SCIP_MAXSTRLEN];
5286  const char* origvarname;
5287 
5288  assert(benders != NULL);
5289  assert(set != NULL);
5290  assert(probnumber >= 0 && probnumber < benders->nsubproblems);
5291 
5292  SCIPverbMessage(set->scip, SCIP_VERBLEVEL_HIGH, NULL, " Benders' decomposition: Infeasibility of subproblem %d can't "
5293  "be resolved. Subproblem %d is being merged into the master problem.\n", probnumber, probnumber);
5294 
5295  /* freeing the subproblem because it will be flagged as independent. Since the subproblem is flagged as independent,
5296  * it will no longer be solved or freed within the solving loop.
5297  */
5298  SCIP_CALL( SCIPbendersFreeSubproblem(benders, set, probnumber) );
5299 
5300  subproblem = SCIPbendersSubproblem(benders, probnumber);
5301 
5302  uselocalvarmap = (varmap == NULL);
5303  uselocalconsmap = (consmap == NULL);
5304 
5305  if( uselocalvarmap )
5306  {
5307  /* create the variable mapping hash map */
5308  SCIP_CALL( SCIPhashmapCreate(&localvarmap,