1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 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   benderscut_opt.c
26   	 * @ingroup OTHER_CFILES
27   	 * @brief  Generates a standard Benders' decomposition optimality cut
28   	 * @author Stephen J. Maher
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include "scip/pub_expr.h"
34   	#include "scip/benderscut_opt.h"
35   	#include "scip/cons_linear.h"
36   	#include "scip/pub_benderscut.h"
37   	#include "scip/pub_benders.h"
38   	#include "scip/pub_lp.h"
39   	#include "scip/pub_nlp.h"
40   	#include "scip/pub_message.h"
41   	#include "scip/pub_misc.h"
42   	#include "scip/pub_misc_linear.h"
43   	#include "scip/pub_var.h"
44   	#include "scip/scip.h"
45   	#include <string.h>
46   	
47   	#define BENDERSCUT_NAME             "optimality"
48   	#define BENDERSCUT_DESC             "Standard Benders' decomposition optimality cut"
49   	#define BENDERSCUT_PRIORITY         5000
50   	#define BENDERSCUT_LPCUT            TRUE
51   	
52   	#define SCIP_DEFAULT_ADDCUTS             FALSE  /** Should cuts be generated, instead of constraints */
53   	#define SCIP_DEFAULT_CALCMIR             TRUE   /** Should the mixed integer rounding procedure be used for the cut */
54   	
55   	/*
56   	 * Data structures
57   	 */
58   	
59   	/** Benders' decomposition cuts data */
60   	struct SCIP_BenderscutData
61   	{
62   	   SCIP_Bool             addcuts;            /**< should cuts be generated instead of constraints */
63   	   SCIP_Bool             calcmir;            /**< should the mixed integer rounding procedure be applied to cuts */
64   	};
65   	
66   	
67   	/*
68   	 * Local methods
69   	 */
70   	
71   	/** in the case of numerical troubles, the LP is resolved with solution polishing activated */
72   	static
73   	SCIP_RETCODE polishSolution(
74   	   SCIP*                 subproblem,         /**< the SCIP data structure */
75   	   SCIP_Bool*            success             /**< TRUE is the resolving of the LP was successful */
76   	   )
77   	{
78   	   int oldpolishing;
79   	   SCIP_Bool lperror;
80   	   SCIP_Bool cutoff;
81   	
82   	   assert(subproblem != NULL);
83   	   assert(SCIPinProbing(subproblem));
84   	
85   	   (*success) = FALSE;
86   	
87   	   /* setting the solution polishing parameter */
88   	   SCIP_CALL( SCIPgetIntParam(subproblem, "lp/solutionpolishing", &oldpolishing) );
89   	   SCIP_CALL( SCIPsetIntParam(subproblem, "lp/solutionpolishing", 2) );
90   	
91   	   /* resolving the probing LP */
92   	   SCIP_CALL( SCIPsolveProbingLP(subproblem, -1, &lperror, &cutoff) );
93   	
94   	   if( SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_OPTIMAL )
95   	      (*success) = TRUE;
96   	
97   	   /* resetting the solution polishing parameter */
98   	   SCIP_CALL( SCIPsetIntParam(subproblem, "lp/solutionpolishing", oldpolishing) );
99   	
100  	   return SCIP_OKAY;
101  	}
102  	
103  	/** verifying the activity of the cut when master variables are within epsilon of their upper or lower bounds
104  	 *
105  	 *  When setting up the Benders' decomposition subproblem, master variables taking values that are within epsilon
106  	 *  greater than their upper bound or less than their lower bound are set to their upper and lower bounds respectively.
107  	 *  As such, there can be a difference between the subproblem dual solution objective and the optimality cut activity,
108  	 *  when computed using the master problem solution directly. This check is to verify whether this difference is an
109  	 *  actual error or due to the violation of the upper and lower bounds when setting up the Benders' decomposition
110  	 *  subproblem.
111  	 */
112  	static
113  	SCIP_RETCODE checkSetupTolerances(
114  	   SCIP*                 masterprob,         /**< the SCIP data structure */
115  	   SCIP_SOL*             sol,                /**< the master problem solution */
116  	   SCIP_VAR**            vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
117  	   SCIP_Real*            vals,               /**< pointer to array of coefficients of the variables in the generated cut */
118  	   SCIP_Real             lhs,                /**< the left hand side of the cut */
119  	   SCIP_Real             checkobj,           /**< the objective of the subproblem computed from the dual solution */
120  	   int                   nvars,              /**< the number of variables in the cut */
121  	   SCIP_Bool*            valid               /**< returns true is the cut is valid */
122  	   )
123  	{
124  	   SCIP_Real verifyobj;
125  	   int i;
126  	
127  	   assert(masterprob != NULL);
128  	   assert(vars != NULL);
129  	   assert(vals != NULL);
130  	
131  	   /* initialising the verify objective with the left hand side of the optimality cut */
132  	   verifyobj = lhs;
133  	
134  	   /* computing the activity of the cut from the master solution and the constraint values */
135  	   for( i = 0; i < nvars; i++ )
136  	   {
137  	      SCIP_Real solval;
138  	
139  	      solval = SCIPgetSolVal(masterprob, sol, vars[i]);
140  	
141  	      /* checking whether the solution value is less than or greater than the variable bounds */
142  	      if( !SCIPisLT(masterprob, solval, SCIPvarGetUbLocal(vars[i])) )
143  	         solval = SCIPvarGetUbLocal(vars[i]);
144  	      else if( !SCIPisGT(masterprob, solval, SCIPvarGetLbLocal(vars[i])) )
145  	         solval = SCIPvarGetLbLocal(vars[i]);
146  	
147  	      verifyobj -= solval*vals[i];
148  	   }
149  	
150  	   (*valid) = SCIPisFeasEQ(masterprob, checkobj, verifyobj);
151  	
152  	   return SCIP_OKAY;
153  	}
154  	
155  	/** when solving NLP subproblems, numerical issues are addressed by tightening the feasibility tolerance */
156  	static
157  	SCIP_RETCODE resolveNLPWithTighterFeastol(
158  	   SCIP*                 subproblem,         /**< the SCIP data structure */
159  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
160  	   SCIP_Real             multiplier,         /**< the amount by which to decrease the tolerance */
161  	   SCIP_Bool*            success             /**< TRUE is the resolving of the LP was successful */
162  	   )
163  	{
164  	   SCIP_NLPSOLSTAT nlpsolstat;
165  	#ifdef SCIP_DEBUG
166  	   SCIP_NLPTERMSTAT nlptermstat;
167  	#endif
168  	   SCIP_NLPPARAM nlpparam = SCIPbendersGetNLPParam(benders);
169  	#ifdef SCIP_MOREDEBUG
170  	   SCIP_SOL* nlpsol;
171  	#endif
172  	
173  	   assert(subproblem != NULL);
174  	   assert(SCIPinProbing(subproblem));
175  	
176  	   (*success) = FALSE;
177  	
178  	   /* reduce the default feasibility and optimality tolerance by given factor (typically 0.01) */
179  	   nlpparam.feastol *= multiplier;
180  	   nlpparam.opttol *= multiplier;
181  	
182  	   SCIP_CALL( SCIPsolveNLPParam(subproblem, nlpparam) );
183  	
184  	   nlpsolstat = SCIPgetNLPSolstat(subproblem);
185  	#ifdef SCIP_DEBUG
186  	   nlptermstat = SCIPgetNLPTermstat(subproblem);
187  	   SCIPdebugMsg(subproblem, "NLP solstat %d termstat %d\n", nlpsolstat, nlptermstat);
188  	#endif
189  	
190  	   if( nlpsolstat == SCIP_NLPSOLSTAT_LOCOPT || nlpsolstat == SCIP_NLPSOLSTAT_GLOBOPT
191  	      || nlpsolstat == SCIP_NLPSOLSTAT_FEASIBLE )
192  	   {
193  	#ifdef SCIP_MOREDEBUG
194  	      SCIP_CALL( SCIPcreateNLPSol(subproblem, &nlpsol, NULL) );
195  	      SCIP_CALL( SCIPprintSol(subproblem, nlpsol, NULL, FALSE) );
196  	      SCIP_CALL( SCIPfreeSol(subproblem, &nlpsol) );
197  	#endif
198  	
199  	      (*success) = TRUE;
200  	   }
201  	
202  	   return SCIP_OKAY;
203  	}
204  	
205  	/** adds a variable and value to the constraint/row arrays */
206  	static
207  	SCIP_RETCODE addVariableToArray(
208  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
209  	   SCIP_VAR***           vars,               /**< pointer to the array of variables in the generated cut with non-zero coefficient */
210  	   SCIP_Real**           vals,               /**< pointer to the array of coefficients of the variables in the generated cut */
211  	   SCIP_VAR*             addvar,             /**< the variable that will be added to the array */
212  	   SCIP_Real             addval,             /**< the value that will be added to the array */
213  	   int*                  nvars,              /**< the number of variables in the variable array */
214  	   int*                  varssize            /**< the length of the variable size */
215  	   )
216  	{
217  	   assert(masterprob != NULL);
218  	   assert(vars != NULL);
219  	   assert(*vars != NULL);
220  	   assert(vals != NULL);
221  	   assert(*vals != NULL);
222  	   assert(addvar != NULL);
223  	   assert(nvars != NULL);
224  	   assert(varssize != NULL);
225  	
226  	   if( *nvars >= *varssize )
227  	   {
228  	      *varssize = SCIPcalcMemGrowSize(masterprob, *varssize + 1);
229  	      SCIP_CALL( SCIPreallocBufferArray(masterprob, vars, *varssize) );
230  	      SCIP_CALL( SCIPreallocBufferArray(masterprob, vals, *varssize) );
231  	   }
232  	   assert(*nvars < *varssize);
233  	
234  	   (*vars)[*nvars] = addvar;
235  	   (*vals)[*nvars] = addval;
236  	   (*nvars)++;
237  	
238  	   return SCIP_OKAY;
239  	}
240  	
241  	/** returns the variable solution either from the NLP or from the primal vals array */
242  	static
243  	SCIP_Real getNlpVarSol(
244  	   SCIP_VAR*             var,                /**< the variable for which the solution is requested */
245  	   SCIP_Real*            primalvals,         /**< the primal solutions for the NLP, can be NULL */
246  	   SCIP_HASHMAP*         var2idx             /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */
247  	   )
248  	{
249  	   SCIP_Real varsol;
250  	   int idx;
251  	
252  	   assert(var != NULL);
253  	   assert((primalvals == NULL && var2idx == NULL) || (primalvals != NULL && var2idx != NULL));
254  	
255  	   if( var2idx != NULL && primalvals != NULL )
256  	   {
257  	      assert(SCIPhashmapExists(var2idx, (void*)var) );
258  	      idx = SCIPhashmapGetImageInt(var2idx, (void*)var);
259  	      varsol = primalvals[idx];
260  	   }
261  	   else
262  	      varsol = SCIPvarGetNLPSol(var);
263  	
264  	   return varsol;
265  	}
266  	
267  	/** calculates a MIR cut from the coefficients of the standard optimality cut */
268  	static
269  	SCIP_RETCODE computeMIRForOptimalityCut(
270  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
271  	   SCIP_SOL*             sol,                /**< primal CIP solution */
272  	   SCIP_VAR**            vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
273  	   SCIP_Real*            vals,               /**< pointer to array of coefficients of the variables in the generated cut */
274  	   SCIP_Real             lhs,                /**< the left hand side of the cut */
275  	   SCIP_Real             rhs,                /**< the right hand side of the cut */
276  	   int                   nvars,              /**< the number of variables in the cut */
277  	   SCIP_Real*            cutcoefs,           /**< the coefficients of the MIR cut */
278  	   int*                  cutinds,            /**< the variable indices of the MIR cut */
279  	   SCIP_Real*            cutrhs,             /**< the RHS of the MIR cut */
280  	   int*                  cutnnz,             /**< the number of non-zeros in the cut */
281  	   SCIP_Bool*            success             /**< was the MIR cut successfully computed? */
282  	   )
283  	{
284  	   SCIP_AGGRROW* aggrrow;
285  	   SCIP_Real* rowvals;
286  	   int* rowinds;
287  	
288  	   SCIP_Real cutefficacy;
289  	   int cutrank;
290  	   SCIP_Bool cutislocal;
291  	
292  	   SCIP_Bool cutsuccess;
293  	
294  	   int i;
295  	
296  	   /* creating the aggregation row. There will be only a single row in this aggregation, since it is only used to
297  	    * compute the MIR coefficients
298  	    */
299  	   SCIP_CALL( SCIPaggrRowCreate(masterprob, &aggrrow) );
300  	
301  	   /* retrieving the indices for the variables in the optimality cut. All of the values must be negated, since the
302  	    * aggregation row requires a RHS, where the optimality cut is computed with an LHS
303  	    */
304  	   SCIP_CALL( SCIPallocBufferArray(masterprob, &rowvals, nvars) );
305  	   SCIP_CALL( SCIPallocBufferArray(masterprob, &rowinds, nvars) );
306  	
307  	   assert(SCIPisInfinity(masterprob, rhs));
308  	   assert(!SCIPisInfinity(masterprob, lhs));
309  	   for( i = 0; i < nvars; i++ )
310  	   {
311  	      rowinds[i] = SCIPvarGetProbindex(vars[i]);
312  	      rowvals[i] = -vals[i];
313  	   }
314  	
315  	   /* adding the optimality cut to the aggregation row */
316  	   SCIP_CALL( SCIPaggrRowAddCustomCons(masterprob, aggrrow, rowinds, rowvals, nvars, -lhs, 1.0, 1, FALSE) );
317  	
318  	   /* calculating a flow cover for the optimality cut */
319  	   SCIP_CALL( SCIPcalcFlowCover(masterprob, sol, TRUE, 0.9999, FALSE, aggrrow, cutcoefs, cutrhs, cutinds, cutnnz,
320  	         &cutefficacy, NULL, &cutislocal, &cutsuccess) );
321  	   (*success) = cutsuccess;
322  	
323  	   /* calculating the MIR coefficients for the optimality cut */
324  	   SCIP_CALL( SCIPcalcMIR(masterprob, sol, TRUE, 0.9999, TRUE, FALSE, FALSE, NULL, NULL, 0.001, 0.999, 1.0, aggrrow,
325  	         cutcoefs, cutrhs, cutinds, cutnnz, &cutefficacy, &cutrank, &cutislocal, &cutsuccess) );
326  	   (*success) = ((*success) || cutsuccess);
327  	
328  	   /* the cut is only successful if the efficacy is high enough */
329  	   (*success) = (*success) && SCIPisEfficacious(masterprob, cutefficacy);
330  	
331  	   /* try to tighten the coefficients of the cut */
332  	   if( (*success) )
333  	   {
334  	      SCIP_Bool redundant;
335  	      int nchgcoefs;
336  	
337  	      redundant = SCIPcutsTightenCoefficients(masterprob, FALSE, cutcoefs, cutrhs, cutinds, cutnnz, &nchgcoefs);
338  	
339  	      (*success) = !redundant;
340  	   }
341  	
342  	   /* freeing the local memory */
343  	   SCIPfreeBufferArray(masterprob, &rowinds);
344  	   SCIPfreeBufferArray(masterprob, &rowvals);
345  	   SCIPaggrRowFree(masterprob, &aggrrow);
346  	
347  	   return SCIP_OKAY;
348  	}
349  	
350  	/** computes a standard Benders' optimality cut from the dual solutions of the LP */
351  	static
352  	SCIP_RETCODE computeStandardLPOptimalityCut(
353  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
354  	   SCIP*                 subproblem,         /**< the SCIP instance of the subproblem */
355  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
356  	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
357  	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
358  	   SCIP_Real*            lhs,                /**< the left hand side of the cut */
359  	   SCIP_Real*            rhs,                /**< the right hand side of the cut */
360  	   int*                  nvars,              /**< the number of variables in the cut */
361  	   int*                  varssize,           /**< the number of variables in the array */
362  	   SCIP_Real*            checkobj,           /**< stores the objective function computed from the dual solution */
363  	   SCIP_Bool*            success             /**< was the cut generation successful? */
364  	   )
365  	{
366  	   SCIP_VAR** subvars;
367  	   SCIP_VAR** fixedvars;
368  	   int nsubvars;
369  	   int nfixedvars;
370  	   SCIP_Real dualsol;
371  	   SCIP_Real addval;
372  	   int nrows;
373  	   int i;
374  	
375  	   (*checkobj) = 0;
376  	
377  	   assert(masterprob != NULL);
378  	   assert(subproblem != NULL);
379  	   assert(benders != NULL);
380  	   assert(vars != NULL);
381  	   assert(*vars != NULL);
382  	   assert(vals != NULL);
383  	   assert(*vals != NULL);
384  	
385  	   (*success) = FALSE;
386  	
387  	   /* looping over all LP rows and setting the coefficients of the cut */
388  	   nrows = SCIPgetNLPRows(subproblem);
389  	   for( i = 0; i < nrows; i++ )
390  	   {
391  	      SCIP_ROW* lprow;
392  	
393  	      lprow = SCIPgetLPRows(subproblem)[i];
394  	      assert(lprow != NULL);
395  	
396  	      dualsol = SCIProwGetDualsol(lprow);
397  	      assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) );
398  	
399  	      if( SCIPisZero(subproblem, dualsol) )
400  	         continue;
401  	
402  	      if( dualsol > 0.0 )
403  	         addval = dualsol*SCIProwGetLhs(lprow);
404  	      else
405  	         addval = dualsol*SCIProwGetRhs(lprow);
406  	
407  	      (*lhs) += addval;
408  	
409  	      /* if the bound becomes infinite, then the cut generation terminates. */
410  	      if( SCIPisInfinity(masterprob, (*lhs)) || SCIPisInfinity(masterprob, -(*lhs))
411  	         || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval))
412  	      {
413  	         (*success) = FALSE;
414  	         SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut. lhs = %g addval = %g.\n", (*lhs), addval);
415  	         return SCIP_OKAY;
416  	      }
417  	   }
418  	
419  	   nsubvars = SCIPgetNVars(subproblem);
420  	   subvars = SCIPgetVars(subproblem);
421  	   nfixedvars = SCIPgetNFixedVars(subproblem);
422  	   fixedvars = SCIPgetFixedVars(subproblem);
423  	
424  	   /* looping over all variables to update the coefficients in the computed cut. */
425  	   for( i = 0; i < nsubvars + nfixedvars; i++ )
426  	   {
427  	      SCIP_VAR* var;
428  	      SCIP_VAR* mastervar;
429  	      SCIP_Real redcost;
430  	
431  	      if( i < nsubvars )
432  	         var = subvars[i];
433  	      else
434  	         var = fixedvars[i - nsubvars];
435  	
436  	      /* retrieving the master problem variable for the given subproblem variable. */
437  	      SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) );
438  	
439  	      redcost = SCIPgetVarRedcost(subproblem, var);
440  	
441  	      (*checkobj) += SCIPvarGetUnchangedObj(var)*SCIPvarGetSol(var, TRUE);
442  	
443  	      /* checking whether the subproblem variable has a corresponding master variable. */
444  	      if( mastervar != NULL )
445  	      {
446  	         SCIP_Real coef;
447  	
448  	         coef = -1.0*(SCIPvarGetObj(var) + redcost);
449  	
450  	         if( !SCIPisZero(masterprob, coef) )
451  	         {
452  	            /* adding the variable to the storage */
453  	            SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) );
454  	         }
455  	      }
456  	      else
457  	      {
458  	         if( !SCIPisZero(subproblem, redcost) )
459  	         {
460  	            addval = 0;
461  	
462  	            if( SCIPisPositive(subproblem, redcost) )
463  	               addval = redcost*SCIPvarGetLbLocal(var);
464  	            else if( SCIPisNegative(subproblem, redcost) )
465  	               addval = redcost*SCIPvarGetUbLocal(var);
466  	
467  	            (*lhs) += addval;
468  	
469  	            /* if the bound becomes infinite, then the cut generation terminates. */
470  	            if( SCIPisInfinity(masterprob, (*lhs)) || SCIPisInfinity(masterprob, -(*lhs))
471  	               || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval))
472  	            {
473  	               (*success) = FALSE;
474  	               SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut.\n");
475  	               return SCIP_OKAY;
476  	            }
477  	         }
478  	      }
479  	   }
480  	
481  	   assert(SCIPisInfinity(masterprob, (*rhs)));
482  	   /* the rhs should be infinite. If it changes, then there is an error */
483  	   if( !SCIPisInfinity(masterprob, (*rhs)) )
484  	   {
485  	      (*success) = FALSE;
486  	      SCIPdebugMsg(masterprob, "RHS is not infinite. rhs = %g.\n", (*rhs));
487  	      return SCIP_OKAY;
488  	   }
489  	
490  	   (*success) = TRUE;
491  	
492  	   return SCIP_OKAY;
493  	}
494  	
495  	/** computes a standard Benders' optimality cut from the dual solutions of the NLP */
496  	static
497  	SCIP_RETCODE computeStandardNLPOptimalityCut(
498  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
499  	   SCIP*                 subproblem,         /**< the SCIP instance of the subproblem */
500  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
501  	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
502  	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
503  	   SCIP_Real*            lhs,                /**< the left hand side of the cut */
504  	   SCIP_Real*            rhs,                /**< the right hand side of the cut */
505  	   int*                  nvars,              /**< the number of variables in the cut */
506  	   int*                  varssize,           /**< the number of variables in the array */
507  	   SCIP_Real             objective,          /**< the objective function of the subproblem */
508  	   SCIP_Real*            primalvals,         /**< the primal solutions for the NLP, can be NULL */
509  	   SCIP_Real*            consdualvals,       /**< dual variables for the constraints, can be NULL */
510  	   SCIP_Real*            varlbdualvals,      /**< the dual variables for the variable lower bounds, can be NULL */
511  	   SCIP_Real*            varubdualvals,      /**< the dual variables for the variable upper bounds, can be NULL */
512  	   SCIP_HASHMAP*         row2idx,            /**< mapping between the row in the subproblem to the index in the dual array, can be NULL */
513  	   SCIP_HASHMAP*         var2idx,            /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */
514  	   SCIP_Real*            checkobj,           /**< stores the objective function computed from the dual solution */
515  	   SCIP_Bool*            success             /**< was the cut generation successful? */
516  	   )
517  	{
518  	   SCIP_VAR** subvars;
519  	   SCIP_VAR** fixedvars;
520  	   int nsubvars;
521  	   int nfixedvars;
522  	   SCIP_Real dirderiv;
523  	   SCIP_Real dualsol;
524  	   int nrows;
525  	   int idx;
526  	   int i;
527  	
528  	   (*checkobj) = 0;
529  	
530  	   assert(masterprob != NULL);
531  	   assert(subproblem != NULL);
532  	   assert(benders != NULL);
533  	   assert(SCIPisNLPConstructed(subproblem));
534  	   assert(SCIPgetNLPSolstat(subproblem) <= SCIP_NLPSOLSTAT_FEASIBLE || consdualvals != NULL);
535  	   assert(SCIPhasNLPSolution(subproblem) || consdualvals != NULL);
536  	
537  	   (*success) = FALSE;
538  	
539  	   if( !(primalvals == NULL && consdualvals == NULL && varlbdualvals == NULL && varubdualvals == NULL && row2idx == NULL && var2idx == NULL)
540  	      && !(primalvals != NULL && consdualvals != NULL && varlbdualvals != NULL && varubdualvals != NULL && row2idx != NULL && var2idx != NULL) ) /*lint !e845*/
541  	   {
542  	      SCIPerrorMessage("The optimality cut must generated from either a SCIP instance or all of the dual solutions and indices must be supplied");
543  	      (*success) = FALSE;
544  	
545  	      return SCIP_ERROR;
546  	   }
547  	
548  	   nsubvars = SCIPgetNNLPVars(subproblem);
549  	   subvars = SCIPgetNLPVars(subproblem);
550  	   nfixedvars = SCIPgetNFixedVars(subproblem);
551  	   fixedvars = SCIPgetFixedVars(subproblem);
552  	
553  	   /* our optimality cut implementation assumes that SCIP did not modify the objective function and sense,
554  	    * that is, that the objective function value of the NLP corresponds to the value of the auxiliary variable
555  	    * if that wouldn't be the case, then the scaling and offset may have to be considered when adding the
556  	    * auxiliary variable to the cut (cons/row)?
557  	    */
558  	   assert(SCIPgetTransObjoffset(subproblem) == 0.0);
559  	   assert(SCIPgetTransObjscale(subproblem) == 1.0);
560  	   assert(SCIPgetObjsense(subproblem) == SCIP_OBJSENSE_MINIMIZE);
561  	
562  	   (*lhs) = objective;
563  	   assert(!SCIPisInfinity(subproblem, REALABS(*lhs)));
564  	
565  	   (*rhs) = SCIPinfinity(masterprob);
566  	
567  	   dirderiv = 0.0;
568  	
569  	   /* looping over all NLP rows and setting the corresponding coefficients of the cut */
570  	   nrows = SCIPgetNNLPNlRows(subproblem);
571  	   for( i = 0; i < nrows; i++ )
572  	   {
573  	      SCIP_NLROW* nlrow;
574  	
575  	      nlrow = SCIPgetNLPNlRows(subproblem)[i];
576  	      assert(nlrow != NULL);
577  	
578  	      if( row2idx != NULL && consdualvals != NULL )
579  	      {
580  	         assert(SCIPhashmapExists(row2idx, (void*)nlrow) );
581  	         idx = SCIPhashmapGetImageInt(row2idx, (void*)nlrow);
582  	         dualsol = consdualvals[idx];
583  	      }
584  	      else
585  	         dualsol = SCIPnlrowGetDualsol(nlrow);
586  	      assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) );
587  	
588  	      if( SCIPisZero(subproblem, dualsol) )
589  	         continue;
590  	
591  	      SCIP_CALL( SCIPaddNlRowGradientBenderscutOpt(masterprob, subproblem, benders, nlrow,
592  	            -dualsol, primalvals, var2idx, &dirderiv, vars, vals, nvars, varssize) );
593  	   }
594  	
595  	   /* looping over sub- and fixed variables to compute checkobj */
596  	   for( i = 0; i < nsubvars; i++ )
597  	      (*checkobj) += SCIPvarGetObj(subvars[i]) * getNlpVarSol(subvars[i], primalvals, var2idx);
598  	
599  	   for( i = 0; i < nfixedvars; i++ )
600  	      *checkobj += SCIPvarGetUnchangedObj(fixedvars[i]) * getNlpVarSol(fixedvars[i], primalvals, var2idx);
601  	
602  	   *lhs += dirderiv;
603  	
604  	   /* if the side became infinite or dirderiv was infinite, then the cut generation terminates. */
605  	   if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs)
606  	      || SCIPisInfinity(masterprob, dirderiv) || SCIPisInfinity(masterprob, -dirderiv))
607  	   {
608  	      (*success) = FALSE;
609  	      SCIPdebugMsg(masterprob, "Infinite bound when generating optimality cut. lhs = %g dirderiv = %g.\n", *lhs, dirderiv);
610  	      return SCIP_OKAY;
611  	   }
612  	
613  	   (*success) = TRUE;
614  	
615  	   return SCIP_OKAY;
616  	}
617  	
618  	
619  	/** Adds the auxiliary variable to the generated cut. If this is the first optimality cut for the subproblem, then the
620  	 *  auxiliary variable is first created and added to the master problem.
621  	 */
622  	static
623  	SCIP_RETCODE addAuxiliaryVariableToCut(
624  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
625  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
626  	   SCIP_VAR**            vars,               /**< the variables in the generated cut with non-zero coefficient */
627  	   SCIP_Real*            vals,               /**< the coefficients of the variables in the generated cut */
628  	   int*                  nvars,              /**< the number of variables in the cut */
629  	   int                   probnumber          /**< the number of the pricing problem */
630  	   )
631  	{
632  	   SCIP_VAR* auxiliaryvar;
633  	
634  	   assert(masterprob != NULL);
635  	   assert(benders != NULL);
636  	   assert(vars != NULL);
637  	   assert(vals != NULL);
638  	
639  	   auxiliaryvar = SCIPbendersGetAuxiliaryVar(benders, probnumber);
640  	
641  	   vars[(*nvars)] = auxiliaryvar;
642  	   vals[(*nvars)] = 1.0;
643  	   (*nvars)++;
644  	
645  	   return SCIP_OKAY;
646  	}
647  	
648  	
649  	/*
650  	 * Callback methods of Benders' decomposition cuts
651  	 */
652  	
653  	/** destructor of Benders' decomposition cuts to free user data (called when SCIP is exiting) */
654  	static
655  	SCIP_DECL_BENDERSCUTFREE(benderscutFreeOpt)
656  	{  /*lint --e{715}*/
657  	   SCIP_BENDERSCUTDATA* benderscutdata;
658  	
659  	   assert( benderscut != NULL );
660  	   assert( strcmp(SCIPbenderscutGetName(benderscut), BENDERSCUT_NAME) == 0 );
661  	
662  	   /* free Benders' cut data */
663  	   benderscutdata = SCIPbenderscutGetData(benderscut);
664  	   assert( benderscutdata != NULL );
665  	
666  	   SCIPfreeBlockMemory(scip, &benderscutdata);
667  	
668  	   SCIPbenderscutSetData(benderscut, NULL);
669  	
670  	   return SCIP_OKAY;
671  	}
672  	
673  	
674  	/** execution method of Benders' decomposition cuts */
675  	static
676  	SCIP_DECL_BENDERSCUTEXEC(benderscutExecOpt)
677  	{  /*lint --e{715}*/
678  	   SCIP* subproblem;
679  	   SCIP_BENDERSCUTDATA* benderscutdata;
680  	   SCIP_Bool nlprelaxation;
681  	   SCIP_Bool addcut;
682  	   char cutname[SCIP_MAXSTRLEN];
683  	
684  	   assert(scip != NULL);
685  	   assert(benders != NULL);
686  	   assert(benderscut != NULL);
687  	   assert(result != NULL);
688  	   assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
689  	
690  	   /* retrieving the Benders' cut data */
691  	   benderscutdata = SCIPbenderscutGetData(benderscut);
692  	
693  	   /* if the cuts are generated prior to the solving stage, then rows can not be generated. So constraints must be
694  	    * added to the master problem.
695  	    */
696  	   if( SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE )
697  	      addcut = FALSE;
698  	   else
699  	      addcut = benderscutdata->addcuts;
700  	
701  	   /* setting the name of the generated cut */
702  	   (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "optimalitycut_%d_%" SCIP_LONGINT_FORMAT, probnumber,
703  	      SCIPbenderscutGetNFound(benderscut) );
704  	
705  	   subproblem = SCIPbendersSubproblem(benders, probnumber);
706  	
707  	   if( subproblem == NULL )
708  	   {
709  	      SCIPdebugMsg(scip, "The subproblem %d is set to NULL. The <%s> Benders' decomposition cut can not be executed.\n",
710  	         probnumber, BENDERSCUT_NAME);
711  	
712  	      (*result) = SCIP_DIDNOTRUN;
713  	      return SCIP_OKAY;
714  	   }
715  	
716  	   /* setting a flag to indicate whether the NLP relaxation should be used to generate cuts */
717  	   nlprelaxation = SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem);
718  	
719  	   /* only generate optimality cuts if the subproblem LP or NLP is optimal,
720  	    * since we use the dual solution of the LP/NLP to construct the optimality cut
721  	    */
722  	   if( SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING &&
723  	      ((!nlprelaxation && SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_OPTIMAL) ||
724  	       (nlprelaxation && SCIPgetNLPSolstat(subproblem) <= SCIP_NLPSOLSTAT_FEASIBLE)) )
725  	   {
726  	      /* generating a cut for a given subproblem */
727  	      SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname,
728  	            SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut,
729  	            FALSE, result) );
730  	
731  	      /* if it was not possible to generate a cut, this could be due to numerical issues. So the solution to the LP is
732  	       * resolved and the generation of the cut is reattempted. For NLPs, we do not have such a polishing yet.
733  	       */
734  	      if( (*result) == SCIP_DIDNOTFIND )
735  	      {
736  	         SCIP_Bool success;
737  	
738  	         SCIPdebugMsg(scip, "Numerical trouble generating optimality cut for subproblem %d.\n", probnumber);
739  	
740  	         if( !nlprelaxation )
741  	         {
742  	            SCIPdebugMsg(scip, "Attempting to polish the LP solution to find an alternative dual extreme point.\n");
743  	
744  	            SCIP_CALL( polishSolution(subproblem, &success) );
745  	
746  	            /* only attempt to generate a cut if the solution polishing was successful */
747  	            if( success )
748  	            {
749  	               SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname,
750  	                     SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut,
751  	                     FALSE, result) );
752  	            }
753  	         }
754  	         else
755  	         {
756  	            SCIP_Real multiplier = 0.01;
757  	
758  	            SCIPdebugMsg(scip, "Attempting to resolve the NLP with a tighter feasibility tolerance to find an "
759  	               "alternative dual extreme point.\n");
760  	
761  	            while( multiplier > 1e-06 && (*result) == SCIP_DIDNOTFIND )
762  	            {
763  	               SCIP_CALL( resolveNLPWithTighterFeastol(subproblem, benders, multiplier, &success) );
764  	
765  	               if( success )
766  	               {
767  	                  SCIP_CALL( SCIPgenerateAndApplyBendersOptCut(scip, subproblem, benders, benderscut, sol, probnumber, cutname,
768  	                        SCIPbendersGetSubproblemObjval(benders, probnumber), NULL, NULL, NULL, NULL, NULL, NULL, type, addcut,
769  	                        FALSE, result) );
770  	               }
771  	
772  	               multiplier *= 0.1;
773  	            }
774  	         }
775  	      }
776  	   }
777  	
778  	   return SCIP_OKAY;
779  	}
780  	
781  	
782  	/*
783  	 * Benders' decomposition cuts specific interface methods
784  	 */
785  	
786  	/** creates the opt Benders' decomposition cuts and includes it in SCIP */
787  	SCIP_RETCODE SCIPincludeBenderscutOpt(
788  	   SCIP*                 scip,               /**< SCIP data structure */
789  	   SCIP_BENDERS*         benders             /**< Benders' decomposition */
790  	   )
791  	{
792  	   SCIP_BENDERSCUTDATA* benderscutdata;
793  	   SCIP_BENDERSCUT* benderscut;
794  	   char paramname[SCIP_MAXSTRLEN];
795  	
796  	   assert(benders != NULL);
797  	
798  	   /* create opt Benders' decomposition cuts data */
799  	   SCIP_CALL( SCIPallocBlockMemory(scip, &benderscutdata) );
800  	
801  	   benderscut = NULL;
802  	
803  	   /* include Benders' decomposition cuts */
804  	   SCIP_CALL( SCIPincludeBenderscutBasic(scip, benders, &benderscut, BENDERSCUT_NAME, BENDERSCUT_DESC,
805  	         BENDERSCUT_PRIORITY, BENDERSCUT_LPCUT, benderscutExecOpt, benderscutdata) );
806  	
807  	   assert(benderscut != NULL);
808  	
809  	   /* setting the non fundamental callbacks via setter functions */
810  	   SCIP_CALL( SCIPsetBenderscutFree(scip, benderscut, benderscutFreeOpt) );
811  	
812  	   /* add opt Benders' decomposition cuts parameters */
813  	   (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/benderscut/%s/addcuts",
814  	      SCIPbendersGetName(benders), BENDERSCUT_NAME);
815  	   SCIP_CALL( SCIPaddBoolParam(scip, paramname,
816  	         "should cuts be generated and added to the cutpool instead of global constraints directly added to the problem.",
817  	         &benderscutdata->addcuts, FALSE, SCIP_DEFAULT_ADDCUTS, NULL, NULL) );
818  	
819  	   (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "benders/%s/benderscut/%s/mir",
820  	      SCIPbendersGetName(benders), BENDERSCUT_NAME);
821  	   SCIP_CALL( SCIPaddBoolParam(scip, paramname,
822  	         "should the mixed integer rounding procedure be applied to cuts",
823  	         &benderscutdata->calcmir, FALSE, SCIP_DEFAULT_CALCMIR, NULL, NULL) );
824  	
825  	   return SCIP_OKAY;
826  	}
827  	
828  	/** Generates a classical Benders' optimality cut using the dual solutions from the subproblem or the input arrays. If
829  	 *  the dual solutions are input as arrays, then a mapping between the array indices and the rows/variables is required.
830  	 *  As a cut strengthening approach, when an optimality cut is being generated (i.e. not for feasibility cuts) a MIR
831  	 *  procedure is performed on the row. This procedure attempts to find a stronger constraint, if this doesn't happen,
832  	 *  then the original constraint is added to SCIP.
833  	 *
834  	 *  This method can also be used to generate a feasibility cut, if a problem to minimise the infeasibilities has been solved
835  	 *  to generate the dual solutions
836  	 */
837  	SCIP_RETCODE SCIPgenerateAndApplyBendersOptCut(
838  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
839  	   SCIP*                 subproblem,         /**< the SCIP instance of the pricing problem */
840  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition */
841  	   SCIP_BENDERSCUT*      benderscut,         /**< the benders' decomposition cut method */
842  	   SCIP_SOL*             sol,                /**< primal CIP solution */
843  	   int                   probnumber,         /**< the number of the pricing problem */
844  	   char*                 cutname,            /**< the name for the cut to be generated */
845  	   SCIP_Real             objective,          /**< the objective function of the subproblem */
846  	   SCIP_Real*            primalvals,         /**< the primal solutions for the NLP, can be NULL */
847  	   SCIP_Real*            consdualvals,       /**< dual variables for the constraints, can be NULL */
848  	   SCIP_Real*            varlbdualvals,      /**< the dual variables for the variable lower bounds, can be NULL */
849  	   SCIP_Real*            varubdualvals,      /**< the dual variables for the variable upper bounds, can be NULL */
850  	   SCIP_HASHMAP*         row2idx,            /**< mapping between the row in the subproblem to the index in the dual array, can be NULL */
851  	   SCIP_HASHMAP*         var2idx,            /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */
852  	   SCIP_BENDERSENFOTYPE  type,               /**< the enforcement type calling this function */
853  	   SCIP_Bool             addcut,             /**< should the Benders' cut be added as a cut or constraint */
854  	   SCIP_Bool             feasibilitycut,     /**< is this called for the generation of a feasibility cut */
855  	   SCIP_RESULT*          result              /**< the result from solving the subproblems */
856  	   )
857  	{
858  	   SCIP_CONSHDLR* consbenders;
859  	   SCIP_CONS* cons;
860  	   SCIP_ROW* row;
861  	   SCIP_VAR** vars;
862  	   SCIP_Real* vals;
863  	   SCIP_Real lhs;
864  	   SCIP_Real rhs;
865  	   int nvars;
866  	   int varssize;
867  	   int nmastervars;
868  	   SCIP_Bool calcmir;
869  	   SCIP_Bool optimal;
870  	   SCIP_Bool success;
871  	   SCIP_Bool mirsuccess;
872  	
873  	   SCIP_Real checkobj;
874  	   SCIP_Real verifyobj;
875  	
876  	   assert(masterprob != NULL);
877  	   assert(subproblem != NULL);
878  	   assert(benders != NULL);
879  	   assert(benderscut != NULL);
880  	   assert(result != NULL);
881  	   assert((primalvals == NULL && consdualvals == NULL && varlbdualvals == NULL && varubdualvals == NULL
882  	         && row2idx == NULL && var2idx == NULL)
883  	      || (primalvals != NULL && consdualvals != NULL && varlbdualvals != NULL && varubdualvals != NULL
884  	         && row2idx != NULL && var2idx != NULL));
885  	
886  	   row = NULL;
887  	   cons = NULL;
888  	
889  	   calcmir = SCIPbenderscutGetData(benderscut)->calcmir && SCIPgetStage(masterprob) >= SCIP_STAGE_INITSOLVE && SCIPgetSubscipDepth(masterprob) == 0;
890  	   success = FALSE;
891  	   mirsuccess = FALSE;
892  	
893  	   /* retrieving the Benders' decomposition constraint handler */
894  	   consbenders = SCIPfindConshdlr(masterprob, "benders");
895  	
896  	   /* checking the optimality of the original problem with a comparison between the auxiliary variable and the
897  	    * objective value of the subproblem */
898  	   if( feasibilitycut )
899  	      optimal = FALSE;
900  	   else
901  	   {
902  	      SCIP_CALL( SCIPcheckBendersSubproblemOptimality(masterprob, benders, sol, probnumber, &optimal) );
903  	   }
904  	
905  	   if( optimal )
906  	   {
907  	      (*result) = SCIP_FEASIBLE;
908  	      SCIPdebugMsg(masterprob, "No cut added for subproblem %d\n", probnumber);
909  	      return SCIP_OKAY;
910  	   }
911  	
912  	   /* allocating memory for the variable and values arrays */
913  	   nmastervars = SCIPgetNVars(masterprob) + SCIPgetNFixedVars(masterprob);
914  	   SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vars, nmastervars) );
915  	   SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vals, nmastervars) );
916  	   lhs = 0.0;
917  	   rhs = SCIPinfinity(masterprob);
918  	   nvars = 0;
919  	   varssize = nmastervars;
920  	
921  	   if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) )
922  	   {
923  	      /* computing the coefficients of the optimality cut */
924  	      SCIP_CALL( computeStandardNLPOptimalityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &rhs, &nvars,
925  	            &varssize, objective, primalvals, consdualvals, varlbdualvals, varubdualvals, row2idx,
926  	            var2idx, &checkobj, &success) );
927  	   }
928  	   else
929  	   {
930  	      /* computing the coefficients of the optimality cut */
931  	      SCIP_CALL( computeStandardLPOptimalityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &rhs, &nvars,
932  	            &varssize, &checkobj, &success) );
933  	   }
934  	
935  	   /* if success is FALSE, then there was an error in generating the optimality cut. No cut will be added to the master
936  	    * problem. Otherwise, the constraint is added to the master problem.
937  	    */
938  	   if( !success )
939  	   {
940  	      (*result) = SCIP_DIDNOTFIND;
941  	      SCIPdebugMsg(masterprob, "Error in generating Benders' optimality cut for problem %d.\n", probnumber);
942  	   }
943  	   else
944  	   {
945  	      /* initially a row/constraint is created for the optimality cut using the master variables and coefficients
946  	       * computed in computeStandardLPOptimalityCut. At this stage, the auxiliary variable is not added since the
947  	       * activity of the row/constraint in its current form is used to determine the validity of the optimality cut.
948  	       */
949  	      if( addcut )
950  	      {
951  	         SCIP_CALL( SCIPcreateEmptyRowConshdlr(masterprob, &row, consbenders, cutname, lhs, rhs, FALSE, FALSE, TRUE) );
952  	         SCIP_CALL( SCIPaddVarsToRow(masterprob, row, nvars, vars, vals) );
953  	      }
954  	      else
955  	      {
956  	         SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cons, cutname, nvars, vars, vals, lhs, rhs) );
957  	         SCIP_CALL( SCIPsetConsDynamic(masterprob, cons, TRUE) );
958  	         SCIP_CALL( SCIPsetConsRemovable(masterprob, cons, TRUE) );
959  	      }
960  	
961  	      /* computing the objective function from the cut activity to verify the accuracy of the constraint */
962  	      verifyobj = 0.0;
963  	      if( addcut )
964  	      {
965  	         verifyobj += SCIProwGetLhs(row) - SCIPgetRowSolActivity(masterprob, row, sol);
966  	      }
967  	      else
968  	      {
969  	         verifyobj += SCIPgetLhsLinear(masterprob, cons) - SCIPgetActivityLinear(masterprob, cons, sol);
970  	      }
971  	
972  	      if( feasibilitycut && verifyobj < SCIPfeastol(masterprob) )
973  	      {
974  	         success = FALSE;
975  	         SCIPdebugMsg(masterprob, "The violation of the feasibility cut (%g) is too small. Skipping feasibility cut.\n", verifyobj);
976  	      }
977  	
978  	      /* it is possible that numerics will cause the generated cut to be invalid. This cut should not be added to the
979  	       * master problem, since its addition could cut off feasible solutions. The success flag is set of false, indicating
980  	       * that the Benders' cut could not find a valid cut.
981  	       */
982  	      if( !feasibilitycut && !SCIPisFeasEQ(masterprob, checkobj, verifyobj) )
983  	      {
984  	         SCIP_Bool valid;
985  	
986  	         /* the difference in the checkobj and verifyobj could be due to the setup tolerances. This is checked, and if
987  	          * so, then the generated cut is still valid
988  	          */
989  	         SCIP_CALL( checkSetupTolerances(masterprob, sol, vars, vals, lhs, checkobj, nvars, &valid) );
990  	
991  	         if( !valid )
992  	         {
993  	            success = FALSE;
994  	            SCIPdebugMsg(masterprob, "The objective function and cut activity are not equal (%g != %g).\n", checkobj,
995  	               verifyobj);
996  	
997  	#ifdef SCIP_DEBUG
998  	            /* we only need to abort if cut strengthen is not used. If cut strengthen has been used in this round and the
999  	             * cut could not be generated, then another subproblem solving round will be executed
1000 	             */
1001 	            if( !SCIPbendersInStrengthenRound(benders) )
1002 	            {
1003 	#ifdef SCIP_MOREDEBUG
1004 	               int i;
1005 	
1006 	               for( i = 0; i < nvars; i++ )
1007 	                  printf("<%s> %g %g\n", SCIPvarGetName(vars[i]), vals[i], SCIPgetSolVal(masterprob, sol, vars[i]));
1008 	#endif
1009 	               SCIPABORT();
1010 	            }
1011 	#endif
1012 	         }
1013 	      }
1014 	
1015 	      if( success )
1016 	      {
1017 	         /* adding the auxiliary variable to the optimality cut. The auxiliary variable is added to the vars and vals
1018 	          * arrays prior to the execution of the MIR procedure. This is necessary because the MIR procedure must be
1019 	          * executed on the complete cut, not just the row/constraint without the auxiliary variable.
1020 	          */
1021 	         if( !feasibilitycut )
1022 	         {
1023 	            SCIP_CALL( addAuxiliaryVariableToCut(masterprob, benders, vars, vals, &nvars, probnumber) );
1024 	         }
1025 	
1026 	         /* performing the MIR procedure. If the procedure is successful, then the vars and vals arrays are no longer
1027 	          * needed for creating the optimality cut. These are superseeded with the cutcoefs and cutinds arrays. In the
1028 	          * case that the MIR procedure is successful, the row/constraint that has been created previously is destroyed
1029 	          * and the MIR cut is added in its place
1030 	          */
1031 	         if( calcmir )
1032 	         {
1033 	            SCIP_Real* cutcoefs;
1034 	            int* cutinds;
1035 	            SCIP_Real cutrhs;
1036 	            int cutnnz;
1037 	
1038 	            /* allocating memory to compute the MIR cut */
1039 	            SCIP_CALL( SCIPallocBufferArray(masterprob, &cutcoefs, nvars) );
1040 	            SCIP_CALL( SCIPallocBufferArray(masterprob, &cutinds, nvars) );
1041 	
1042 	            SCIP_CALL( computeMIRForOptimalityCut(masterprob, sol, vars, vals, lhs, rhs, nvars, cutcoefs,
1043 	                  cutinds, &cutrhs, &cutnnz, &mirsuccess) );
1044 	
1045 	            /* if the MIR cut was computed successfully, then the current row/constraint needs to be destroyed and
1046 	             * replaced with the updated coefficients
1047 	             */
1048 	            if( mirsuccess )
1049 	            {
1050 	               SCIP_VAR** mastervars;
1051 	               int i;
1052 	
1053 	               mastervars = SCIPgetVars(masterprob);
1054 	
1055 	               if( addcut )
1056 	               {
1057 	                  SCIP_CALL( SCIPreleaseRow(masterprob, &row) );
1058 	
1059 	                  SCIP_CALL( SCIPcreateEmptyRowConshdlr(masterprob, &row, consbenders, cutname,
1060 	                        -SCIPinfinity(masterprob), cutrhs, FALSE, FALSE, TRUE) );
1061 	
1062 	                  for( i = 0; i < cutnnz; i++)
1063 	                  {
1064 	                     SCIP_CALL( SCIPaddVarToRow(masterprob, row, mastervars[cutinds[i]], cutcoefs[i]) );
1065 	                  }
1066 	               }
1067 	               else
1068 	               {
1069 	                  SCIP_CALL( SCIPreleaseCons(masterprob, &cons) );
1070 	
1071 	                  SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cons, cutname, 0, NULL, NULL,
1072 	                        -SCIPinfinity(masterprob), cutrhs) );
1073 	                  SCIP_CALL( SCIPsetConsDynamic(masterprob, cons, TRUE) );
1074 	                  SCIP_CALL( SCIPsetConsRemovable(masterprob, cons, TRUE) );
1075 	
1076 	                  for( i = 0; i < cutnnz; i++ )
1077 	                  {
1078 	                     SCIP_CALL( SCIPaddCoefLinear(masterprob, cons, mastervars[cutinds[i]], cutcoefs[i]) );
1079 	                  }
1080 	               }
1081 	            }
1082 	
1083 	            /* freeing the memory required to compute the MIR cut */
1084 	            SCIPfreeBufferArray(masterprob, &cutinds);
1085 	            SCIPfreeBufferArray(masterprob, &cutcoefs);
1086 	         }
1087 	
1088 	         /* adding the constraint to the master problem */
1089 	         if( addcut )
1090 	         {
1091 	            SCIP_Bool infeasible;
1092 	
1093 	            /* adding the auxiliary variable coefficient to the row. This is only added if the MIR procedure is not
1094 	             * successful. If the MIR procedure was successful, then the auxiliary variable is already included in the
1095 	             * row
1096 	             */
1097 	            if( !feasibilitycut && !mirsuccess )
1098 	            {
1099 	               SCIP_CALL( SCIPaddVarToRow(masterprob, row, vars[nvars - 1], vals[nvars - 1]) );
1100 	            }
1101 	
1102 	            if( type == SCIP_BENDERSENFOTYPE_LP || type == SCIP_BENDERSENFOTYPE_RELAX )
1103 	            {
1104 	               SCIP_CALL( SCIPaddRow(masterprob, row, FALSE, &infeasible) );
1105 	               assert(!infeasible);
1106 	            }
1107 	            else
1108 	            {
1109 	               assert(type == SCIP_BENDERSENFOTYPE_CHECK || type == SCIP_BENDERSENFOTYPE_PSEUDO);
1110 	               SCIP_CALL( SCIPaddPoolCut(masterprob, row) );
1111 	            }
1112 	
1113 	            (*result) = SCIP_SEPARATED;
1114 	         }
1115 	         else
1116 	         {
1117 	            /* adding the auxiliary variable coefficient to the row. This is only added if the MIR procedure is not
1118 	             * successful. If the MIR procedure was successful, then the auxiliary variable is already included in the
1119 	             * constraint.
1120 	             */
1121 	            if( !feasibilitycut && !mirsuccess )
1122 	            {
1123 	               SCIP_CALL( SCIPaddCoefLinear(masterprob, cons, vars[nvars - 1], vals[nvars - 1]) );
1124 	            }
1125 	
1126 	            SCIPdebugPrintCons(masterprob, cons, NULL);
1127 	
1128 	            SCIP_CALL( SCIPaddCons(masterprob, cons) );
1129 	
1130 	            (*result) = SCIP_CONSADDED;
1131 	         }
1132 	
1133 	         /* storing the data that is used to create the cut */
1134 	         SCIP_CALL( SCIPstoreBendersCut(masterprob, benders, vars, vals, lhs, rhs, nvars) );
1135 	      }
1136 	      else
1137 	      {
1138 	         (*result) = SCIP_DIDNOTFIND;
1139 	         SCIPdebugMsg(masterprob, "Error in generating Benders' %s cut for problem %d.\n", feasibilitycut ? "feasibility" : "optimality", probnumber);
1140 	      }
1141 	
1142 	      /* releasing the row or constraint */
1143 	      if( addcut )
1144 	      {
1145 	         /* release the row */
1146 	         SCIP_CALL( SCIPreleaseRow(masterprob, &row) );
1147 	      }
1148 	      else
1149 	      {
1150 	         /* release the constraint */
1151 	         SCIP_CALL( SCIPreleaseCons(masterprob, &cons) );
1152 	      }
1153 	   }
1154 	
1155 	   SCIPfreeBufferArray(masterprob, &vals);
1156 	   SCIPfreeBufferArray(masterprob, &vars);
1157 	
1158 	   return SCIP_OKAY;
1159 	}
1160 	
1161 	
1162 	/** adds the gradient of a nonlinear row in the current NLP solution of a subproblem to a linear row or constraint in the master problem
1163 	 *
1164 	 * Only computes gradient w.r.t. master problem variables.
1165 	 * Computes also the directional derivative, that is, mult times gradient times solution.
1166 	 */
1167 	SCIP_RETCODE SCIPaddNlRowGradientBenderscutOpt(
1168 	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
1169 	   SCIP*                 subproblem,         /**< the SCIP instance of the subproblem */
1170 	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
1171 	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
1172 	   SCIP_Real             mult,               /**< multiplier */
1173 	   SCIP_Real*            primalvals,         /**< the primal solutions for the NLP, can be NULL */
1174 	   SCIP_HASHMAP*         var2idx,            /**< mapping from variable of the subproblem to the index in the dual arrays, can be NULL */
1175 	   SCIP_Real*            dirderiv,           /**< storage to add directional derivative */
1176 	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
1177 	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
1178 	   int*                  nvars,              /**< the number of variables in the cut */
1179 	   int*                  varssize            /**< the number of variables in the array */
1180 	   )
1181 	{
1182 	   SCIP_EXPR* expr;
1183 	   SCIP_VAR* var;
1184 	   SCIP_VAR* mastervar;
1185 	   SCIP_Real coef;
1186 	   int i;
1187 	
1188 	   assert(masterprob != NULL);
1189 	   assert(subproblem != NULL);
1190 	   assert(benders != NULL);
1191 	   assert(nlrow != NULL);
1192 	   assert((primalvals == NULL && var2idx == NULL) || (primalvals != NULL && var2idx != NULL));
1193 	   assert(mult != 0.0);
1194 	   assert(dirderiv != NULL);
1195 	   assert(vars != NULL);
1196 	   assert(vals != NULL);
1197 	
1198 	   /* linear part */
1199 	   for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ )
1200 	   {
1201 	      var = SCIPnlrowGetLinearVars(nlrow)[i];
1202 	      assert(var != NULL);
1203 	
1204 	      /* retrieving the master problem variable for the given subproblem variable. */
1205 	      SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) );
1206 	      if( mastervar == NULL )
1207 	         continue;
1208 	
1209 	      coef = mult * SCIPnlrowGetLinearCoefs(nlrow)[i];
1210 	
1211 	      /* adding the variable to the storage */
1212 	      SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) );
1213 	
1214 	      *dirderiv += coef * getNlpVarSol(var, primalvals, var2idx);
1215 	   }
1216 	
1217 	   /* expression part */
1218 	   expr = SCIPnlrowGetExpr(nlrow);
1219 	   if( expr != NULL )
1220 	   {
1221 	      SCIP_SOL* primalsol;
1222 	      SCIP_EXPRITER* it;
1223 	
1224 	      /* create primalsol, either from primalvals, or pointing to NLP solution */
1225 	      if( primalvals != NULL )
1226 	      {
1227 	         SCIP_CALL( SCIPcreateSol(subproblem, &primalsol, NULL) );
1228 	
1229 	         /* TODO would be better to change primalvals to a SCIP_SOL and do this once for the whole NLP instead of repeating it for each expr */
1230 	         for( i = 0; i < SCIPhashmapGetNEntries(var2idx); ++i )
1231 	         {
1232 	            SCIP_HASHMAPENTRY* entry;
1233 	            entry = SCIPhashmapGetEntry(var2idx, i);
1234 	            if( entry == NULL )
1235 	               continue;
1236 	            SCIP_CALL( SCIPsetSolVal(subproblem, primalsol, (SCIP_VAR*) SCIPhashmapEntryGetOrigin(entry), primalvals[SCIPhashmapEntryGetImageInt(entry)]) );
1237 	         }
1238 	      }
1239 	      else
1240 	      {
1241 	         SCIP_CALL( SCIPcreateNLPSol(subproblem, &primalsol, NULL) );
1242 	      }
1243 	
1244 	      /* eval gradient */
1245 	      SCIP_CALL( SCIPevalExprGradient(subproblem, expr, primalsol, 0L) );
1246 	
1247 	      assert(SCIPexprGetDerivative(expr) != SCIP_INVALID);  /* TODO this should be a proper check&abort */ /*lint !e777*/
1248 	
1249 	      SCIP_CALL( SCIPfreeSol(subproblem, &primalsol) );
1250 	
1251 	      /* update corresponding gradient entry */
1252 	      SCIP_CALL( SCIPcreateExpriter(subproblem, &it) );
1253 	      SCIP_CALL( SCIPexpriterInit(it, expr, SCIP_EXPRITER_DFS, FALSE) );
1254 	      for( ; !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )  /*lint !e441*/ /*lint !e440*/
1255 	      {
1256 	         if( !SCIPisExprVar(subproblem, expr) )
1257 	            continue;
1258 	
1259 	         var = SCIPgetVarExprVar(expr);
1260 	         assert(var != NULL);
1261 	
1262 	         /* retrieving the master problem variable for the given subproblem variable. */
1263 	         SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) );
1264 	         if( mastervar == NULL )
1265 	            continue;
1266 	
1267 	         assert(SCIPexprGetDerivative(expr) != SCIP_INVALID);  /*lint !e777*/
1268 	         coef = mult * SCIPexprGetDerivative(expr);
1269 	
1270 	         /* adding the variable to the storage */
1271 	         SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, coef, nvars, varssize) );
1272 	
1273 	         *dirderiv += coef * getNlpVarSol(var, primalvals, var2idx);
1274 	      }
1275 	      SCIPfreeExpriter(&it);
1276 	   }
1277 	
1278 	   return SCIP_OKAY;
1279 	}
1280