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