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   heur_mpec.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  mpec primal heuristic
28   	 * @author Felipe Serrano
29   	 * @author Benjamin Mueller
30   	 */
31   	
32   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
33   	
34   	#include "blockmemshell/memory.h"
35   	#include "scip/pub_expr.h"
36   	#include "scip/expr_var.h"
37   	#include "scip/expr_sum.h"
38   	#include "scip/expr_pow.h"
39   	#include "scip/heur_mpec.h"
40   	#include "scip/heur_subnlp.h"
41   	#include "scip/pub_cons.h"
42   	#include "scip/pub_heur.h"
43   	#include "scip/pub_message.h"
44   	#include "scip/pub_misc.h"
45   	#include "scip/pub_nlp.h"
46   	#include "scip/pub_var.h"
47   	#include "scip/scip_cons.h"
48   	#include "scip/scip_general.h"
49   	#include "scip/scip_heur.h"
50   	#include "scip/scip_mem.h"
51   	#include "scip/scip_message.h"
52   	#include "scip/scip_nlp.h"
53   	#include "scip/scip_nlpi.h"
54   	#include "scip/scip_numerics.h"
55   	#include "scip/scip_param.h"
56   	#include "scip/scip_prob.h"
57   	#include "scip/scip_sol.h"
58   	#include "scip/scip_solvingstats.h"
59   	#include "scip/scip_timing.h"
60   	#include <string.h>
61   	
62   	
63   	#define HEUR_NAME             "mpec"
64   	#define HEUR_DESC             "regularization heuristic for convex and nonconvex MINLPs"
65   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_DIVING
66   	#define HEUR_PRIORITY         -2050000
67   	#define HEUR_FREQ             50
68   	#define HEUR_FREQOFS          0
69   	#define HEUR_MAXDEPTH         -1
70   	#define HEUR_TIMING           SCIP_HEURTIMING_AFTERLPNODE
71   	#define HEUR_USESSUBSCIP      TRUE           /**< disable the heuristic in sub-SCIPs, even though it does not use any */
72   	
73   	#define DEFAULT_INITTHETA     0.125          /**< default initial regularization right-hand side value (< 0.25) */
74   	#define DEFAULT_SIGMA         0.5            /**< default regularization update factor (< 1) */
75   	#define DEFAULT_MAXITER       100            /**< default maximum number of iterations of the MPEC loop */
76   	#define DEFAULT_MAXNLPITER    500            /**< default maximum number of NLP iterations per solve */
77   	#define DEFAULT_MINGAPLEFT    0.05           /**< default minimum amount of gap left in order to call the heuristic */
78   	#define DEFAULT_SUBNLPTRIGGER 1e-3           /**< default maximum integrality violation before triggering a sub-NLP call */
79   	#define DEFAULT_MAXNLPCOST    1e+8           /**< default maximum cost available for solving NLPs per call of the heuristic */
80   	#define DEFAULT_MINIMPROVE    0.01           /**< default factor by which heuristic should at least improve the incumbent */
81   	#define DEFAULT_MAXNUNSUCC    10             /**< default maximum number of consecutive calls for which the heuristic did not find an improving solution */
82   	
83   	/*
84   	 * Data structures
85   	 */
86   	
87   	/** primal heuristic data */
88   	struct SCIP_HeurData
89   	{
90   	   SCIP_NLPI*            nlpi;               /**< nlpi used to create the nlpi problem */
91   	   SCIP_NLPIPROBLEM*     nlpiprob;           /**< nlpi problem representing the NLP relaxation */
92   	   SCIP_HASHMAP*         var2idx;            /**< mapping between variables and nlpi indices */
93   	   SCIP_HEUR*            subnlp;             /**< sub-NLP heuristic */
94   	
95   	   SCIP_Real             inittheta;          /**< initial regularization right-hand side value */
96   	   SCIP_Real             sigma;              /**< regularization update factor */
97   	   SCIP_Real             subnlptrigger;      /**< maximum number of NLP iterations per solve */
98   	   SCIP_Real             maxnlpcost;         /**< maximum cost available for solving NLPs per call of the heuristic */
99   	   SCIP_Real             minimprove;         /**< factor by which heuristic should at least improve the incumbent */
100  	   SCIP_Real             mingapleft;         /**< minimum amount of gap left in order to call the heuristic */
101  	   int                   maxiter;            /**< maximum number of iterations of the MPEC loop */
102  	   int                   maxnlpiter;         /**< maximum number of NLP iterations per solve */
103  	   int                   nunsucc;             /**< number of consecutive calls for which the heuristic did not find an
104  	                                              * improving solution */
105  	   int                   maxnunsucc;         /**< maximum number of consecutive calls for which the heuristic did not
106  	                                              * find an improving solution */
107  	};
108  	
109  	
110  	/*
111  	 * Local methods
112  	 */
113  	
114  	/** creates the data structure for generating the current NLP relaxation */
115  	static
116  	SCIP_RETCODE createNLP(
117  	   SCIP*                 scip,               /**< SCIP data structure */
118  	   SCIP_HEURDATA*        heurdata            /**< heuristic data */
119  	   )
120  	{
121  	   SCIP_Real cutoff = SCIPinfinity(scip);
122  	
123  	   assert(heurdata != NULL);
124  	   assert(heurdata->nlpi != NULL);
125  	
126  	   /* NLP has been already created */
127  	   if( heurdata->nlpiprob != NULL )
128  	      return SCIP_OKAY;
129  	
130  	   /* compute cutoff value to ensure minimum improvement */
131  	   if( SCIPgetNSols(scip) > 0 )
132  	   {
133  	      SCIP_Real upperbound = SCIPgetUpperbound(scip) - SCIPsumepsilon(scip);
134  	
135  	      assert( !SCIPisInfinity(scip, SCIPgetUpperbound(scip)) );
136  	
137  	      if( !SCIPisInfinity(scip, -SCIPgetLowerbound(scip)) )
138  	      {
139  	         cutoff = (1.0 - heurdata->minimprove) * SCIPgetUpperbound(scip)
140  	            + heurdata->minimprove * SCIPgetLowerbound(scip);
141  	      }
142  	      else
143  	      {
144  	         if( SCIPgetUpperbound(scip) >= 0.0 )
145  	            cutoff = ( 1.0 - heurdata->minimprove ) * SCIPgetUpperbound(scip);
146  	         else
147  	            cutoff = ( 1.0 + heurdata->minimprove ) * SCIPgetUpperbound(scip);
148  	      }
149  	      cutoff = MIN(upperbound, cutoff);
150  	      SCIPdebugMsg(scip, "set objective limit %g in [%g,%g]\n", cutoff, SCIPgetLowerbound(scip),
151  	         SCIPgetUpperbound(scip));
152  	   }
153  	
154  	   SCIP_CALL( SCIPhashmapCreate(&heurdata->var2idx, SCIPblkmem(scip), SCIPgetNVars(scip)) );
155  	   SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, heurdata->nlpi, &heurdata->nlpiprob, "MPEC-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip),
156  	         heurdata->var2idx, NULL, NULL, cutoff, TRUE, FALSE) );
157  	
158  	   return SCIP_OKAY;
159  	}
160  	
161  	/** frees the data structures for the NLP relaxation */
162  	static
163  	SCIP_RETCODE freeNLP(
164  	   SCIP*                 scip,               /**< SCIP data structure */
165  	   SCIP_HEURDATA*        heurdata            /**< heuristic data */
166  	   )
167  	{
168  	   assert(heurdata != NULL);
169  	
170  	   /* NLP has not been created yet */
171  	   if( heurdata->nlpiprob == NULL )
172  	      return SCIP_OKAY;
173  	
174  	   assert(heurdata->nlpi != NULL);
175  	   assert(heurdata->var2idx != NULL);
176  	
177  	   SCIPhashmapFree(&heurdata->var2idx);
178  	   SCIP_CALL( SCIPfreeNlpiProblem(scip, heurdata->nlpi, &heurdata->nlpiprob) );
179  	
180  	   return SCIP_OKAY;
181  	}
182  	
183  	/** adds or updates the regularization constraints to the NLP; for a given parameter theta we add for each non-fixed
184  	 *  binary variable z the constraint z*(1-z) <= theta; if these constraint are already present we update the theta on
185  	 *  the right-hand side
186  	 */
187  	static
188  	SCIP_RETCODE addRegularScholtes(
189  	   SCIP*                 scip,               /**< SCIP data structure */
190  	   SCIP_HEURDATA*        heurdata,           /**< heuristic data */
191  	   SCIP_VAR**            binvars,            /**< array containing all non-fixed binary variables */
192  	   int                   nbinvars,           /**< total number of non-fixed binary variables */
193  	   SCIP_Real             theta,              /**< regularization parameter */
194  	   SCIP_Bool             update              /**< should the regularization constraints be added or updated? */
195  	   )
196  	{
197  	   int i;
198  	
199  	   assert(binvars != NULL);
200  	   assert(nbinvars > 0);
201  	
202  	   /* add or update regularization for each non-fixed binary variables */
203  	   if( !update )
204  	   {
205  	      SCIP_NLROW** nlrows;
206  	
207  	      SCIP_CALL( SCIPallocBufferArray(scip, &nlrows, nbinvars) );
208  	
209  	      for( i = 0; i < nbinvars; ++i )
210  	      {
211  	         SCIP_Real one = 1.0;
212  	         SCIP_Real minusone = -1.0;
213  	         SCIP_EXPR* varexpr;
214  	         SCIP_EXPR* powexpr;
215  	         SCIP_EXPR* sumexpr;
216  	         char name[SCIP_MAXSTRLEN];
217  	
218  	         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_reg", SCIPvarGetName(binvars[i]));
219  	
220  	         /* -binvars[i]^2 */
221  	         SCIP_CALL( SCIPcreateExprVar(scip, &varexpr, binvars[i], NULL, NULL) );
222  	         SCIP_CALL( SCIPcreateExprPow(scip, &powexpr, varexpr, 2.0, NULL, NULL) );
223  	         SCIP_CALL( SCIPcreateExprSum(scip, &sumexpr, 1, &powexpr, &minusone, 0.0, NULL, NULL) );
224  	
225  	         /* binvars[i] - binvars[i]^2 <= theta */
226  	         SCIP_CALL( SCIPcreateNlRow(scip, &nlrows[i], name, 0.0, 1, &binvars[i], &one, sumexpr, -SCIPinfinity(scip), theta, SCIP_EXPRCURV_CONCAVE) );
227  	
228  	         SCIP_CALL( SCIPreleaseExpr(scip, &sumexpr) );
229  	         SCIP_CALL( SCIPreleaseExpr(scip, &powexpr) );
230  	         SCIP_CALL( SCIPreleaseExpr(scip, &varexpr) );
231  	      }
232  	
233  	      SCIP_CALL( SCIPaddNlpiProblemNlRows(scip, heurdata->nlpi, heurdata->nlpiprob, heurdata->var2idx, nlrows, nbinvars) );
234  	
235  	      for( i = nbinvars-1; i >= 0; --i )
236  	      {
237  	         SCIP_CALL( SCIPreleaseNlRow(scip, &nlrows[i]) );
238  	      }
239  	
240  	      SCIPfreeBufferArray(scip, &nlrows);
241  	   }
242  	   else
243  	   {
244  	      int startidx = SCIPgetNNLPNlRows(scip) + 1; /* the cutoff is a separate constraint */
245  	      SCIP_Real* lhss;
246  	      SCIP_Real* rhss;
247  	      int* indices;
248  	
249  	      SCIP_CALL( SCIPallocBufferArray(scip, &lhss, nbinvars) );
250  	      SCIP_CALL( SCIPallocBufferArray(scip, &rhss, nbinvars) );
251  	      SCIP_CALL( SCIPallocBufferArray(scip, &indices, nbinvars) );
252  	
253  	      for( i = 0; i < nbinvars; ++i )
254  	      {
255  	         lhss[i] = -SCIPinfinity(scip);
256  	         rhss[i] = theta;
257  	         indices[i] = startidx + i;
258  	      }
259  	
260  	      SCIP_CALL( SCIPchgNlpiConsSides(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lhss, rhss) );
261  	
262  	      SCIPfreeBufferArray(scip, &indices);
263  	      SCIPfreeBufferArray(scip, &rhss);
264  	      SCIPfreeBufferArray(scip, &lhss);
265  	   }
266  	
267  	   return SCIP_OKAY;
268  	}
269  	
270  	/** recursive helper function to count the number of nodes in a sub-expr */
271  	static
272  	int getExprSize(
273  	   SCIP_EXPR*            expr                /**< expression */
274  	   )
275  	{
276  	   int sum;
277  	   int i;
278  	
279  	   assert(expr != NULL);
280  	
281  	   sum = 0;
282  	   for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
283  	   {
284  	      SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
285  	      sum += getExprSize(child);
286  	   }
287  	   return 1 + sum;
288  	}
289  	
290  	/** main execution function of the MPEC heuristic */
291  	static
292  	SCIP_RETCODE heurExec(
293  	   SCIP*                 scip,               /**< SCIP data structure */
294  	   SCIP_HEUR*            heur,               /**< MPEC heuristic */
295  	   SCIP_HEURDATA*        heurdata,           /**< heuristic data */
296  	   SCIP_RESULT*          result              /**< pointer to store the result */
297  	   )
298  	{
299  	   SCIP_NLPSTATISTICS nlpstatistics;
300  	   SCIP_NLPPARAM nlpparam = SCIP_NLPPARAM_DEFAULT(scip);  /*lint !e446*/
301  	   SCIP_VAR** binvars = NULL;
302  	   SCIP_Real* initguess = NULL;
303  	   SCIP_Real* ubs = NULL;
304  	   SCIP_Real* lbs = NULL;
305  	   int* indices = NULL;
306  	   SCIP_Real theta = heurdata->inittheta;
307  	   SCIP_Real nlpcostperiter = 0.0;
308  	   SCIP_Real nlpcostleft = heurdata->maxnlpcost;
309  	   SCIP_Bool reinit = TRUE;
310  	   SCIP_Bool fixed = FALSE;
311  	   SCIP_Bool subnlpcalled = FALSE;
312  	   int nbinvars = 0;
313  	   int i;
314  	
315  	   assert(heurdata->nlpiprob != NULL);
316  	   assert(heurdata->var2idx != NULL);
317  	   assert(heurdata->nlpi != NULL);
318  	   assert(result != NULL);
319  	
320  	   SCIP_CALL( SCIPallocBufferArray(scip, &binvars, SCIPgetNBinVars(scip)) );
321  	
322  	   /* collect all non-fixed binary variables */
323  	   for( i = 0; i < SCIPgetNBinVars(scip); ++i )
324  	   {
325  	      SCIP_VAR* var = SCIPgetVars(scip)[i];
326  	      assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY);
327  	
328  	      if( !SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)) )
329  	         binvars[nbinvars++] = var;
330  	   }
331  	
332  	   /* all binary variables are fixed */
333  	   SCIPdebugMsg(scip, "nbinvars %d\n", nbinvars);
334  	   if( nbinvars == 0 )
335  	      goto TERMINATE;
336  	
337  	   SCIP_CALL( SCIPallocBufferArray(scip, &initguess, SCIPgetNVars(scip)) );
338  	   SCIP_CALL( SCIPallocBufferArray(scip, &lbs, nbinvars) );
339  	   SCIP_CALL( SCIPallocBufferArray(scip, &ubs, nbinvars) );
340  	   SCIP_CALL( SCIPallocBufferArray(scip, &indices, nbinvars) );
341  	
342  	   /* compute estimate cost for each NLP iteration */
343  	   for( i = 0; i < SCIPgetNNLPNlRows(scip); ++i )
344  	   {
345  	      SCIP_NLROW* nlrow = SCIPgetNLPNlRows(scip)[i];
346  	      assert(nlrow != NULL);
347  	
348  	      nlpcostperiter += 1.0 * SCIPnlrowGetNLinearVars(nlrow);
349  	
350  	      if( SCIPnlrowGetExpr(nlrow) != NULL )
351  	         nlpcostperiter += 3.0 * getExprSize(SCIPnlrowGetExpr(nlrow));
352  	   }
353  	
354  	   /* set initial guess */
355  	   for( i = 0; i < SCIPgetNVars(scip); ++i )
356  	   {
357  	      SCIP_VAR* var = SCIPgetVars(scip)[i];
358  	      initguess[i] = SCIPgetSolVal(scip, NULL, var);
359  	      /* SCIPdebugMsg(scip, "set initial value for %s to %g\n", SCIPvarGetName(var), initguess[i]); */
360  	   }
361  	   SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, initguess, NULL, NULL, NULL) );
362  	
363  	   /* set parameters of NLP solver */
364  	   nlpparam.feastol /= 10.0;
365  	   nlpparam.opttol /= 10.0;
366  	   nlpparam.iterlimit = heurdata->maxnlpiter;
367  	
368  	   /* main loop */
369  	   for( i = 0; i < heurdata->maxiter && *result != SCIP_FOUNDSOL && nlpcostleft > 0.0 && !SCIPisStopped(scip); ++i )
370  	   {
371  	      SCIP_Real* primal = NULL;
372  	      SCIP_Bool binaryfeasible;
373  	      SCIP_Bool regularfeasible;
374  	      SCIP_NLPSOLSTAT solstat;
375  	      SCIP_Real maxviolbin = 0.0;
376  	      SCIP_Real maxviolreg = 0.0;
377  	      int j;
378  	
379  	      /* add or update regularization */
380  	      SCIP_CALL( addRegularScholtes(scip, heurdata, binvars, nbinvars, theta, i > 0) );
381  	
382  	      /* solve NLP */
383  	      SCIP_CALL( SCIPsolveNlpiParam(scip, heurdata->nlpi, heurdata->nlpiprob, nlpparam) );
384  	      solstat = SCIPgetNlpiSolstat(scip, heurdata->nlpi, heurdata->nlpiprob);
385  	
386  	      /* give up if an error occurred or no primal values are accessible */
387  	      if( solstat > SCIP_NLPSOLSTAT_LOCINFEASIBLE )
388  	      {
389  	         SCIPdebugMsg(scip, "error occured during NLP solve -> stop!\n");
390  	         break;
391  	      }
392  	
393  	      /* update nlpcostleft */
394  	      SCIP_CALL( SCIPgetNlpiStatistics(scip, heurdata->nlpi, heurdata->nlpiprob, &nlpstatistics) );
395  	      nlpcostleft -= nlpstatistics.niterations * nlpcostperiter * nbinvars;
396  	      SCIPdebugMsg(scip, "nlpcostleft = %e\n", nlpcostleft);
397  	
398  	      SCIP_CALL( SCIPgetNlpiSolution(scip, heurdata->nlpi, heurdata->nlpiprob, &primal, NULL, NULL, NULL, NULL) );
399  	      assert(primal != NULL);
400  	
401  	      /* check for binary feasibility */
402  	      binaryfeasible = TRUE;
403  	      regularfeasible = TRUE;
404  	      for( j = 0; j < nbinvars; ++j )
405  	      {
406  	         int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]);
407  	         binaryfeasible = binaryfeasible && SCIPisFeasIntegral(scip, primal[idx]);
408  	         regularfeasible = regularfeasible && SCIPisLE(scip, primal[idx] - SQR(primal[idx]), theta);
409  	
410  	         maxviolreg = MAX(maxviolreg, primal[idx] - SQR(primal[idx]) - theta);
411  	         maxviolbin = MAX(maxviolbin, MIN(primal[idx], 1.0-primal[idx]));
412  	      }
413  	      SCIPdebugMsg(scip, "maxviol-regularization %g maxviol-integrality %g\n", maxviolreg, maxviolbin);
414  	
415  	      /* call sub-NLP heuristic when the maximum binary infeasibility is small enough (or this is the last iteration
416  	       * because we reached the nlpcost limit)
417  	       */
418  	      if( !subnlpcalled && heurdata->subnlp != NULL
419  	         && (SCIPisLE(scip, maxviolbin, heurdata->subnlptrigger) || nlpcostleft <= 0.0)
420  	         && !SCIPisStopped(scip) )
421  	      {
422  	         SCIP_SOL* refpoint;
423  	         SCIP_RESULT subnlpresult;
424  	
425  	         SCIPdebugMsg(scip, "call sub-NLP heuristic because binary infeasibility is small enough\n");
426  	         SCIP_CALL( SCIPcreateSol(scip, &refpoint, heur) );
427  	
428  	         for( j = 0; j < SCIPgetNVars(scip); ++j )
429  	         {
430  	            SCIP_VAR* var = SCIPgetVars(scip)[j];
431  	            SCIP_Real val = SCIPvarIsBinary(var) ? SCIPfeasRound(scip, primal[j]) : primal[j];
432  	            SCIP_CALL( SCIPsetSolVal(scip, refpoint, var, val) );
433  	         }
434  	
435  	         SCIP_CALL( SCIPapplyHeurSubNlp(scip, heurdata->subnlp, &subnlpresult, refpoint, NULL) );
436  	         SCIP_CALL( SCIPfreeSol(scip, &refpoint) );
437  	         SCIPdebugMsg(scip, "result of sub-NLP call: %d\n", subnlpresult);
438  	
439  	         /* stop MPEC heuristic when the sub-NLP heuristic has found a feasible solution */
440  	         if( subnlpresult == SCIP_FOUNDSOL )
441  	         {
442  	            SCIPdebugMsg(scip, "sub-NLP found a feasible solution -> stop!\n");
443  	            break;
444  	         }
445  	
446  	         subnlpcalled = TRUE;
447  	      }
448  	
449  	      /* NLP feasible + binary feasible -> add solution and stop */
450  	      if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE && binaryfeasible )
451  	      {
452  	         SCIP_SOL* sol;
453  	         SCIP_Bool stored;
454  	
455  	         SCIP_CALL( SCIPcreateSol(scip, &sol, heur) );
456  	
457  	         for( j = 0; j < SCIPgetNVars(scip); ++j )
458  	         {
459  	            SCIP_VAR* var = SCIPgetVars(scip)[j];
460  	            assert(j == SCIPhashmapGetImageInt(heurdata->var2idx, (void*)var));
461  	            SCIP_CALL( SCIPsetSolVal(scip, sol, var, primal[j]) );
462  	         }
463  	
464  	#ifdef SCIP_DEBUG
465  	         SCIP_CALL( SCIPtrySolFree(scip, &sol, TRUE, TRUE, TRUE, TRUE, FALSE, &stored) );
466  	#else
467  	         SCIP_CALL( SCIPtrySolFree(scip, &sol, FALSE, FALSE, TRUE, TRUE, FALSE, &stored) );
468  	#endif
469  	         SCIPdebugMsg(scip, "found a solution (stored = %u)\n", stored);
470  	
471  	         if( stored )
472  	            *result = SCIP_FOUNDSOL;
473  	         break;
474  	      }
475  	
476  	      /* NLP feasible + binary infeasible -> reduce theta */
477  	      else if( solstat <= SCIP_NLPSOLSTAT_FEASIBLE && !binaryfeasible )
478  	      {
479  	         BMScopyMemoryArray(initguess, primal, SCIPgetNVars(scip));
480  	         SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, primal, NULL, NULL, NULL) );
481  	         SCIPdebugMsg(scip, "update theta from %g -> %g\n", theta, theta*heurdata->sigma);
482  	
483  	         if( !reinit )
484  	         {
485  	            SCIPdebugMsg(scip, "reinit fixed the infeasibility\n");
486  	            reinit = TRUE;
487  	         }
488  	
489  	         theta *= heurdata->sigma;
490  	
491  	         /* unfix binary variables */
492  	         if( fixed )
493  	         {
494  	            SCIPdebugMsg(scip, "unfixing binary variables\n");
495  	            for( j = 0; j < nbinvars; ++j )
496  	            {
497  	               lbs[j] = 0.0;
498  	               ubs[j] = 1.0;
499  	               indices[j] = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]);
500  	            }
501  	            SCIP_CALL( SCIPchgNlpiVarBounds(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lbs, ubs) );
502  	            fixed = FALSE;
503  	         }
504  	      }
505  	
506  	      /* NLP infeasible + regularization feasible -> stop (give up) */
507  	      else if( solstat > SCIP_NLPSOLSTAT_FEASIBLE && regularfeasible )
508  	      {
509  	         SCIPdebugMsg(scip, "NLP is infeasible but regularization constraints are satisfied -> stop!\n");
510  	         break;
511  	      }
512  	
513  	      /* NLP infeasible + binary infeasible -> set initial point / fix binary variables */
514  	      else
515  	      {
516  	         assert(solstat > SCIP_NLPSOLSTAT_FEASIBLE && !regularfeasible);
517  	
518  	         SCIPdebugMsg(scip, "NLP solution is not feasible for the NLP and the binary variables\n");
519  	
520  	         /* stop if fixing did not resolve the infeasibility */
521  	         if( fixed )
522  	         {
523  	            SCIPdebugMsg(scip, "fixing variables did not resolve infeasibility -> stop!\n");
524  	            break;
525  	         }
526  	
527  	         /* fix variables if reinit is FALSE; otherwise set another initial point */
528  	         if( !reinit )
529  	         {
530  	            int nfixedvars = 0;
531  	
532  	            /* fix binary variables */
533  	            for( j = 0; j < nbinvars; ++j )
534  	            {
535  	               int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]);
536  	               indices[j] = idx;
537  	
538  	               if( SCIPisFeasLE(scip, primal[idx] - SQR(primal[idx]), theta) )
539  	               {
540  	                  lbs[j] = 0.0;
541  	                  ubs[j] = 1.0;
542  	               }
543  	               else
544  	               {
545  	                  lbs[j] = primal[idx] >= 0.5 ? 0.0 : 1.0;
546  	                  ubs[j] = primal[idx] >= 0.5 ? 0.0 : 1.0;
547  	                  ++nfixedvars;
548  	                  /* SCIPdebugMsg(scip, "fix binary variable %s = %g\n", SCIPvarGetName(binvars[j]), ubs[j]); */
549  	               }
550  	            }
551  	            SCIPdebugMsg(scip, "fixed %d binary variables\n", nfixedvars);
552  	            SCIP_CALL( SCIPchgNlpiVarBounds(scip, heurdata->nlpi, heurdata->nlpiprob, nbinvars, indices, lbs, ubs) );
553  	            fixed = TRUE;
554  	         }
555  	         else
556  	         {
557  	            SCIPdebugMsg(scip, "update initial guess\n");
558  	
559  	            /* set initial point */
560  	            for( j = 0; j < nbinvars; ++j )
561  	            {
562  	               int idx = SCIPhashmapGetImageInt(heurdata->var2idx, (void*)binvars[j]);
563  	               initguess[idx] = primal[idx] >= 0.5 ? 0.0 : 1.0;
564  	               /* SCIPdebugMsg(scip, "update init guess for %s to %g\n", SCIPvarGetName(binvars[j]), initguess[idx]); */
565  	            }
566  	            SCIP_CALL( SCIPsetNlpiInitialGuess(scip, heurdata->nlpi, heurdata->nlpiprob, initguess, NULL, NULL, NULL) );
567  	            reinit = FALSE;
568  	         }
569  	      }
570  	   }
571  	
572  	TERMINATE:
573  	   SCIPfreeBufferArrayNull(scip, &indices);
574  	   SCIPfreeBufferArrayNull(scip, &ubs);
575  	   SCIPfreeBufferArrayNull(scip, &lbs);
576  	   SCIPfreeBufferArrayNull(scip, &initguess);
577  	   SCIPfreeBufferArray(scip, &binvars);
578  	
579  	   return SCIP_OKAY;
580  	}
581  	
582  	
583  	/*
584  	 * Callback methods of primal heuristic
585  	 */
586  	
587  	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
588  	static
589  	SCIP_DECL_HEURCOPY(heurCopyMpec)
590  	{  /*lint --e{715}*/
591  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
592  	
593  	   /* call inclusion method of primal heuristic */
594  	   SCIP_CALL( SCIPincludeHeurMpec(scip) );
595  	
596  	   return SCIP_OKAY;
597  	}
598  	
599  	/** destructor of primal heuristic to free user data (called when SCIP is exiting) */
600  	static
601  	SCIP_DECL_HEURFREE(heurFreeMpec)
602  	{  /*lint --e{715}*/
603  	   SCIP_HEURDATA* heurdata = SCIPheurGetData(heur);
604  	   assert(heurdata != NULL);
605  	
606  	   SCIPfreeBlockMemory(scip, &heurdata);
607  	   SCIPheurSetData(heur, NULL);
608  	
609  	   return SCIP_OKAY;
610  	}
611  	
612  	/** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
613  	static
614  	SCIP_DECL_HEURINITSOL(heurInitsolMpec)
615  	{  /*lint --e{715}*/
616  	   SCIP_HEURDATA* heurdata = SCIPheurGetData(heur);
617  	
618  	   assert(heurdata != NULL);
619  	   assert(heurdata->nlpi == NULL);
620  	
621  	   if( SCIPgetNNlpis(scip) > 0 )
622  	   {
623  	      heurdata->nlpi = SCIPgetNlpis(scip)[0];
624  	      heurdata->subnlp = SCIPfindHeur(scip, "subnlp");
625  	   }
626  	
627  	   return SCIP_OKAY;
628  	}
629  	
630  	
631  	/** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */
632  	static
633  	SCIP_DECL_HEUREXITSOL(heurExitsolMpec)
634  	{  /*lint --e{715}*/
635  	   SCIP_HEURDATA* heurdata = SCIPheurGetData(heur);
636  	
637  	   assert(heurdata != NULL);
638  	   heurdata->nlpi = NULL;
639  	
640  	   return SCIP_OKAY;
641  	}
642  	
643  	/** execution method of primal heuristic */
644  	static
645  	SCIP_DECL_HEUREXEC(heurExecMpec)
646  	{  /*lint --e{715}*/
647  	   SCIP_HEURDATA* heurdata = SCIPheurGetData(heur);
648  	   SCIP_CONSHDLR* sosonehdlr = SCIPfindConshdlr(scip, "SOS1");
649  	   SCIP_CONSHDLR* sostwohdlr = SCIPfindConshdlr(scip, "SOS2");
650  	
651  	   assert(heurdata != NULL);
652  	
653  	   *result = SCIP_DIDNOTRUN;
654  	
655  	   if( SCIPgetNIntVars(scip) > 0 || SCIPgetNBinVars(scip) == 0
656  	      || heurdata->nlpi == NULL || !SCIPisNLPConstructed(scip)
657  	      || heurdata->mingapleft > SCIPgetGap(scip)
658  	      || heurdata->nunsucc > heurdata->maxnunsucc )
659  	      return SCIP_OKAY;
660  	
661  	   /* skip heuristic if constraints without a nonlinear representation are present */
662  	   if( (sosonehdlr != NULL && SCIPconshdlrGetNConss(sosonehdlr) > 0) ||
663  	      (sostwohdlr != NULL && SCIPconshdlrGetNConss(sostwohdlr) > 0) )
664  	   {
665  	      return SCIP_OKAY;
666  	   }
667  	
668  	   *result = SCIP_DIDNOTFIND;
669  	
670  	   /* call MPEC method */
671  	   SCIP_CALL( createNLP(scip, heurdata) );
672  	   SCIP_CALL( heurExec(scip, heur, heurdata, result) );
673  	   SCIP_CALL( freeNLP(scip, heurdata) );
674  	
675  	   /* update number of unsuccessful calls */
676  	   heurdata->nunsucc = (*result == SCIP_FOUNDSOL) ? 0 : heurdata->nunsucc + 1;
677  	
678  	   return SCIP_OKAY;
679  	}
680  	
681  	
682  	/*
683  	 * primal heuristic specific interface methods
684  	 */
685  	
686  	/** creates the mpec primal heuristic and includes it in SCIP */
687  	SCIP_RETCODE SCIPincludeHeurMpec(
688  	   SCIP*                 scip                /**< SCIP data structure */
689  	   )
690  	{
691  	   SCIP_HEURDATA* heurdata = NULL;
692  	   SCIP_HEUR* heur = NULL;
693  	
694  	   /* create mpec primal heuristic data */
695  	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
696  	   BMSclearMemory(heurdata);
697  	
698  	   /* include primal heuristic */
699  	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
700  	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
701  	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecMpec, heurdata) );
702  	
703  	   assert(heur != NULL);
704  	
705  	   /* set non fundamental callbacks via setter functions */
706  	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyMpec) );
707  	   SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeMpec) );
708  	   SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolMpec) );
709  	   SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolMpec) );
710  	
711  	   /* add mpec primal heuristic parameters */
712  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/inittheta",
713  	         "initial regularization right-hand side value",
714  	         &heurdata->inittheta, FALSE, DEFAULT_INITTHETA, 0.0, 0.25, NULL, NULL) );
715  	
716  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/sigma",
717  	         "regularization update factor",
718  	         &heurdata->sigma, FALSE, DEFAULT_SIGMA, 0.0, 1.0, NULL, NULL) );
719  	
720  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/subnlptrigger",
721  	         "maximum number of NLP iterations per solve",
722  	         &heurdata->subnlptrigger, FALSE, DEFAULT_SUBNLPTRIGGER, 0.0, 1.0, NULL, NULL) );
723  	
724  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxnlpcost",
725  	         "maximum cost available for solving NLPs per call of the heuristic",
726  	         &heurdata->maxnlpcost, FALSE, DEFAULT_MAXNLPCOST, 0.0, SCIPinfinity(scip), NULL, NULL) );
727  	
728  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprove",
729  	         "factor by which heuristic should at least improve the incumbent",
730  	         &heurdata->minimprove, FALSE, DEFAULT_MINIMPROVE, 0.0, 1.0, NULL, NULL) );
731  	
732  	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/mingapleft",
733  	         "minimum amount of gap left in order to call the heuristic",
734  	         &heurdata->mingapleft, FALSE, DEFAULT_MINGAPLEFT, 0.0, SCIPinfinity(scip), NULL, NULL) );
735  	
736  	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiter",
737  	         "maximum number of iterations of the MPEC loop",
738  	         &heurdata->maxiter, FALSE, DEFAULT_MAXITER, 0, INT_MAX, NULL, NULL) );
739  	
740  	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxnlpiter",
741  	         "maximum number of NLP iterations per solve",
742  	         &heurdata->maxnlpiter, FALSE, DEFAULT_MAXNLPITER, 0, INT_MAX, NULL, NULL) );
743  	
744  	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxnunsucc",
745  	         "maximum number of consecutive calls for which the heuristic did not find an improving solution",
746  	         &heurdata->maxnunsucc, FALSE, DEFAULT_MAXNUNSUCC, 0, INT_MAX, NULL, NULL) );
747  	
748  	   return SCIP_OKAY;
749  	}
750