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_feas.c
26   	 * @ingroup OTHER_CFILES
27   	 * @brief  Standard feasibility cuts 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 "scip/pub_expr.h"
34   	#include "scip/benderscut_feas.h"
35   	#include "scip/benderscut_opt.h"
36   	#include "scip/cons_linear.h"
37   	#include "scip/pub_benderscut.h"
38   	#include "scip/pub_benders.h"
39   	#include "scip/pub_lp.h"
40   	#include "scip/pub_message.h"
41   	#include "scip/pub_misc.h"
42   	#include "scip/pub_misc_linear.h"
43   	#include "scip/pub_nlp.h"
44   	#include "scip/pub_var.h"
45   	#include "scip/scip_benders.h"
46   	#include "scip/scip_cons.h"
47   	#include "scip/scip_general.h"
48   	#include "scip/scip_lp.h"
49   	#include "scip/scip_mem.h"
50   	#include "scip/scip_message.h"
51   	#include "scip/scip_nlp.h"
52   	#include "scip/scip_nlpi.h"
53   	#include "scip/scip_numerics.h"
54   	#include "scip/scip_prob.h"
55   	#include "scip/scip_solvingstats.h"
56   	#include "scip/scip_var.h"
57   	
58   	#define BENDERSCUT_NAME             "feas"
59   	#define BENDERSCUT_DESC             "Standard feasibility cuts for Benders' decomposition"
60   	#define BENDERSCUT_PRIORITY     10000
61   	#define BENDERSCUT_LPCUT         TRUE
62   	
63   	/*
64   	 * Local methods
65   	 */
66   	
67   	/** adds a variable and value to the constraint/row arrays */
68   	static
69   	SCIP_RETCODE addVariableToArray(
70   	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
71   	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
72   	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
73   	   SCIP_VAR*             addvar,             /**< the variable that will be added to the array */
74   	   SCIP_Real             addval,             /**< the value that will be added to the array */
75   	   int*                  nvars,              /**< the number of variables in the variable array */
76   	   int*                  varssize            /**< the length of the variable size */
77   	   )
78   	{
79   	   assert(masterprob != NULL);
80   	   assert(vars != NULL);
81   	   assert(*vars != NULL);
82   	   assert(vals != NULL);
83   	   assert(*vals != NULL);
84   	   assert(addvar != NULL);
85   	   assert(nvars != NULL);
86   	   assert(varssize != NULL);
87   	
88   	   if( *nvars >= *varssize )
89   	   {
90   	      *varssize = SCIPcalcMemGrowSize(masterprob, *varssize + 1);
91   	      SCIP_CALL( SCIPreallocBufferArray(masterprob, vars, *varssize) );
92   	      SCIP_CALL( SCIPreallocBufferArray(masterprob, vals, *varssize) );
93   	   }
94   	   assert(*nvars < *varssize);
95   	
96   	   (*vars)[*nvars] = addvar;
97   	   (*vals)[*nvars] = addval;
98   	   (*nvars)++;
99   	
100  	   return SCIP_OKAY;
101  	}
102  	
103  	/** computing as standard Benders' feasibility cut from the dual solutions of the LP */
104  	static
105  	SCIP_RETCODE computeStandardLPFeasibilityCut(
106  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
107  	   SCIP*                 subproblem,         /**< the SCIP instance of the pricing problem */
108  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
109  	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
110  	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
111  	   SCIP_Real*            lhs,                /**< the left hand side of the cut */
112  	   int*                  nvars,              /**< the number of variables in the cut */
113  	   int*                  varssize,           /**< the number of variables in the array */
114  	   SCIP_Bool*            success             /**< was the cut generation successful? */
115  	   )
116  	{
117  	   SCIP_VAR** subvars;
118  	   int nsubvars;
119  	   int nrows;
120  	   SCIP_Real dualsol;
121  	   SCIP_Real addval;    /* the value that must be added to the lhs */
122  	   int i;
123  	
124  	   assert(masterprob != NULL);
125  	   assert(subproblem != NULL);
126  	   assert(benders != NULL);
127  	   assert(SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_INFEASIBLE);
128  	
129  	   (*success) = FALSE;
130  	
131  	   /* looping over all LP rows and setting the coefficients of the cut */
132  	   nrows = SCIPgetNLPRows(subproblem);
133  	   for( i = 0; i < nrows; i++ )
134  	   {
135  	      SCIP_ROW* lprow;
136  	
137  	      lprow = SCIPgetLPRows(subproblem)[i];
138  	      assert(lprow != NULL);
139  	
140  	      dualsol = SCIProwGetDualfarkas(lprow);
141  	      assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) );
142  	
143  	      if( SCIPisDualfeasZero(subproblem, dualsol) )
144  	         continue;
145  	
146  	      if( dualsol > 0.0 )
147  	         addval = dualsol*SCIProwGetLhs(lprow);
148  	      else
149  	         addval = dualsol*SCIProwGetRhs(lprow);
150  	
151  	      *lhs += addval;
152  	
153  	      /* if the bound becomes infinite, then the cut generation terminates. */
154  	      if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs)
155  	         || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval))
156  	      {
157  	         (*success) = FALSE;
158  	         SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut.\n");
159  	         return SCIP_OKAY;
160  	      }
161  	   }
162  	
163  	   nsubvars = SCIPgetNVars(subproblem);
164  	   subvars = SCIPgetVars(subproblem);
165  	
166  	   /* looping over all variables to update the coefficients in the computed cut. */
167  	   for( i = 0; i < nsubvars; i++ )
168  	   {
169  	      SCIP_VAR* var;
170  	      SCIP_VAR* mastervar;
171  	
172  	      var = subvars[i];
173  	
174  	      /* retrieving the master problem variable for the given subproblem variable. */
175  	      SCIP_CALL( SCIPgetBendersMasterVar(masterprob, benders, var, &mastervar) );
176  	
177  	      dualsol = SCIPgetVarFarkasCoef(subproblem, var);
178  	
179  	      if( SCIPisZero(subproblem, dualsol) )
180  	         continue;
181  	
182  	      /* checking whether the original variable is a linking variable.
183  	       * If this is the case, then the corresponding master variable is added to the generated cut.
184  	       * If the pricing variable is not a linking variable, then the farkas dual value is added to the lhs
185  	       */
186  	      if( mastervar != NULL )
187  	      {
188  	         SCIPdebugMsg(masterprob ,"Adding coeffs to feasibility cut: <%s> dualsol %g\n", SCIPvarGetName(mastervar), dualsol);
189  	
190  	         /* adding the variable to the storage */
191  	         SCIP_CALL( addVariableToArray(masterprob, vars, vals, mastervar, dualsol, nvars, varssize) );
192  	      }
193  	      else
194  	      {
195  	         addval = 0;
196  	
197  	         if( SCIPisPositive(subproblem, dualsol) )
198  	            addval = dualsol*SCIPvarGetUbGlobal(var);
199  	         else if( SCIPisNegative(subproblem, dualsol) )
200  	            addval = dualsol*SCIPvarGetLbGlobal(var);
201  	
202  	         *lhs -= addval;
203  	
204  	         /* if the bound becomes infinite, then the cut generation terminates. */
205  	         if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs)
206  	            || SCIPisInfinity(masterprob, addval) || SCIPisInfinity(masterprob, -addval))
207  	         {
208  	            (*success) = FALSE;
209  	            SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut.\n");
210  	            return SCIP_OKAY;
211  	         }
212  	      }
213  	   }
214  	
215  	   (*success) = TRUE;
216  	
217  	   return SCIP_OKAY;
218  	}
219  	
220  	
221  	/** computing as standard Benders' feasibility cut from the dual solutions of the NLP
222  	 *
223  	 *  NOTE: The cut must be created before being passed to this function
224  	 */
225  	static
226  	SCIP_RETCODE computeStandardNLPFeasibilityCut(
227  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
228  	   SCIP*                 subproblem,         /**< the SCIP instance of the pricing problem */
229  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition structure */
230  	   SCIP_VAR***           vars,               /**< pointer to array of variables in the generated cut with non-zero coefficient */
231  	   SCIP_Real**           vals,               /**< pointer to array of coefficients of the variables in the generated cut */
232  	   SCIP_Real*            lhs,                /**< the left hand side of the cut */
233  	   int*                  nvars,              /**< the number of variables in the cut */
234  	   int*                  varssize,           /**< the number of variables in the array */
235  	   SCIP_Bool*            success             /**< was the cut generation successful? */
236  	   )
237  	{
238  	   int nrows;
239  	   SCIP_Real activity;
240  	   SCIP_Real dirderiv;
241  	   SCIP_Real dualsol;
242  	   int i;
243  	
244  	   assert(masterprob != NULL);
245  	   assert(subproblem != NULL);
246  	   assert(benders != NULL);
247  	   assert(SCIPisNLPConstructed(subproblem));
248  	   assert(SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_LOCINFEASIBLE || SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_GLOBINFEASIBLE);
249  	
250  	   (*success) = FALSE;
251  	
252  	   *lhs = 0.0;
253  	   dirderiv = 0.0;
254  	
255  	   /* looping over all NLP rows and setting the corresponding coefficients of the cut */
256  	   nrows = SCIPgetNNLPNlRows(subproblem);
257  	   for( i = 0; i < nrows; i++ )
258  	   {
259  	      SCIP_NLROW* nlrow;
260  	
261  	      nlrow = SCIPgetNLPNlRows(subproblem)[i];
262  	      assert(nlrow != NULL);
263  	
264  	      dualsol = SCIPnlrowGetDualsol(nlrow);
265  	      assert( !SCIPisInfinity(subproblem, dualsol) && !SCIPisInfinity(subproblem, -dualsol) );
266  	
267  	      if( SCIPisZero(subproblem, dualsol) )
268  	         continue;
269  	
270  	      SCIP_CALL( SCIPaddNlRowGradientBenderscutOpt(masterprob, subproblem, benders, nlrow, -dualsol,
271  	            NULL, NULL, &dirderiv, vars, vals, nvars, varssize) );
272  	
273  	      SCIP_CALL( SCIPgetNlRowActivity(subproblem, nlrow, &activity) );
274  	
275  	      if( dualsol > 0.0 )
276  	      {
277  	         assert(!SCIPisInfinity(subproblem, SCIPnlrowGetRhs(nlrow)));
278  	         *lhs += dualsol * (activity - SCIPnlrowGetRhs(nlrow));
279  	      }
280  	      else
281  	      {
282  	         assert(!SCIPisInfinity(subproblem, -SCIPnlrowGetLhs(nlrow)));
283  	         *lhs += dualsol * (activity - SCIPnlrowGetLhs(nlrow));
284  	      }
285  	   }
286  	
287  	   *lhs += dirderiv;
288  	
289  	   /* if the side became infinite or dirderiv was infinite, then the cut generation terminates. */
290  	   if( SCIPisInfinity(masterprob, *lhs) || SCIPisInfinity(masterprob, -*lhs)
291  	      || SCIPisInfinity(masterprob, dirderiv) || SCIPisInfinity(masterprob, -dirderiv))
292  	   {
293  	      (*success) = FALSE;
294  	      SCIPdebugMsg(masterprob, "Infinite bound when generating feasibility cut. lhs = %g dirderiv = %g.\n", *lhs, dirderiv);
295  	      return SCIP_OKAY;
296  	   }
297  	
298  	   (*success) = TRUE;
299  	
300  	   return SCIP_OKAY;
301  	}
302  	
303  	/** generates and applies Benders' cuts */
304  	static
305  	SCIP_RETCODE generateAndApplyBendersCuts(
306  	   SCIP*                 masterprob,         /**< the SCIP instance of the master problem */
307  	   SCIP*                 subproblem,         /**< the SCIP instance of the pricing problem */
308  	   SCIP_BENDERS*         benders,            /**< the benders' decomposition */
309  	   SCIP_BENDERSCUT*      benderscut,         /**< the benders' decomposition cut method */
310  	   SCIP_SOL*             sol,                /**< primal CIP solution */
311  	   int                   probnumber,         /**< the number of the pricing problem */
312  	   SCIP_RESULT*          result              /**< the result from solving the subproblems */
313  	   )
314  	{
315  	   SCIP_CONS* cut;
316  	   SCIP_VAR** vars;
317  	   SCIP_Real* vals;
318  	   SCIP_Real lhs;
319  	   SCIP_Real activity;
320  	   int nvars;
321  	   int varssize;
322  	   int nmastervars;
323  	   char cutname[SCIP_MAXSTRLEN];
324  	   SCIP_Bool success;
325  	
326  	   assert(masterprob != NULL);
327  	   assert(subproblem != NULL);
328  	   assert(benders != NULL);
329  	   assert(result != NULL);
330  	
331  	   /* allocating memory for the variable and values arrays */
332  	   nmastervars = SCIPgetNVars(masterprob) + SCIPgetNFixedVars(masterprob);
333  	   SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vars, nmastervars) );
334  	   SCIP_CALL( SCIPallocClearBufferArray(masterprob, &vals, nmastervars) );
335  	   lhs = 0.0;
336  	   nvars = 0;
337  	   varssize = nmastervars;
338  	
339  	   /* setting the name of the generated cut */
340  	   (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "feasibilitycut_%d_%" SCIP_LONGINT_FORMAT, probnumber,
341  	      SCIPbenderscutGetNFound(benderscut) );
342  	
343  	   if( SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem) )
344  	   {
345  	      /* computing the coefficients of the feasibility cut from the NLP */
346  	      SCIP_CALL( computeStandardNLPFeasibilityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &nvars, &varssize,
347  	            &success) );
348  	   }
349  	   else
350  	   {
351  	      if( SCIPgetNLPIterations(subproblem) == 0 )
352  	      {
353  	         SCIPverbMessage(masterprob, SCIP_VERBLEVEL_FULL, NULL, "There were no iterations in pricing problem %d. "
354  	           "A Benders' decomposition feasibility cut will be generated from the presolved LP data.\n", probnumber);
355  	      }
356  	
357  	      /* computing the coefficients of the feasibility cut from the LP */
358  	      SCIP_CALL( computeStandardLPFeasibilityCut(masterprob, subproblem, benders, &vars, &vals, &lhs, &nvars, &varssize,
359  	            &success) );
360  	   }
361  	
362  	   /* if success is FALSE, then there was an error in generating the feasibility cut. No cut will be added to the master
363  	    * problem. Otherwise, the constraint is added to the master problem.
364  	    */
365  	   if( !success )
366  	   {
367  	      (*result) = SCIP_DIDNOTFIND;
368  	      SCIPdebugMsg(masterprob, "Error in generating Benders' feasibility cut for problem %d.\n", probnumber);
369  	   }
370  	   else
371  	   {
372  	      /* creating a constraint with the variables and coefficients previously stored */
373  	      SCIP_CALL( SCIPcreateConsBasicLinear(masterprob, &cut, cutname, nvars, vars, vals, lhs, SCIPinfinity(masterprob)) );
374  	      SCIP_CALL( SCIPsetConsDynamic(masterprob, cut, TRUE) );
375  	      SCIP_CALL( SCIPsetConsRemovable(masterprob, cut, TRUE) );
376  	
377  	      assert(SCIPisInfinity(masterprob, SCIPgetRhsLinear(masterprob, cut)));
378  	
379  	      /* the activity of the cut should be less than the lhs. This will ensure that the evaluated solution will be cut off.
380  	       * It is possible that the activity is greater than the lhs. This could be caused by numerical difficulties. In this
381  	       * case, no cut will be generated.
382  	       */
383  	      lhs = SCIPgetLhsLinear(masterprob, cut);
384  	      activity = SCIPgetActivityLinear(masterprob, cut, sol);
385  	      if( SCIPisGE(masterprob, activity, lhs) )
386  	      {
387  	         success = FALSE;
388  	         SCIPdebugMsg(masterprob ,"Invalid feasibility cut - activity is greater than lhs %g >= %g.\n", activity, lhs);
389  	#ifdef SCIP_DEBUG
390  	         SCIPABORT();
391  	#endif
392  	      }
393  	
394  	      assert(cut != NULL);
395  	
396  	      if( success )
397  	      {
398  	         /* adding the constraint to the master problem */
399  	         SCIP_CALL( SCIPaddCons(masterprob, cut) );
400  	
401  	         SCIPdebugPrintCons(masterprob, cut, NULL);
402  	
403  	         (*result) = SCIP_CONSADDED;
404  	      }
405  	
406  	      SCIP_CALL( SCIPreleaseCons(masterprob, &cut) );
407  	   }
408  	
409  	   SCIPfreeBufferArray(masterprob, &vals);
410  	   SCIPfreeBufferArray(masterprob, &vars);
411  	
412  	   return SCIP_OKAY;
413  	}
414  	
415  	/*
416  	 * Callback methods of Benders' decomposition cuts
417  	 */
418  	
419  	/** execution method of Benders' decomposition cuts */
420  	static
421  	SCIP_DECL_BENDERSCUTEXEC(benderscutExecFeas)
422  	{  /*lint --e{715}*/
423  	   SCIP* subproblem;
424  	   SCIP_Bool nlprelaxation;
425  	
426  	   assert(scip != NULL);
427  	   assert(benders != NULL);
428  	   assert(benderscut != NULL);
429  	   assert(result != NULL);
430  	   assert(probnumber >= 0 && probnumber < SCIPbendersGetNSubproblems(benders));
431  	
432  	   subproblem = SCIPbendersSubproblem(benders, probnumber);
433  	
434  	   if( subproblem == NULL )
435  	   {
436  	      SCIPdebugMsg(scip, "The subproblem %d is set to NULL. The <%s> Benders' decomposition cut can not be executed.\n",
437  	         probnumber, BENDERSCUT_NAME);
438  	
439  	      (*result) = SCIP_DIDNOTRUN;
440  	      return SCIP_OKAY;
441  	   }
442  	
443  	   /* setting a flag to indicate whether the NLP relaxation should be used to generate cuts */
444  	   nlprelaxation = SCIPisNLPConstructed(subproblem) && SCIPgetNNlpis(subproblem);
445  	
446  	   /* only generate feasibility cuts if the subproblem LP or NLP is infeasible,
447  	    * since we use the farkas proof from the LP or the dual solution of the NLP to construct the feasibility cut
448  	    */
449  	   if( SCIPgetStage(subproblem) == SCIP_STAGE_SOLVING &&
450  	      ((!nlprelaxation && SCIPgetLPSolstat(subproblem) == SCIP_LPSOLSTAT_INFEASIBLE) ||
451  	       (nlprelaxation && (SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_LOCINFEASIBLE || SCIPgetNLPSolstat(subproblem) == SCIP_NLPSOLSTAT_GLOBINFEASIBLE))) )
452  	   {
453  	      /* generating a cut for a given subproblem */
454  	      SCIP_CALL( generateAndApplyBendersCuts(scip, subproblem, benders, benderscut,
455  	            sol, probnumber, result) );
456  	   }
457  	
458  	   return SCIP_OKAY;
459  	}
460  	
461  	
462  	/*
463  	 * Benders' decomposition cuts specific interface methods
464  	 */
465  	
466  	/** creates the Standard Feasibility Benders' decomposition cuts and includes it in SCIP */
467  	SCIP_RETCODE SCIPincludeBenderscutFeas(
468  	   SCIP*                 scip,               /**< SCIP data structure */
469  	   SCIP_BENDERS*         benders             /**< Benders' decomposition */
470  	   )
471  	{
472  	   SCIP_BENDERSCUT* benderscut;
473  	
474  	   assert(benders != NULL);
475  	
476  	   benderscut = NULL;
477  	
478  	   /* include Benders' decomposition cuts */
479  	   SCIP_CALL( SCIPincludeBenderscutBasic(scip, benders, &benderscut, BENDERSCUT_NAME, BENDERSCUT_DESC,
480  	         BENDERSCUT_PRIORITY, BENDERSCUT_LPCUT, benderscutExecFeas, NULL) );
481  	
482  	   assert(benderscut != NULL);
483  	
484  	   return SCIP_OKAY;
485  	}
486