Scippy

SCIP

Solving Constraint Integer Programs

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