1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file heuristics.c
17 * @ingroup OTHER_CFILES
18 * @brief methods commonly used by primal heuristics
19 * @author Gregor Hendel
20 */
21
22 #include "scip/heuristics.h"
23 #include "scip/cons_linear.h"
24 #include "scip/scipdefplugins.h"
25
26 #include "scip/pub_heur.h"
27
28 /* the indicator and SOS1 constraint handlers are included for the diving algorithm SCIPperformGenericDivingAlgorithm() */
29 #include "scip/cons_indicator.h"
30 #include "scip/cons_sos1.h"
31
32 #define MINLPITER 10000 /**< minimal number of LP iterations allowed in each LP solving call */
33
34
35 /** solve probing LP */
36 static
solveLP(SCIP * scip,SCIP_DIVESET * diveset,SCIP_Longint maxnlpiterations,SCIP_DIVECONTEXT divecontext,SCIP_Bool * lperror,SCIP_Bool * cutoff)37 SCIP_RETCODE solveLP(
38 SCIP* scip, /**< SCIP data structure */
39 SCIP_DIVESET* diveset, /**< diving settings */
40 SCIP_Longint maxnlpiterations, /**< maximum number of allowed LP iterations */
41 SCIP_DIVECONTEXT divecontext, /**< context for diving statistics */
42 SCIP_Bool* lperror, /**< pointer to store if an internal LP error occurred */
43 SCIP_Bool* cutoff /**< pointer to store whether the LP was infeasible */
44 )
45 {
46 int lpiterationlimit;
47 SCIP_RETCODE retstat;
48 SCIP_Longint nlpiterations;
49
50 assert(lperror != NULL);
51 assert(cutoff != NULL);
52
53 nlpiterations = SCIPgetNLPIterations(scip);
54
55 /* allow at least MINLPITER more iterations so as not to run out of LP iterations during this solve */
56 lpiterationlimit = (int)(maxnlpiterations - SCIPdivesetGetNLPIterations(diveset, divecontext));
57 lpiterationlimit = MAX(lpiterationlimit, MINLPITER);
58
59 retstat = SCIPsolveProbingLP(scip, lpiterationlimit, lperror, cutoff);
60
61 /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic.
62 * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop.
63 */
64 #ifdef NDEBUG
65 if( retstat != SCIP_OKAY )
66 {
67 SCIPwarningMessage(scip, "Error while solving LP in %s diving heuristic; LP solve terminated with code <%d>.\n", SCIPdivesetGetName(diveset), retstat);
68 }
69 #else
70 SCIP_CALL( retstat );
71 #endif
72
73 /* update iteration count */
74 SCIPupdateDivesetLPStats(scip, diveset, SCIPgetNLPIterations(scip) - nlpiterations, divecontext);
75
76 return SCIP_OKAY;
77 }
78
79 /** select the next variable and type of diving */
80 static
selectNextDiving(SCIP * scip,SCIP_DIVESET * diveset,SCIP_SOL * worksol,SCIP_Bool onlylpbranchcands,SCIP_Bool storelpcandscores,SCIP_VAR ** lpcands,SCIP_Real * lpcandssol,SCIP_Real * lpcandsfrac,SCIP_Real * lpcandsscores,SCIP_Bool * lpcandroundup,int * nviollpcands,int nlpcands,SCIP_Bool * enfosuccess,SCIP_Bool * infeasible)81 SCIP_RETCODE selectNextDiving(
82 SCIP* scip, /**< SCIP data structure */
83 SCIP_DIVESET* diveset, /**< dive set */
84 SCIP_SOL* worksol, /**< current working solution */
85 SCIP_Bool onlylpbranchcands, /**< should only LP branching candidates be considered? */
86 SCIP_Bool storelpcandscores, /**< should the scores of the LP candidates be updated? */
87 SCIP_VAR** lpcands, /**< LP branching candidates, or NULL if not needed */
88 SCIP_Real * lpcandssol, /**< solution values LP branching candidates, or NULL if not needed */
89 SCIP_Real* lpcandsfrac, /**< fractionalities of LP branching candidates, or NULL if not needed*/
90 SCIP_Real* lpcandsscores, /**< array with LP branching candidate scores, or NULL */
91 SCIP_Bool* lpcandroundup, /**< array to remember whether the preferred branching direction is upwards */
92 int* nviollpcands, /**< pointer to store the number of LP candidates whose solution value already violates local bounds */
93 int nlpcands, /**< number of current LP cands */
94 SCIP_Bool* enfosuccess, /**< pointer to store whether a candidate was sucessfully found */
95 SCIP_Bool* infeasible /**< pointer to store whether the diving can be immediately aborted because it is infeasible */
96 )
97 {
98 assert(scip != NULL);
99 assert(worksol != NULL);
100 assert(!onlylpbranchcands || lpcandsscores != NULL);
101 assert(!onlylpbranchcands || lpcandroundup != NULL);
102 assert(enfosuccess != NULL);
103 assert(infeasible != NULL);
104 assert(nviollpcands != NULL);
105
106 *nviollpcands = 0;
107 /* we use diving solution enforcement provided by the constraint handlers */
108 if( !onlylpbranchcands )
109 {
110 SCIP_CALL( SCIPgetDiveBoundChanges(scip, diveset, worksol, enfosuccess, infeasible) );
111 }
112 else
113 {
114 int c;
115 int bestcandidx;
116 SCIP_Real bestscore;
117 SCIP_Real score;
118
119 bestscore = SCIP_REAL_MIN;
120 bestcandidx = -1;
121
122 SCIPclearDiveBoundChanges(scip);
123
124 /* search for the candidate that maximizes the dive set score function and whose solution value is still feasible */
125 for( c = 0; c < nlpcands; ++c )
126 {
127 assert(SCIPgetSolVal(scip, worksol, lpcands[c]) == lpcandssol[c]); /*lint !e777 doesn't like comparing floats for equality */
128
129 /* scores are kept in arrays for faster reuse */
130 if( storelpcandscores )
131 {
132 SCIP_CALL( SCIPgetDivesetScore(scip, diveset, SCIP_DIVETYPE_INTEGRALITY, lpcands[c], lpcandssol[c],
133 lpcandsfrac[c], &lpcandsscores[c], &lpcandroundup[c]) );
134 }
135
136 score = lpcandsscores[c];
137 /* update the best candidate if it has a higher score and a solution value which does not violate one of the local bounds */
138 if( SCIPisLE(scip, SCIPvarGetLbLocal(lpcands[c]), lpcandssol[c]) && SCIPisGE(scip, SCIPvarGetUbLocal(lpcands[c]), lpcandssol[c]) )
139 {
140 if( score > bestscore )
141 {
142 bestcandidx = c;
143 bestscore = score;
144 }
145 }
146 else
147 ++(*nviollpcands);
148 }
149
150 /* there is no guarantee that a candidate is found since local bounds might render all solution values infeasible */
151 *enfosuccess = (bestcandidx >= 0);
152 if( *enfosuccess )
153 {
154 /* if we want to round up the best candidate, it is added as the preferred bound change */
155 SCIP_CALL( SCIPaddDiveBoundChange(scip, lpcands[bestcandidx], SCIP_BRANCHDIR_UPWARDS,
156 SCIPceil(scip, lpcandssol[bestcandidx]), lpcandroundup[bestcandidx]) );
157 SCIP_CALL( SCIPaddDiveBoundChange(scip, lpcands[bestcandidx], SCIP_BRANCHDIR_DOWNWARDS,
158 SCIPfloor(scip, lpcandssol[bestcandidx]), ! lpcandroundup[bestcandidx]) );
159 }
160 }
161 return SCIP_OKAY;
162 }
163
164 /** return the LP iteration budget suggestion for this dive set */
165 static
getDivesetIterLimit(SCIP * scip,SCIP_DIVESET * diveset,SCIP_DIVECONTEXT divecontext)166 SCIP_Longint getDivesetIterLimit(
167 SCIP* scip, /**< SCIP data structure */
168 SCIP_DIVESET* diveset, /**< dive set data structure */
169 SCIP_DIVECONTEXT divecontext /**< context for diving statistics */
170 )
171 {
172 SCIP_Longint iterlimit;
173 /*todo another factor of 10, REALLY? */
174 iterlimit = (SCIP_Longint)((1.0 + 10*(SCIPdivesetGetNSols(diveset, divecontext)+1.0)/(SCIPdivesetGetNCalls(diveset, divecontext)+1.0)) * SCIPdivesetGetMaxLPIterQuot(diveset) * SCIPgetNNodeLPIterations(scip));
175 iterlimit += SCIPdivesetGetMaxLPIterOffset(diveset);
176 iterlimit -= SCIPdivesetGetNLPIterations(diveset, divecontext);
177
178 return iterlimit;
179 }
180
181 /** performs a diving within the limits of the @p diveset parameters
182 *
183 * This method performs a diving according to the settings defined by the diving settings @p diveset; Contrary to the
184 * name, SCIP enters probing mode (not diving mode) and dives along a path into the tree. Domain propagation
185 * is applied at every node in the tree, whereas probing LPs might be solved less frequently.
186 *
187 * Starting from the current LP solution, the algorithm selects candidates which maximize the
188 * score defined by the @p diveset and whose solution value has not yet been rendered infeasible by propagation,
189 * and propagates the bound change on this candidate.
190 *
191 * The algorithm iteratively selects the the next (unfixed) candidate in the list, until either enough domain changes
192 * or the resolve frequency of the LP trigger an LP resolve (and hence, the set of potential candidates changes),
193 * or the last node is proven to be infeasible. It optionally backtracks and tries the
194 * other branching direction.
195 *
196 * After the set of remaining candidates is empty or the targeted depth is reached, the node LP is
197 * solved, and the old candidates are replaced by the new LP candidates.
198 *
199 * @see heur_guideddiving.c for an example implementation of a dive set controlling the diving algorithm.
200 *
201 * @note the node from where the algorithm is called is checked for a basic LP solution. If the solution
202 * is non-basic, e.g., when barrier without crossover is used, the method returns without performing a dive.
203 *
204 * @note currently, when multiple diving heuristics call this method and solve an LP at the same node, only the first
205 * call will be executed, see SCIPgetLastDiveNode()
206 *
207 * @todo generalize method to work correctly with pseudo or external branching/diving candidates
208 */
SCIPperformGenericDivingAlgorithm(SCIP * scip,SCIP_DIVESET * diveset,SCIP_SOL * worksol,SCIP_HEUR * heur,SCIP_RESULT * result,SCIP_Bool nodeinfeasible,SCIP_Longint iterlim,SCIP_DIVECONTEXT divecontext)209 SCIP_RETCODE SCIPperformGenericDivingAlgorithm(
210 SCIP* scip, /**< SCIP data structure */
211 SCIP_DIVESET* diveset, /**< settings for diving */
212 SCIP_SOL* worksol, /**< non-NULL working solution */
213 SCIP_HEUR* heur, /**< the calling primal heuristic */
214 SCIP_RESULT* result, /**< SCIP result pointer */
215 SCIP_Bool nodeinfeasible, /**< is the current node known to be infeasible? */
216 SCIP_Longint iterlim, /**< nonnegative iteration limit for the LP solves, or -1 for dynamic setting */
217 SCIP_DIVECONTEXT divecontext /**< context for diving statistics */
218 )
219 {
220 SCIP_CONSHDLR* indconshdlr; /* constraint handler for indicator constraints */
221 SCIP_CONSHDLR* sos1conshdlr; /* constraint handler for SOS1 constraints */
222 SCIP_VAR** lpcands;
223 SCIP_Real* lpcandssol;
224
225 SCIP_VAR** previouscands;
226 SCIP_Real* lpcandsscores;
227 SCIP_Real* previousvals;
228 SCIP_Real* lpcandsfrac;
229 SCIP_Bool* lpcandroundup;
230 SCIP_Real searchubbound;
231 SCIP_Real searchavgbound;
232 SCIP_Real searchbound;
233 SCIP_Real ubquot;
234 SCIP_Real avgquot;
235 SCIP_Longint maxnlpiterations;
236 SCIP_Longint domreds;
237 int startndivecands;
238 int depth;
239 int maxdepth;
240 int maxdivedepth;
241 int totalnbacktracks;
242 int totalnprobingnodes;
243 int lastlpdepth;
244 int previouscandssize;
245 int lpcandsscoressize;
246 int nviollpcands;
247 SCIP_Longint oldnsolsfound;
248 SCIP_Longint oldnbestsolsfound;
249 SCIP_Longint oldnconflictsfound;
250
251 SCIP_Bool success;
252 SCIP_Bool leafsol;
253 SCIP_Bool enfosuccess;
254 SCIP_Bool lperror;
255 SCIP_Bool cutoff;
256 SCIP_Bool backtracked;
257 SCIP_Bool backtrack;
258 SCIP_Bool onlylpbranchcands;
259
260 int nlpcands;
261 int lpsolvefreq;
262
263 assert(scip != NULL);
264 assert(heur != NULL);
265 assert(result != NULL);
266 assert(worksol != NULL);
267
268 *result = SCIP_DELAYED;
269
270 /* do not call heuristic in node that was already detected to be infeasible */
271 if( nodeinfeasible )
272 return SCIP_OKAY;
273
274 /* only call heuristic, if an optimal LP solution is at hand */
275 if( !SCIPhasCurrentNodeLP(scip) || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
276 return SCIP_OKAY;
277
278 /* only call heuristic, if the LP objective value is smaller than the cutoff bound */
279 if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) )
280 return SCIP_OKAY;
281
282 /* only call heuristic, if the LP solution is basic (which allows fast resolve in diving) */
283 if( !SCIPisLPSolBasic(scip) )
284 return SCIP_OKAY;
285
286 /* don't dive two times at the same node */
287 if( SCIPgetLastDivenode(scip) == SCIPgetNNodes(scip) && SCIPgetDepth(scip) > 0 )
288 return SCIP_OKAY;
289
290 *result = SCIP_DIDNOTRUN;
291
292 /* only try to dive, if we are in the correct part of the tree, given by minreldepth and maxreldepth */
293 depth = SCIPgetDepth(scip);
294 maxdepth = SCIPgetMaxDepth(scip);
295 maxdepth = MAX(maxdepth, 30);
296 if( depth < SCIPdivesetGetMinRelDepth(diveset) * maxdepth || depth > SCIPdivesetGetMaxRelDepth(diveset) * maxdepth )
297 return SCIP_OKAY;
298
299 /* calculate the maximal number of LP iterations */
300 if( iterlim < 0 )
301 {
302 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + getDivesetIterLimit(scip, diveset, divecontext);
303 }
304 else
305 {
306 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + iterlim;
307 }
308
309 /* don't try to dive, if we took too many LP iterations during diving */
310 if( SCIPdivesetGetNLPIterations(diveset, divecontext) >= maxnlpiterations )
311 return SCIP_OKAY;
312
313 /* allow at least a certain number of LP iterations in this dive */
314 if( SCIPdivesetGetNLPIterations(diveset, divecontext) + MINLPITER > maxnlpiterations )
315 maxnlpiterations = SCIPdivesetGetNLPIterations(diveset, divecontext) + MINLPITER;
316
317 /* these constraint handlers are required for polishing an LP relaxation solution beyond rounding */
318 indconshdlr = SCIPfindConshdlr(scip, "indicator");
319 sos1conshdlr = SCIPfindConshdlr(scip, "SOS1");
320
321 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, &lpcandsfrac, &nlpcands, NULL, NULL) );
322
323 onlylpbranchcands = SCIPdivesetUseOnlyLPBranchcands(diveset);
324 /* don't try to dive, if there are no diving candidates */
325 if( onlylpbranchcands && nlpcands == 0 )
326 return SCIP_OKAY;
327
328 /* calculate the objective search bound */
329 if( SCIPgetNSolsFound(scip) == 0 )
330 {
331 ubquot = SCIPdivesetGetUbQuotNoSol(diveset);
332 avgquot = SCIPdivesetGetAvgQuotNoSol(diveset);
333 }
334 else
335 {
336 ubquot = SCIPdivesetGetUbQuot(diveset);
337 avgquot = SCIPdivesetGetAvgQuot(diveset);
338 }
339
340 if( ubquot > 0.0 )
341 searchubbound = SCIPgetLowerbound(scip) + ubquot * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip));
342 else
343 searchubbound = SCIPinfinity(scip);
344
345 if( avgquot > 0.0 )
346 searchavgbound = SCIPgetLowerbound(scip) + avgquot * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip));
347 else
348 searchavgbound = SCIPinfinity(scip);
349
350 searchbound = MIN(searchubbound, searchavgbound);
351
352 if( SCIPisObjIntegral(scip) )
353 searchbound = SCIPceil(scip, searchbound);
354
355 /* calculate the maximal diving depth: 10 * min{number of integer variables, max depth} */
356 maxdivedepth = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
357 if ( sos1conshdlr != NULL )
358 maxdivedepth += SCIPgetNSOS1Vars(sos1conshdlr);
359 maxdivedepth = MIN(maxdivedepth, maxdepth);
360 maxdivedepth *= 10;
361
362 *result = SCIP_DIDNOTFIND;
363
364 /* start probing mode */
365 SCIP_CALL( SCIPstartProbing(scip) );
366
367 /* enables collection of variable statistics during probing */
368 SCIPenableVarHistory(scip);
369
370 SCIPdebugMsg(scip, "(node %" SCIP_LONGINT_FORMAT ") executing %s heuristic: depth=%d, %d fractionals, dualbound=%g, avgbound=%g, cutoffbound=%g, searchbound=%g\n",
371 SCIPgetNNodes(scip), SCIPheurGetName(heur), SCIPgetDepth(scip), nlpcands, SCIPgetDualbound(scip), SCIPgetAvgDualbound(scip),
372 SCIPretransformObj(scip, SCIPgetCutoffbound(scip)), SCIPretransformObj(scip, searchbound));
373
374 /* allocate buffer storage for previous candidates and their branching values for pseudo cost updates */
375 lpsolvefreq = SCIPdivesetGetLPSolveFreq(diveset);
376 previouscandssize = MAX(1, lpsolvefreq);
377 SCIP_CALL( SCIPallocBufferArray(scip, &previouscands, previouscandssize) );
378 SCIP_CALL( SCIPallocBufferArray(scip, &previousvals, previouscandssize) );
379
380 /* keep some statistics */
381 lperror = FALSE;
382 cutoff = FALSE;
383 lastlpdepth = -1;
384 startndivecands = nlpcands;
385 totalnbacktracks = 0;
386 totalnprobingnodes = 0;
387 oldnsolsfound = SCIPgetNSolsFound(scip);
388 oldnbestsolsfound = SCIPgetNBestSolsFound(scip);
389 oldnconflictsfound = SCIPgetNConflictConssFound(scip);
390
391 /* link the working solution to the dive set */
392 SCIPdivesetSetWorkSolution(diveset, worksol);
393
394 if( onlylpbranchcands )
395 {
396 SCIP_CALL( SCIPallocBufferArray(scip, &lpcandsscores, nlpcands) );
397 SCIP_CALL( SCIPallocBufferArray(scip, &lpcandroundup, nlpcands) );
398
399 lpcandsscoressize = nlpcands;
400 }
401 else
402 {
403 lpcandroundup = NULL;
404 lpcandsscores = NULL;
405 lpcandsscoressize = 0;
406 }
407
408 enfosuccess = TRUE;
409 leafsol = FALSE;
410
411 /* LP loop; every time a new LP was solved, conditions are checked
412 * dive as long we are in the given objective, depth and iteration limits and fractional variables exist, but
413 * - if possible, we dive at least with the depth 10
414 * - if the number of fractional variables decreased at least with 1 variable per 2 dive depths, we continue diving
415 */
416 while( !lperror && !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL && enfosuccess
417 && (SCIPgetProbingDepth(scip) < 10
418 || nlpcands <= startndivecands - SCIPgetProbingDepth(scip) / 2
419 || (SCIPgetProbingDepth(scip) < maxdivedepth && SCIPdivesetGetNLPIterations(diveset, divecontext) < maxnlpiterations && SCIPgetLPObjval(scip) < searchbound))
420 && !SCIPisStopped(scip) )
421 {
422 SCIP_Real lastlpobjval;
423 int c;
424 SCIP_Bool allroundable;
425 SCIP_Bool infeasible;
426 int prevcandsinsertidx;
427
428 /* remember the last LP depth */
429 assert(lastlpdepth < SCIPgetProbingDepth(scip));
430 lastlpdepth = SCIPgetProbingDepth(scip);
431 domreds = 0;
432
433 SCIPdebugMsg(scip, "%s heuristic continues diving at depth %d, %d candidates left\n",
434 SCIPdivesetGetName(diveset), lastlpdepth, nlpcands);
435
436 /* loop over candidates and determine if they are roundable */
437 allroundable = TRUE;
438 c = 0;
439 while( allroundable && c < nlpcands )
440 {
441 if( SCIPvarMayRoundDown(lpcands[c]) || SCIPvarMayRoundUp(lpcands[c]) )
442 allroundable = TRUE;
443 else
444 allroundable = FALSE;
445 ++c;
446 }
447
448 /* if all candidates are roundable, try to round the solution */
449 if( allroundable )
450 {
451 success = FALSE;
452
453 /* working solution must be linked to LP solution */
454 SCIP_CALL( SCIPlinkLPSol(scip, worksol) );
455 /* create solution from diving LP and try to round it */
456 SCIP_CALL( SCIProundSol(scip, worksol, &success) );
457
458 /* adjust SOS1 constraints */
459 if( success && sos1conshdlr != NULL )
460 {
461 SCIP_Bool changed;
462 SCIP_CALL( SCIPmakeSOS1sFeasible(scip, sos1conshdlr, worksol, &changed, &success) );
463 }
464
465 /* succesfully rounded solutions are tried for primal feasibility */
466 if( success )
467 {
468 SCIP_Bool changed = FALSE;
469 SCIPdebugMsg(scip, "%s found roundable primal solution: obj=%g\n", SCIPdivesetGetName(diveset), SCIPgetSolOrigObj(scip, worksol));
470
471 /* adjust indicator constraints */
472 if( indconshdlr != NULL )
473 {
474 SCIP_CALL( SCIPmakeIndicatorsFeasible(scip, indconshdlr, worksol, &changed) );
475 }
476
477 success = FALSE;
478 /* try to add solution to SCIP */
479 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, changed, &success) );
480
481 /* check, if solution was feasible and good enough */
482 if( success )
483 {
484 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n");
485 *result = SCIP_FOUNDSOL;
486 leafsol = (nlpcands == 0);
487
488 /* the rounded solution found above led to a cutoff of the node LP solution */
489 if( SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OBJLIMIT )
490 {
491 cutoff = TRUE;
492 break;
493 }
494 }
495 }
496 }
497
498 /* working solution must be linked to LP solution */
499 assert(SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL);
500 lastlpobjval = SCIPgetLPObjval(scip);
501 SCIP_CALL( SCIPlinkLPSol(scip, worksol) );
502
503 /* we must explicitly store the solution values by unlinking the solution, otherwise,
504 * the working solution may contain wrong entries, if, e.g., a backtrack occurred after an
505 * intermediate LP resolve or the LP was resolved during conflict analysis.
506 */
507 SCIP_CALL( SCIPunlinkSol(scip, worksol) );
508
509 /* ensure array sizes for the diving on the fractional variables */
510 if( onlylpbranchcands && nlpcands > lpcandsscoressize )
511 {
512 assert(nlpcands > 0);
513 assert(lpcandsscores != NULL);
514 assert(lpcandroundup != NULL);
515
516 SCIP_CALL( SCIPreallocBufferArray(scip, &lpcandsscores, nlpcands) );
517 SCIP_CALL( SCIPreallocBufferArray(scip, &lpcandroundup, nlpcands) );
518
519 lpcandsscoressize = nlpcands;
520 }
521
522 enfosuccess = FALSE;
523 /* select the next diving action by selecting appropriate dive bound changes for the preferred and alternative child */
524 SCIP_CALL( selectNextDiving(scip, diveset, worksol, onlylpbranchcands, SCIPgetProbingDepth(scip) == lastlpdepth,
525 lpcands, lpcandssol, lpcandsfrac, lpcandsscores, lpcandroundup, &nviollpcands, nlpcands,
526 &enfosuccess, &infeasible) );
527
528 /* if we did not succeed finding an enforcement, the solution is potentially feasible and we break immediately */
529 if( ! enfosuccess )
530 break;
531
532 prevcandsinsertidx = -1;
533
534 /* start propagating candidate variables
535 * - until the desired targetdepth is reached,
536 * - or there is no further candidate variable left because of intermediate bound changes,
537 * - or a cutoff is detected
538 */
539 do
540 {
541 SCIP_VAR* bdchgvar;
542 SCIP_Real bdchgvalue;
543 SCIP_Longint localdomreds;
544 SCIP_BRANCHDIR bdchgdir;
545 int nbdchanges;
546
547 /* ensure that a new candidate was successfully determined (usually at the end of the previous loop iteration) */
548 assert(enfosuccess);
549 bdchgvar = NULL;
550 bdchgvalue = SCIP_INVALID;
551 bdchgdir = SCIP_BRANCHDIR_AUTO;
552
553 backtracked = FALSE;
554 do
555 {
556 int d;
557 SCIP_VAR** bdchgvars;
558 SCIP_BRANCHDIR* bdchgdirs;
559 SCIP_Real* bdchgvals;
560
561 nbdchanges = 0;
562 /* get the bound change information stored in the dive set */
563 SCIPgetDiveBoundChangeData(scip, &bdchgvars, &bdchgdirs, &bdchgvals, &nbdchanges, !backtracked);
564
565 assert(nbdchanges > 0);
566 assert(bdchgvars != NULL);
567 assert(bdchgdirs != NULL);
568 assert(bdchgvals != NULL);
569
570 /* return if we reached the depth limit */
571 if( SCIP_MAXTREEDEPTH <= SCIPgetDepth(scip) )
572 {
573 SCIPdebugMsg(scip, "depth limit reached, we stop diving immediately.\n");
574 goto TERMINATE;
575 }
576
577 /* dive deeper into the tree */
578 SCIP_CALL( SCIPnewProbingNode(scip) );
579 ++totalnprobingnodes;
580
581 /* apply all suggested domain changes of the variables */
582 for( d = 0; d < nbdchanges; ++d )
583 {
584 SCIP_Real lblocal;
585 SCIP_Real ublocal;
586 SCIP_Bool infeasbdchange;
587
588 /* get the next bound change data: variable, direction, and value */
589 bdchgvar = bdchgvars[d];
590 bdchgvalue = bdchgvals[d];
591 bdchgdir = bdchgdirs[d];
592
593 assert(bdchgvar != NULL);
594 lblocal = SCIPvarGetLbLocal(bdchgvar);
595 ublocal = SCIPvarGetUbLocal(bdchgvar);
596
597 SCIPdebugMsg(scip, " dive %d/%d, LP iter %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ": var <%s>, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
598 SCIPgetProbingDepth(scip), maxdivedepth, SCIPdivesetGetNLPIterations(diveset, divecontext), maxnlpiterations,
599 SCIPvarGetName(bdchgvar),
600 lblocal, ublocal,
601 bdchgdir == SCIP_BRANCHDIR_DOWNWARDS ? lblocal : bdchgvalue,
602 bdchgdir == SCIP_BRANCHDIR_UPWARDS ? ublocal : bdchgvalue);
603
604 infeasbdchange = FALSE;
605 /* tighten the lower and/or upper bound depending on the bound change type */
606 switch( bdchgdir )
607 {
608 case SCIP_BRANCHDIR_UPWARDS:
609 /* test if bound change is possible, otherwise propagation might have deduced the same
610 * bound already or numerical troubles might have occurred */
611 if( SCIPisFeasLE(scip, bdchgvalue, lblocal) || SCIPisFeasGT(scip, bdchgvalue, ublocal) )
612 infeasbdchange = TRUE;
613 else
614 {
615 /* round variable up */
616 SCIP_CALL( SCIPchgVarLbProbing(scip, bdchgvar, bdchgvalue) );
617 }
618 break;
619 case SCIP_BRANCHDIR_DOWNWARDS:
620 /* test if bound change is possible, otherwise propagation might have deduced the same
621 * bound already or numerical troubles might have occurred */
622 if( SCIPisFeasGE(scip, bdchgvalue, ublocal) || SCIPisFeasLT(scip, bdchgvalue, lblocal) )
623 infeasbdchange = TRUE;
624 else
625 {
626 /* round variable down */
627 SCIP_CALL( SCIPchgVarUbProbing(scip, bdchgvar, bdchgvalue) );
628 }
629 break;
630 case SCIP_BRANCHDIR_FIXED:
631 /* test if bound change is possible, otherwise propagation might have deduced the same
632 * bound already or numerical troubles might have occurred */
633 if( SCIPisFeasLT(scip, bdchgvalue, lblocal) || SCIPisFeasGT(scip, bdchgvalue, ublocal) ||
634 (SCIPisFeasEQ(scip, lblocal, ublocal) && nbdchanges < 2) )
635 infeasbdchange = TRUE;
636 else
637 {
638 /* fix variable to the bound change value */
639 if( SCIPisFeasLT(scip, lblocal, bdchgvalue) )
640 {
641 SCIP_CALL( SCIPchgVarLbProbing(scip, bdchgvar, bdchgvalue) );
642 }
643 if( SCIPisFeasGT(scip, ublocal, bdchgvalue) )
644 {
645 SCIP_CALL( SCIPchgVarUbProbing(scip, bdchgvar, bdchgvalue) );
646 }
647 }
648 break;
649 case SCIP_BRANCHDIR_AUTO:
650 default:
651 SCIPerrorMessage("Error: Unsupported bound change direction <%d> specified for diving, aborting\n",bdchgdirs[d]);
652 SCIPABORT();
653 return SCIP_INVALIDDATA; /*lint !e527*/
654 }
655 /* if the variable domain has been shrunk in the meantime, numerical troubles may have
656 * occured or variable was fixed by propagation while backtracking => Abort diving!
657 */
658 if( infeasbdchange )
659 {
660 SCIPdebugMsg(scip, "\nSelected variable <%s> domain already [%g,%g] as least as tight as desired bound change, diving aborted \n",
661 SCIPvarGetName(bdchgvar), SCIPvarGetLbLocal(bdchgvar), SCIPvarGetUbLocal(bdchgvar));
662 cutoff = TRUE;
663 break;
664 }
665 }
666 /* break loop immediately if we detected a cutoff */
667 if( cutoff )
668 break;
669
670 /* apply domain propagation */
671 localdomreds = 0;
672 SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, &localdomreds) );
673
674 /* add the number of bound changes we applied by ourselves after propagation, otherwise the counter would have been reset */
675 localdomreds += nbdchanges;
676
677 SCIPdebugMsg(scip, "domain reductions: %" SCIP_LONGINT_FORMAT " (total: %" SCIP_LONGINT_FORMAT ")\n",
678 localdomreds, domreds + localdomreds);
679
680 /* resolve the diving LP if the diving resolve frequency is reached or a sufficient number of intermediate bound changes
681 * was reached
682 */
683 if( ! cutoff
684 && ((lpsolvefreq > 0 && ((SCIPgetProbingDepth(scip) - lastlpdepth) % lpsolvefreq) == 0)
685 || (domreds + localdomreds > SCIPdivesetGetLPResolveDomChgQuot(diveset) * SCIPgetNVars(scip))
686 || (onlylpbranchcands && nviollpcands > (int)(SCIPdivesetGetLPResolveDomChgQuot(diveset) * nlpcands))) )
687 {
688 SCIP_CALL( solveLP(scip, diveset, maxnlpiterations, divecontext, &lperror, &cutoff) );
689
690 /* lp errors lead to early termination */
691 if( lperror )
692 {
693 cutoff = TRUE;
694 break;
695 }
696 }
697
698 /* perform backtracking if a cutoff was detected */
699 if( cutoff && !backtracked && SCIPdivesetUseBacktrack(diveset) )
700 {
701 SCIPdebugMsg(scip, " *** cutoff detected at level %d - backtracking\n", SCIPgetProbingDepth(scip));
702 SCIP_CALL( SCIPbacktrackProbing(scip, SCIPgetProbingDepth(scip) - 1) );
703 ++totalnbacktracks;
704 backtracked = TRUE;
705 backtrack = TRUE;
706 cutoff = FALSE;
707 }
708 else
709 backtrack = FALSE;
710 }
711 while( backtrack );
712
713 /* we add the domain reductions from the last evaluated node */
714 domreds += localdomreds; /*lint !e771 lint thinks localdomreds has not been initialized */
715
716 /* store candidate for pseudo cost update and choose next candidate only if no cutoff was detected */
717 if( ! cutoff )
718 {
719 if( nbdchanges == 1 && (bdchgdir == SCIP_BRANCHDIR_UPWARDS || bdchgdir == SCIP_BRANCHDIR_DOWNWARDS) )
720 {
721 ++prevcandsinsertidx;
722 assert(prevcandsinsertidx <= SCIPgetProbingDepth(scip) - lastlpdepth - 1);
723 assert(SCIPgetProbingDepth(scip) > 0);
724 assert(bdchgvar != NULL);
725 assert(bdchgvalue != SCIP_INVALID); /*lint !e777 doesn't like comparing floats for equality */
726
727 /* extend array in case of a dynamic, domain change based LP resolve strategy */
728 if( prevcandsinsertidx >= previouscandssize )
729 {
730 previouscandssize *= 2;
731 SCIP_CALL( SCIPreallocBufferArray(scip, &previouscands, previouscandssize) );
732 SCIP_CALL( SCIPreallocBufferArray(scip, &previousvals, previouscandssize) );
733 }
734 assert(previouscandssize > prevcandsinsertidx);
735
736 /* store candidate for pseudo cost update */
737 previouscands[prevcandsinsertidx] = bdchgvar;
738 previousvals[prevcandsinsertidx] = bdchgvalue;
739 }
740
741 /* choose next candidate variable and resolve the LP if none is found. */
742 if( SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_NOTSOLVED )
743 {
744 assert(SCIPgetProbingDepth(scip) > lastlpdepth);
745 enfosuccess = FALSE;
746
747 /* select the next diving action */
748 SCIP_CALL( selectNextDiving(scip, diveset, worksol, onlylpbranchcands, SCIPgetProbingDepth(scip) == lastlpdepth,
749 lpcands, lpcandssol, lpcandsfrac, lpcandsscores, lpcandroundup, &nviollpcands, nlpcands,
750 &enfosuccess, &infeasible) );
751
752 /* in case of an unsuccesful candidate search, we solve the node LP */
753 if( !enfosuccess )
754 {
755 SCIP_CALL( solveLP(scip, diveset, maxnlpiterations, divecontext, &lperror, &cutoff) );
756
757 /* check for an LP error and terminate in this case, cutoffs lead to termination anyway */
758 if( lperror )
759 cutoff = TRUE;
760
761 /* enfosuccess must be set to TRUE for entering the main LP loop again */
762 enfosuccess = TRUE;
763 }
764 }
765 }
766 }
767 while( !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_NOTSOLVED );
768
769 assert(cutoff || lperror || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_NOTSOLVED);
770
771 assert(cutoff || (SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OBJLIMIT && SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_INFEASIBLE &&
772 (SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL || SCIPisLT(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)))));
773
774 /* check new LP candidates and use the LP Objective gain to update pseudo cost information */
775 if( ! cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
776 {
777 int v;
778 SCIP_Real gain;
779
780 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL) );
781
782 SCIPdebugMsg(scip, " -> lpsolstat=%d, objval=%g/%g, nfrac=%d\n", SCIPgetLPSolstat(scip), SCIPgetLPObjval(scip), searchbound, nlpcands);
783 /* distribute the gain equally over all variables that we rounded since the last LP */
784 gain = SCIPgetLPObjval(scip) - lastlpobjval;
785 gain = MAX(gain, 0.0);
786 gain /= (1.0 * (SCIPgetProbingDepth(scip) - lastlpdepth));
787
788 /* loop over previously fixed candidates and share gain improvement */
789 for( v = 0; v <= prevcandsinsertidx; ++v )
790 {
791 SCIP_VAR* cand = previouscands[v];
792 SCIP_Real val = previousvals[v];
793 SCIP_Real solval = SCIPgetSolVal(scip, worksol, cand);
794
795 /* todo: should the gain be shared with a smaller weight, instead of dividing the gain itself? */
796 /* it may happen that a variable had an integral solution value beforehand, e.g., for indicator variables */
797 if( ! SCIPisZero(scip, val - solval) )
798 {
799 SCIP_CALL( SCIPupdateVarPseudocost(scip, cand, val - solval, gain, 1.0) );
800 }
801 }
802 }
803 else
804 nlpcands = 0;
805 }
806
807 success = FALSE;
808 /* check if a solution has been found */
809 if( !enfosuccess && !lperror && !cutoff && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
810 {
811 /* create solution from diving LP */
812 SCIP_CALL( SCIPlinkLPSol(scip, worksol) );
813 SCIPdebugMsg(scip, "%s found primal solution: obj=%g\n", SCIPdivesetGetName(diveset), SCIPgetSolOrigObj(scip, worksol));
814
815 /* try to add solution to SCIP */
816 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, FALSE, &success) );
817
818 /* check, if solution was feasible and good enough */
819 if( success )
820 {
821 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n");
822 *result = SCIP_FOUNDSOL;
823 leafsol = TRUE;
824 }
825 }
826
827 SCIPupdateDivesetStats(scip, diveset, totalnprobingnodes, totalnbacktracks, SCIPgetNSolsFound(scip) - oldnsolsfound,
828 SCIPgetNBestSolsFound(scip) - oldnbestsolsfound, SCIPgetNConflictConssFound(scip) - oldnconflictsfound, leafsol, divecontext);
829
830 SCIPdebugMsg(scip, "(node %" SCIP_LONGINT_FORMAT ") finished %s diveset (%s heur): %d fractionals, dive %d/%d, LP iter %" SCIP_LONGINT_FORMAT "/%" SCIP_LONGINT_FORMAT ", objval=%g/%g, lpsolstat=%d, cutoff=%u\n",
831 SCIPgetNNodes(scip), SCIPdivesetGetName(diveset), SCIPheurGetName(heur), nlpcands, SCIPgetProbingDepth(scip), maxdivedepth, SCIPdivesetGetNLPIterations(diveset, divecontext), maxnlpiterations,
832 SCIPretransformObj(scip, SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ? SCIPgetLPObjval(scip) : SCIPinfinity(scip)), SCIPretransformObj(scip, searchbound), SCIPgetLPSolstat(scip), cutoff);
833
834 TERMINATE:
835 /* end probing mode */
836 SCIP_CALL( SCIPendProbing(scip) );
837
838 SCIPdivesetSetWorkSolution(diveset, NULL);
839
840 if( onlylpbranchcands )
841 {
842 SCIPfreeBufferArray(scip, &lpcandroundup);
843 SCIPfreeBufferArray(scip, &lpcandsscores);
844 }
845 SCIPfreeBufferArray(scip, &previousvals);
846 SCIPfreeBufferArray(scip, &previouscands);
847
848 return SCIP_OKAY;
849 }
850
851
852 /** creates the rows of the subproblem */
853 static
createRows(SCIP * scip,SCIP * subscip,SCIP_HASHMAP * varmap)854 SCIP_RETCODE createRows(
855 SCIP* scip, /**< original SCIP data structure */
856 SCIP* subscip, /**< SCIP data structure for the subproblem */
857 SCIP_HASHMAP* varmap /**< a hashmap to store the mapping of source variables to the corresponding
858 * target variables */
859 )
860 {
861 SCIP_ROW** rows; /* original scip rows */
862 SCIP_CONS* cons; /* new constraint */
863 SCIP_VAR** consvars; /* new constraint's variables */
864 SCIP_COL** cols; /* original row's columns */
865
866 SCIP_Real constant; /* constant added to the row */
867 SCIP_Real lhs; /* left hand side of the row */
868 SCIP_Real rhs; /* left right side of the row */
869 SCIP_Real* vals; /* variables' coefficient values of the row */
870
871 int nrows;
872 int nnonz;
873 int i;
874 int j;
875
876 /* get the rows and their number */
877 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
878
879 /* copy all rows to linear constraints */
880 for( i = 0; i < nrows; i++ )
881 {
882 /* ignore rows that are only locally valid */
883 if( SCIProwIsLocal(rows[i]) )
884 continue;
885
886 /* get the row's data */
887 constant = SCIProwGetConstant(rows[i]);
888 lhs = SCIProwGetLhs(rows[i]) - constant;
889 rhs = SCIProwGetRhs(rows[i]) - constant;
890 vals = SCIProwGetVals(rows[i]);
891 nnonz = SCIProwGetNNonz(rows[i]);
892 cols = SCIProwGetCols(rows[i]);
893
894 assert(lhs <= rhs);
895
896 /* allocate memory array to be filled with the corresponding subproblem variables */
897 SCIP_CALL( SCIPallocBufferArray(scip, &consvars, nnonz) );
898 for( j = 0; j < nnonz; j++ )
899 consvars[j] = (SCIP_VAR*) SCIPhashmapGetImage(varmap, (SCIPcolGetVar(cols[j])));
900
901 /* create a new linear constraint and add it to the subproblem */
902 SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, SCIProwGetName(rows[i]), nnonz, consvars, vals, lhs, rhs,
903 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
904 SCIP_CALL( SCIPaddCons(subscip, cons) );
905 SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
906
907 /* free temporary memory */
908 SCIPfreeBufferArray(scip, &consvars);
909 }
910
911 return SCIP_OKAY;
912 }
913
914
915 /** get a sub-SCIP copy of the transformed problem */
SCIPcopyLargeNeighborhoodSearch(SCIP * sourcescip,SCIP * subscip,SCIP_HASHMAP * varmap,const char * suffix,SCIP_VAR ** fixedvars,SCIP_Real * fixedvals,int nfixedvars,SCIP_Bool uselprows,SCIP_Bool copycuts,SCIP_Bool * success,SCIP_Bool * valid)916 SCIP_RETCODE SCIPcopyLargeNeighborhoodSearch(
917 SCIP* sourcescip, /**< source SCIP data structure */
918 SCIP* subscip, /**< sub-SCIP used by the heuristic */
919 SCIP_HASHMAP* varmap, /**< a hashmap to store the mapping of source variables to the corresponding
920 * target variables */
921 const char* suffix, /**< suffix for the problem name */
922 SCIP_VAR** fixedvars, /**< source variables whose copies should be fixed in the target SCIP environment, or NULL */
923 SCIP_Real* fixedvals, /**< array of fixing values for target SCIP variables, or NULL */
924 int nfixedvars, /**< number of source variables whose copies should be fixed in the target SCIP environment, or NULL */
925 SCIP_Bool uselprows, /**< should the linear relaxation of the problem defined by LP rows be copied? */
926 SCIP_Bool copycuts, /**< should cuts be copied (only if uselprows == FALSE) */
927 SCIP_Bool* success, /**< was the copying successful? */
928 SCIP_Bool* valid /**< pointer to store whether the copying was valid, or NULL */
929 )
930 {
931 assert(sourcescip != NULL);
932 assert(suffix != NULL);
933 assert(subscip != NULL);
934 assert(varmap != NULL);
935 assert(success != NULL);
936
937 if( uselprows )
938 {
939 char probname[SCIP_MAXSTRLEN];
940
941 /* copy all plugins */
942 SCIP_CALL( SCIPincludeDefaultPlugins(subscip) );
943
944 /* set name to the original problem name and possibly add a suffix */
945 (void) SCIPsnprintf(probname, SCIP_MAXSTRLEN, "%s_%s", SCIPgetProbName(sourcescip), suffix);
946
947 /* create the subproblem */
948 SCIP_CALL( SCIPcreateProb(subscip, probname, NULL, NULL, NULL, NULL, NULL, NULL, NULL) );
949
950 /* copy all variables */
951 SCIP_CALL( SCIPcopyVars(sourcescip, subscip, varmap, NULL, fixedvars, fixedvals, nfixedvars, TRUE) );
952
953 /* copy parameter settings */
954 SCIP_CALL( SCIPcopyParamSettings(sourcescip, subscip) );
955
956 /* create linear constraints from LP rows of the source problem */
957 SCIP_CALL( createRows(sourcescip, subscip, varmap) );
958 }
959 else
960 {
961 SCIP_CALL( SCIPcopyConsCompression(sourcescip, subscip, varmap, NULL, suffix, fixedvars, fixedvals, nfixedvars,
962 TRUE, FALSE, FALSE, TRUE, valid) );
963
964 if( copycuts )
965 {
966 /* copies all active cuts from cutpool of sourcescip to linear constraints in targetscip */
967 SCIP_CALL( SCIPcopyCuts(sourcescip, subscip, varmap, NULL, TRUE, NULL) );
968 }
969 }
970
971 *success = TRUE;
972
973 return SCIP_OKAY;
974 }
975
976 /** adds a trust region neighborhood constraint to the @p targetscip
977 *
978 * a trust region constraint measures the deviation from the current incumbent solution \f$x^*\f$ by an auxiliary
979 * continuous variable \f$v \geq 0\f$:
980 * \f[
981 * \sum\limits_{j\in B} |x_j^* - x_j| = v
982 * \f]
983 * Only binary variables are taken into account. The deviation is penalized in the objective function using
984 * a positive \p violpenalty.
985 *
986 * @note: the trust region constraint creates an auxiliary variable to penalize the deviation from
987 * the current incumbent solution. This variable can afterwards be accessed using SCIPfindVar() by its name
988 * 'trustregion_violationvar'
989 */
SCIPaddTrustregionNeighborhoodConstraint(SCIP * sourcescip,SCIP * targetscip,SCIP_VAR ** subvars,SCIP_Real violpenalty)990 SCIP_RETCODE SCIPaddTrustregionNeighborhoodConstraint(
991 SCIP* sourcescip, /**< the data structure for the main SCIP instance */
992 SCIP* targetscip, /**< SCIP data structure of the subproblem */
993 SCIP_VAR** subvars, /**< variables of the subproblem, NULL entries are ignored */
994 SCIP_Real violpenalty /**< the penalty for violating the trust region */
995 )
996 {
997 SCIP_VAR* violvar;
998 SCIP_CONS* trustregioncons;
999 SCIP_VAR** consvars;
1000 SCIP_VAR** vars;
1001 SCIP_SOL* bestsol;
1002
1003 int nvars;
1004 int nbinvars;
1005 int nconsvars;
1006 int i;
1007 SCIP_Real rhs;
1008 SCIP_Real* consvals;
1009 char name[SCIP_MAXSTRLEN];
1010
1011 /* get the data of the variables and the best solution */
1012 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, &nvars, &nbinvars, NULL, NULL, NULL) );
1013 bestsol = SCIPgetBestSol(sourcescip);
1014 assert(bestsol != NULL);
1015 /* otherwise, this neighborhood would not be active in the first place */
1016 assert(nbinvars > 0);
1017
1018 /* memory allocation */
1019 SCIP_CALL( SCIPallocBufferArray(sourcescip, &consvars, nbinvars + 1) );
1020 SCIP_CALL( SCIPallocBufferArray(sourcescip, &consvals, nbinvars + 1) );
1021 nconsvars = 0;
1022
1023 /* set initial left and right hand sides of trust region constraint */
1024 rhs = 0.0;
1025
1026 /* create the distance (to incumbent) function of the binary variables */
1027 for( i = 0; i < nbinvars; i++ )
1028 {
1029 SCIP_Real solval;
1030
1031 if( subvars[i] == NULL )
1032 continue;
1033
1034 solval = SCIPgetSolVal(sourcescip, bestsol, vars[i]);
1035 assert( SCIPisFeasIntegral(sourcescip,solval) );
1036
1037 /* is variable i part of the binary support of bestsol? */
1038 if( SCIPisFeasEQ(sourcescip, solval, 1.0) )
1039 {
1040 consvals[nconsvars] = -1.0;
1041 rhs -= 1.0;
1042 }
1043 else
1044 consvals[nconsvars] = 1.0;
1045 consvars[nconsvars] = subvars[i];
1046 assert(SCIPvarGetType(consvars[nconsvars]) == SCIP_VARTYPE_BINARY);
1047 ++nconsvars;
1048 }
1049
1050 /* adding the violation variable */
1051 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_trustregionviolvar", SCIPgetProbName(sourcescip));
1052 SCIP_CALL( SCIPcreateVarBasic(targetscip, &violvar, name, 0.0, SCIPinfinity(targetscip), violpenalty, SCIP_VARTYPE_CONTINUOUS) );
1053 SCIP_CALL( SCIPaddVar(targetscip, violvar) );
1054 consvars[nconsvars] = violvar;
1055 consvals[nconsvars] = -1.0;
1056 ++nconsvars;
1057
1058 /* creates trustregion constraint and adds it to subscip */
1059 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "%s_trustregioncons", SCIPgetProbName(sourcescip));
1060
1061 SCIP_CALL( SCIPcreateConsLinear(targetscip, &trustregioncons, name, nconsvars, consvars, consvals,
1062 rhs, rhs, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
1063 SCIP_CALL( SCIPaddCons(targetscip, trustregioncons) );
1064 SCIP_CALL( SCIPreleaseCons(targetscip, &trustregioncons) );
1065
1066 /* releasing the violation variable */
1067 SCIP_CALL( SCIPreleaseVar(targetscip, &violvar) );
1068
1069 /* free local memory */
1070 SCIPfreeBufferArray(sourcescip, &consvals);
1071 SCIPfreeBufferArray(sourcescip, &consvars);
1072
1073 return SCIP_OKAY;
1074 }
1075