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