1
2 /*
3 Mixed integer programming optimization drivers for lp_solve v5.0+
4 ----------------------------------------------------------------------------------
5 Author: Michel Berkelaar (to lp_solve v3.2)
6 Kjell Eikland (v4.0 and forward)
7 Contact:
8 License terms: LGPL.
9
10 Requires: string.h, float.h, commonlib.h, lp_lib.h, lp_report.h,
11 lp_simplex.h
12
13 Release notes:
14 v5.0.0 31 January 2004 New unit isolating B&B routines.
15 v5.0.1 01 February 2004 Complete rewrite into non-recursive version.
16 v5.0.2 05 April 2004 Expanded pseudocosting with options for MIP fraction
17 counts and "cost/benefit" ratio (KE special!).
18 Added GUB functionality based on SOS structures.
19 v5.0.3 1 May 2004 Changed routine names to be more intuitive.
20 v5.0.4 15 May 2004 Added functinality to pack bounds in order to
21 conserve memory in B&B-processing large MIP models.
22 v5.1.0 25 July 2004 Added functions for dynamic cut generation.
23 v5.2.0 15 December 2004 Added functions for reduced cost variable fixing
24 and converted to delta-model of B&B bound storage.
25 ----------------------------------------------------------------------------------
26 */
27
28 #include <string.h>
29 #include <float.h>
30 #include "commonlib.h"
31 #include "lp_lib.h"
32 #include "lp_scale.h"
33 #include "lp_report.h"
34 #include "lp_simplex.h"
35 #include "lp_mipbb.h"
36
37 #ifdef FORTIFY
38 # include "lp_fortify.h"
39 #endif
40
41
42 /* Allocation routine for the BB record structure */
create_BB(lprec * lp,BBrec * parentBB,MYBOOL dofullcopy)43 STATIC BBrec *create_BB(lprec *lp, BBrec *parentBB, MYBOOL dofullcopy)
44 {
45 BBrec *newBB;
46
47 newBB = (BBrec *) calloc(1, sizeof(*newBB));
48 if(newBB != NULL) {
49
50 if(parentBB == NULL) {
51 allocREAL(lp, &newBB->upbo, lp->sum + 1, FALSE);
52 allocREAL(lp, &newBB->lowbo, lp->sum + 1, FALSE);
53 MEMCOPY(newBB->upbo, lp->orig_upbo, lp->sum + 1);
54 MEMCOPY(newBB->lowbo, lp->orig_lowbo, lp->sum + 1);
55 }
56 else if(dofullcopy) {
57 allocREAL(lp, &newBB->upbo, lp->sum + 1, FALSE);
58 allocREAL(lp, &newBB->lowbo, lp->sum + 1, FALSE);
59 MEMCOPY(newBB->upbo, parentBB->upbo, lp->sum + 1);
60 MEMCOPY(newBB->lowbo, parentBB->lowbo, lp->sum + 1);
61 }
62 else {
63 newBB->upbo = parentBB->upbo;
64 newBB->lowbo = parentBB->lowbo;
65 }
66 newBB->contentmode = dofullcopy;
67
68 newBB->lp = lp;
69
70 /* Set parent by default, but not child */
71 newBB->parent = parentBB;
72
73 }
74 return( newBB );
75 }
76
77
78 /* Pushing and popping routines for the B&B structure */
79
push_BB(lprec * lp,BBrec * parentBB,int varno,int vartype,int varcus)80 STATIC BBrec *push_BB(lprec *lp, BBrec *parentBB, int varno, int vartype, int varcus)
81 /* Push ingoing bounds and B&B data onto the stack */
82 {
83 BBrec *newBB;
84
85 /* Do initialization and updates */
86 if(parentBB == NULL)
87 parentBB = lp->bb_bounds;
88 newBB = create_BB(lp, parentBB, FALSE);
89 if(newBB != NULL) {
90
91 newBB->varno = varno;
92 newBB->vartype = vartype;
93 newBB->lastvarcus = varcus;
94 incrementUndoLadder(lp->bb_lowerchange);
95 newBB->LBtrack++;
96 incrementUndoLadder(lp->bb_upperchange);
97 newBB->UBtrack++;
98
99 /* Adjust variable fixing/bound tightening based on the last reduced cost */
100 if((parentBB != NULL) && (parentBB->lastrcf > 0)) {
101 MYBOOL isINT;
102 int k, ii, nfixed = 0, ntighten = 0;
103 REAL deltaUL;
104
105 for(k = 1; k <= lp->nzdrow[0]; k++) {
106 ii = lp->nzdrow[k];
107 #ifdef UseMilpSlacksRCF /* Check if we should include ranged constraints */
108 isINT = FALSE;
109 #else
110 if(ii <= lp->rows)
111 continue;
112 isINT = is_int(lp, ii-lp->rows);
113 #endif
114 #ifndef UseMilpExpandedRCF /* Don't include non-integers if it is not defined */
115 if(!isINT)
116 continue;
117 #endif
118 switch(abs(rcfbound_BB(newBB, ii, isINT, &deltaUL, NULL))) {
119 case LE: SETMIN(deltaUL, newBB->upbo[ii]);
120 SETMAX(deltaUL, newBB->lowbo[ii]);
121 modifyUndoLadder(lp->bb_upperchange, ii, newBB->upbo, deltaUL);
122 break;
123 case GE: SETMAX(deltaUL, newBB->lowbo[ii]);
124 SETMIN(deltaUL, newBB->upbo[ii]);
125 modifyUndoLadder(lp->bb_lowerchange, ii, newBB->lowbo, deltaUL);
126 break;
127 default: continue;
128 }
129 if(newBB->upbo[ii] == newBB->lowbo[ii])
130 nfixed++;
131 else
132 ntighten++;
133 }
134 if(lp->bb_trace) {
135 report(lp, DETAILED,
136 "push_BB: Used reduced cost to fix %d variables and tighten %d bounds\n",
137 nfixed, ntighten);
138 }
139 }
140
141 /* Handle case where we are pushing at the end */
142 if(parentBB == lp->bb_bounds)
143 lp->bb_bounds = newBB;
144 /* Handle case where we are pushing in the middle */
145 else
146 newBB->child = parentBB->child;
147 if(parentBB != NULL)
148 parentBB->child = newBB;
149
150 lp->bb_level++;
151 if(lp->bb_level > lp->bb_maxlevel)
152 lp->bb_maxlevel = lp->bb_level;
153
154 if(!initbranches_BB(newBB))
155 newBB = pop_BB(newBB);
156 else if(MIP_count(lp) > 0) {
157 if( (lp->bb_level <= 1) && (lp->bb_varactive == NULL) &&
158 (!allocINT(lp, &lp->bb_varactive, lp->columns+1, TRUE) ||
159 !initcuts_BB(lp)) )
160 newBB = pop_BB(newBB);
161 if(varno > 0) {
162 lp->bb_varactive[varno-lp->rows]++;
163 }
164 }
165 }
166 return( newBB );
167 }
168
free_BB(BBrec ** BB)169 STATIC MYBOOL free_BB(BBrec **BB)
170 {
171 MYBOOL parentreturned = FALSE;
172
173 if((BB != NULL) && (*BB != NULL)) {
174 BBrec *parent = (*BB)->parent;
175
176 if((parent == NULL) || (*BB)->contentmode) {
177 FREE((*BB)->upbo);
178 FREE((*BB)->lowbo);
179 }
180 FREE((*BB)->varmanaged);
181 FREE(*BB);
182
183 parentreturned = (MYBOOL) (parent != NULL);
184 if(parentreturned)
185 *BB = parent;
186
187 }
188 return( parentreturned );
189 }
190
pop_BB(BBrec * BB)191 STATIC BBrec *pop_BB(BBrec *BB)
192 /* Pop / free the previously "pushed" / saved bounds */
193 {
194 int k;
195 BBrec *parentBB;
196 lprec *lp = BB->lp;
197
198 if(BB == NULL)
199 return( BB );
200
201 /* Handle case where we are popping the end of the chain */
202 parentBB = BB->parent;
203 if(BB == lp->bb_bounds) {
204 lp->bb_bounds = parentBB;
205 if(parentBB != NULL)
206 parentBB->child = NULL;
207 }
208 /* Handle case where we are popping inside or at the beginning of the chain */
209 else {
210 if(parentBB != NULL)
211 parentBB->child = BB->child;
212 if(BB->child != NULL)
213 BB->child->parent = parentBB;
214 }
215
216 /* Unwind other variables */
217 if(lp->bb_upperchange != NULL) {
218 restoreUndoLadder(lp->bb_upperchange, BB->upbo);
219 for(; BB->UBtrack > 0; BB->UBtrack--) {
220 decrementUndoLadder(lp->bb_upperchange);
221 restoreUndoLadder(lp->bb_upperchange, BB->upbo);
222 }
223 }
224 if(lp->bb_lowerchange != NULL) {
225 restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
226 for(; BB->LBtrack > 0; BB->LBtrack--) {
227 decrementUndoLadder(lp->bb_lowerchange);
228 restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
229 }
230 }
231 lp->bb_level--;
232 k = BB->varno - lp->rows;
233 if(lp->bb_level == 0) {
234 if(lp->bb_varactive != NULL) {
235 FREE(lp->bb_varactive);
236 freecuts_BB(lp);
237 }
238 if(lp->int_vars+lp->sc_vars > 0)
239 free_pseudocost(lp);
240 pop_basis(lp, FALSE);
241 lp->rootbounds = NULL;
242 }
243 else
244 lp->bb_varactive[k]--;
245
246 /* Undo SOS/GUB markers */
247 if(BB->isSOS && (BB->vartype != BB_INT))
248 SOS_unmark(lp->SOS, 0, k);
249 else if(BB->isGUB)
250 SOS_unmark(lp->GUB, 0, k);
251
252 /* Undo the SC marker */
253 if(BB->sc_canset)
254 lp->sc_lobound[k] *= -1;
255
256 /* Pop the associated basis */
257 #if 1
258 /* Original version that does not restore previous basis */
259 pop_basis(lp, FALSE);
260 #else
261 /* Experimental version that restores previous basis */
262 pop_basis(lp, BB->isSOS);
263 #endif
264
265 /* Finally free the B&B object */
266 free_BB(&BB);
267
268 /* Return the parent BB */
269 return( parentBB );
270 }
271
272 /* Here are heuristic routines to see if we need bother with branching further
273
274 1. A probing routine to see of the best OF can be better than incumbent
275 2. A presolve routine to fix other variables and detect infeasibility
276
277 THIS IS INACTIVE CODE, PLACEHOLDERS FOR FUTURE DEVELOPMENT!!! */
probe_BB(BBrec * BB)278 STATIC REAL probe_BB(BBrec *BB)
279 {
280 int i, ii;
281 REAL coefOF, sum = 0;
282 lprec *lp = BB->lp;
283
284 /* Loop over all ints to see if the best possible solution
285 stands any chance of being better than the incumbent solution */
286 if(lp->solutioncount == 0)
287 return( lp->infinite );
288 for(i = 1; i <= lp->columns; i++) {
289 if(!is_int(lp, i))
290 continue;
291 ii = lp->rows + i;
292 coefOF = lp->obj[i];
293 if(coefOF < 0) {
294 if(is_infinite(lp, BB->lowbo[ii]))
295 return( lp->infinite );
296 sum += coefOF * (lp->solution[ii]-BB->lowbo[ii]);
297 }
298 else {
299 if(is_infinite(lp, BB->upbo[ii]))
300 return( lp->infinite );
301 sum += coefOF * (BB->upbo[ii] - lp->solution[ii]);
302 }
303 }
304 return( sum );
305 }
306
presolve_BB(BBrec * BB)307 STATIC REAL presolve_BB(BBrec *BB)
308 {
309 return( 0 );
310 }
311
312 /* Node and branch management routines */
initbranches_BB(BBrec * BB)313 STATIC MYBOOL initbranches_BB(BBrec *BB)
314 {
315 REAL new_bound, temp;
316 int k;
317 lprec *lp = BB->lp;
318
319 /* Create and initialize local bounds and basis */
320 BB->nodestatus = NOTRUN;
321 BB->noderesult = lp->infinite;
322 push_basis(lp, NULL, NULL, NULL);
323
324 /* Set default number of branches at the current B&B branch */
325 if(BB->vartype == BB_REAL)
326 BB->nodesleft = 1;
327
328 else {
329 /* The default is a binary up-low branching */
330 BB->nodesleft = 2;
331
332 /* Initialize the MIP status code pair and set reference values */
333 k = BB->varno - lp->rows;
334 BB->lastsolution = lp->solution[BB->varno];
335
336 /* Determine if we must process in the B&B SOS mode */
337 BB->isSOS = (MYBOOL) ((BB->vartype == BB_SOS) || SOS_is_member(lp->SOS, 0, k));
338 #ifdef Paranoia
339 if((BB->vartype == BB_SOS) && !SOS_is_member(lp->SOS, 0, k))
340 report(lp, SEVERE, "initbranches_BB: Inconsistent identification of SOS variable %s (%d)\n",
341 get_col_name(lp, k), k);
342 #endif
343
344 /* Check if we have a GUB-member variable that needs a triple-branch */
345 BB->isGUB = (MYBOOL) ((BB->vartype == BB_INT) && SOS_can_activate(lp->GUB, 0, k));
346 if(BB->isGUB) {
347 /* Obtain variable index list from applicable GUB - now the first GUB is used,
348 but we could also consider selecting the longest */
349 BB->varmanaged = SOS_get_candidates(lp->GUB, -1, k, TRUE, BB->upbo, BB->lowbo);
350 BB->nodesleft++;
351 }
352
353
354 /* Set local pruning info, automatic, or user-defined strategy */
355 if(BB->vartype == BB_SOS) {
356 if(!SOS_can_activate(lp->SOS, 0, k)) {
357 BB->nodesleft--;
358 BB->isfloor = TRUE;
359 }
360 else
361 BB->isfloor = (MYBOOL) (BB->lastsolution == 0);
362 }
363
364 /* First check if the user wishes to select the branching direction */
365 else if(lp->bb_usebranch != NULL)
366 BB->isfloor = (MYBOOL) lp->bb_usebranch(lp, lp->bb_branchhandle, k);
367
368 /* Otherwise check if we should do automatic branching */
369 else if(get_var_branch(lp, k) == BRANCH_AUTOMATIC) {
370 new_bound = modf(BB->lastsolution/get_pseudorange(lp->bb_PseudoCost, k, BB->vartype), &temp);
371 if(isnan(new_bound))
372 new_bound = 0;
373 else if(new_bound < 0)
374 new_bound += 1.0;
375 BB->isfloor = (MYBOOL) (new_bound <= 0.5);
376
377 /* Set direction by OF value; note that a zero-value in
378 the OF gives priority to floor_first = TRUE */
379 if(is_bb_mode(lp, NODE_GREEDYMODE)) {
380 if(is_bb_mode(lp, NODE_PSEUDOCOSTMODE))
381 BB->sc_bound = get_pseudonodecost(lp->bb_PseudoCost, k, BB->vartype, BB->lastsolution);
382 else
383 BB->sc_bound = mat_getitem(lp->matA, 0, k);
384 new_bound -= 0.5;
385 BB->sc_bound *= new_bound;
386 BB->isfloor = (MYBOOL) (BB->sc_bound > 0);
387 }
388 /* Set direction by pseudocost (normally used in tandem with NODE_PSEUDOxxxSELECT) */
389 else if(is_bb_mode(lp, NODE_PSEUDOCOSTMODE)) {
390 BB->isfloor = (MYBOOL) (get_pseudobranchcost(lp->bb_PseudoCost, k, TRUE) >
391 get_pseudobranchcost(lp->bb_PseudoCost, k, FALSE));
392 if(is_maxim(lp))
393 BB->isfloor = !BB->isfloor;
394 }
395
396 /* Check for reversal */
397 if(is_bb_mode(lp, NODE_BRANCHREVERSEMODE))
398 BB->isfloor = !BB->isfloor;
399 }
400 else
401 BB->isfloor = (MYBOOL) (get_var_branch(lp, k) == BRANCH_FLOOR);
402
403 /* SC logic: If the current SC variable value is in the [0..NZLOBOUND> range, then
404
405 UP: Set lower bound to NZLOBOUND, upper bound is the original
406 LO: Fix the variable by setting upper/lower bound to zero
407
408 ... indicate that the variable is B&B-active by reversing sign of sc_lobound[]. */
409 new_bound = fabs(lp->sc_lobound[k]);
410 BB->sc_bound = new_bound;
411 BB->sc_canset = (MYBOOL) (new_bound != 0);
412
413 /* Must make sure that we handle fractional lower bounds properly;
414 also to ensure that we do a full binary tree search */
415 new_bound = unscaled_value(lp, new_bound, BB->varno);
416 if(is_int(lp, k) && ((new_bound > 0) &&
417 (BB->lastsolution > floor(new_bound)))) {
418 if(BB->lastsolution < ceil(new_bound))
419 BB->lastsolution += 1;
420 modifyUndoLadder(lp->bb_lowerchange, BB->varno, BB->lowbo,
421 scaled_floor(lp, BB->varno, BB->lastsolution, 1));
422 }
423 }
424
425 /* Now initialize the brances and set to first */
426 return( fillbranches_BB(BB) );
427 }
428
fillbranches_BB(BBrec * BB)429 STATIC MYBOOL fillbranches_BB(BBrec *BB)
430 {
431 int K, k;
432 REAL ult_upbo, ult_lowbo;
433 REAL new_bound, SC_bound, intmargin = BB->lp->epsprimal;
434 lprec *lp = BB->lp;
435 MYBOOL OKstatus = FALSE;
436
437 if(lp->bb_break || userabort(lp, MSG_MILPSTRATEGY))
438 return( OKstatus );
439
440 K = BB->varno;
441 if(K > 0) {
442
443 /* Shortcut variables */
444 k = BB->varno - lp->rows;
445 ult_upbo = lp->orig_upbo[K];
446 ult_lowbo = lp->orig_lowbo[K];
447 SC_bound = unscaled_value(lp, BB->sc_bound, K);
448
449 /* First, establish the upper bound to be applied (when isfloor == TRUE)
450 --------------------------------------------------------------------- */
451 /*SetUB:*/
452 BB->UPbound = lp->infinite;
453
454 /* Handle SC-variables for the [0-LoBound> range */
455 if((SC_bound > 0) && (fabs(BB->lastsolution) < SC_bound-intmargin)) {
456 new_bound = 0;
457 }
458 /* Handle pure integers (non-SOS, non-SC) */
459 else if(BB->vartype == BB_INT) {
460 #if 1
461 if(((ult_lowbo >= 0) &&
462 ((floor(BB->lastsolution) < /* Skip cases where the lower bound becomes violated */
463 unscaled_value(lp, MAX(ult_lowbo, fabs(lp->sc_lobound[k])), K)-intmargin))) ||
464 ((ult_upbo <= 0) && /* Was ((ult_lowbo < 0) && */
465 ((floor(BB->lastsolution) > /* Skip cases where the upper bound becomes violated */
466 unscaled_value(lp, MIN(ult_upbo, -fabs(lp->sc_lobound[k])), K)-intmargin)))) {
467 #else
468 if((floor(BB->lastsolution) < /* Skip cases where the lower bound becomes violated */
469 unscaled_value(lp, MAX(ult_lowbo, fabs(lp->sc_lobound[k])), K)-intmargin)) {
470 #endif
471 BB->nodesleft--;
472 goto SetLB;
473 }
474 new_bound = scaled_floor(lp, K, BB->lastsolution, 1);
475 }
476 else if(BB->isSOS) { /* Handle all SOS variants */
477 new_bound = ult_lowbo;
478 if(is_int(lp, k))
479 new_bound = scaled_ceil(lp, K, unscaled_value(lp, new_bound, K), -1);
480 }
481 else /* Handle all other variable incarnations */
482 new_bound = BB->sc_bound;
483
484 /* Check if the new bound might conflict and possibly make adjustments */
485 if(new_bound < BB->lowbo[K])
486 new_bound = BB->lowbo[K] - my_avoidtiny(new_bound-BB->lowbo[K], intmargin);
487 if(new_bound < BB->lowbo[K]) {
488 #ifdef Paranoia
489 debug_print(lp,
490 "fillbranches_BB: New upper bound value %g conflicts with old lower bound %g\n",
491 new_bound, BB->lowbo[K]);
492 #endif
493 BB->nodesleft--;
494 goto SetLB;
495 }
496 #ifdef Paranoia
497 /* Do additional consistency checking */
498 else if(!check_if_less(lp, new_bound, BB->upbo[K], K)) {
499 BB->nodesleft--;
500 goto SetLB;
501 }
502 #endif
503 /* Bound (at least near) feasible */
504 else {
505 /* Makes a difference with models like QUEEN
506 (note consistent use of epsint for scaled integer variables) */
507 if(fabs(new_bound - BB->lowbo[K]) < intmargin*SCALEDINTFIXRANGE)
508 new_bound = BB->lowbo[K];
509 }
510
511 BB->UPbound = new_bound;
512
513
514 /* Next, establish the lower bound to be applied (when isfloor == FALSE)
515 --------------------------------------------------------------------- */
516 SetLB:
517 BB->LObound = -lp->infinite;
518
519 /* Handle SC-variables for the [0-LoBound> range */
520 if((SC_bound > 0) && (fabs(BB->lastsolution) < SC_bound)) {
521 if(is_int(lp, k))
522 new_bound = scaled_ceil(lp, K, SC_bound, 1);
523 else
524 new_bound = BB->sc_bound;
525 }
526 /* Handle pure integers (non-SOS, non-SC, but Ok for GUB!) */
527 else if((BB->vartype == BB_INT)) {
528 if(((ceil(BB->lastsolution) == BB->lastsolution)) || /* Skip branch 0 if the current solution is integer */
529 (ceil(BB->lastsolution) > /* Skip cases where the upper bound becomes violated */
530 unscaled_value(lp, ult_upbo, K)+intmargin) ||
531 (BB->isSOS && (BB->lastsolution == 0))) { /* Don't branch 0 since this is handled in SOS logic */
532 BB->nodesleft--;
533 goto Finish;
534 }
535 new_bound = scaled_ceil(lp, K, BB->lastsolution, 1);
536 }
537 else if(BB->isSOS) { /* Handle all SOS variants */
538 if(SOS_is_member_of_type(lp->SOS, k, SOS3))
539 new_bound = scaled_floor(lp, K, 1, 1);
540 else {
541 new_bound = ult_lowbo;
542 if(is_int(lp, k))
543 new_bound = scaled_floor(lp, K, unscaled_value(lp, new_bound, K), 1);
544 /* If we have a high-order SOS (SOS3+) and this variable is "intermediate"
545 between members previously lower-bounded at a non-zero level, then we should
546 set this and similar neighbouring variables at non-zero lowbo-values (remember
547 that SOS3+ members are all either integers or semi-continuous). Flag this
548 situation and prune tree, since we cannot lower-bound. */
549 if((lp->SOS->maxorder > 2) && (BB->lastsolution == 0) &&
550 SOS_is_member_of_type(lp->SOS, k, SOSn)) {
551 BB->isSOS = AUTOMATIC;
552 }
553 }
554 }
555 else /* Handle all other variable incarnations */
556 new_bound = BB->sc_bound;
557
558 /* Check if the new bound might conflict and possibly make adjustments */
559 if(new_bound > BB->upbo[K])
560 new_bound = BB->upbo[K] + my_avoidtiny(new_bound-BB->upbo[K], intmargin);
561 if(new_bound > BB->upbo[K]) {
562 #ifdef Paranoia
563 debug_print(lp,
564 "fillbranches_BB: New lower bound value %g conflicts with old upper bound %g\n",
565 new_bound, BB->upbo[K]);
566 #endif
567 BB->nodesleft--;
568 goto Finish;
569 }
570 #ifdef Paranoia
571 /* Do additional consistency checking */
572 else if(!check_if_less(lp, BB->lowbo[K], new_bound, K)) {
573 BB->nodesleft--;
574 goto Finish;
575 }
576 #endif
577 /* Bound (at least near-)feasible */
578 else {
579 /* Makes a difference with models like QUEEN
580 (note consistent use of lp->epsprimal for scaled integer variables) */
581 if(fabs(BB->upbo[K]-new_bound) < intmargin*SCALEDINTFIXRANGE)
582 new_bound = BB->upbo[K];
583 }
584
585 BB->LObound = new_bound;
586
587 /* Prepare for the first branch by making sure we are pointing correctly */
588 Finish:
589 if(BB->nodesleft > 0) {
590
591 /* Make sure the change tracker levels are "clean" for the B&B */
592 if(countsUndoLadder(lp->bb_upperchange) > 0) {
593 incrementUndoLadder(lp->bb_upperchange);
594 BB->UBtrack++;
595 }
596 if(countsUndoLadder(lp->bb_lowerchange) > 0) {
597 incrementUndoLadder(lp->bb_lowerchange);
598 BB->LBtrack++;
599 }
600
601 /* Do adjustments */
602 if((BB->vartype != BB_SOS) && (fabs(BB->LObound-BB->UPbound) < intmargin)) {
603 BB->nodesleft--;
604 if(fabs(BB->lowbo[K]-BB->LObound) < intmargin)
605 BB->isfloor = FALSE;
606 else if(fabs(BB->upbo[K]-BB->UPbound) < intmargin)
607 BB->isfloor = TRUE;
608 else {
609 BB->isfloor = TRUE;
610 report(BB->lp, IMPORTANT, "fillbranches_BB: Inconsistent equal-valued bounds for %s\n",
611 get_col_name(BB->lp, k));
612 }
613 }
614 if((BB->nodesleft == 1) &&
615 ((BB->isfloor && (BB->UPbound >= lp->infinite)) ||
616 (!BB->isfloor && (BB->LObound <= -lp->infinite))))
617 BB->isfloor = !BB->isfloor;
618 /* Header initialization */
619 BB->isfloor = !BB->isfloor;
620 while(!OKstatus && /* !userabort(lp, -1) */ lp->spx_status != TIMEOUT && !lp->bb_break && (BB->nodesleft > 0))
621 OKstatus = nextbranch_BB( BB );
622 }
623
624 /* Set an SC variable active, if necessary */
625 if(BB->sc_canset)
626 lp->sc_lobound[k] *= -1;
627
628 }
629 else {
630 BB->nodesleft--;
631 OKstatus = TRUE;
632 }
633
634 return( OKstatus );
635 }
636
637 STATIC MYBOOL nextbranch_BB(BBrec *BB)
638 {
639 int k;
640 lprec *lp = BB->lp;
641 MYBOOL OKstatus = FALSE;
642
643 /* Undo the most recently imposed B&B bounds using the data
644 in the last level of change tracker; this code handles changes
645 to both upper and lower bounds */
646 if(BB->nodessolved > 0) {
647 restoreUndoLadder(lp->bb_upperchange, BB->upbo);
648 restoreUndoLadder(lp->bb_lowerchange, BB->lowbo);
649 }
650
651 if(lp->bb_break || userabort(lp, MSG_MILPSTRATEGY)) {
652 /* Handle the special case of B&B restart;
653 (typically used with the restart after pseudocost initialization) */
654 if((lp->bb_level == 1) && (lp->bb_break == AUTOMATIC)) {
655 lp->bb_break = FALSE;
656 OKstatus = TRUE;
657 }
658 return( OKstatus );
659 }
660
661 if(BB->nodesleft > 0) {
662
663 /* Step and update remaining branch count */
664 k = BB->varno - lp->rows;
665 BB->isfloor = !BB->isfloor;
666 BB->nodesleft--;
667
668 /* Special SOS handling:
669 1) Undo and set new marker for k,
670 2) In case that previous branch was ceiling restore upper bounds of the
671 non-k variables outside of the SOS window set to 0 */
672 if(BB->isSOS && (BB->vartype != BB_INT)) {
673
674 /* First undo previous marker */
675 if((BB->nodessolved > 0) || ((BB->nodessolved == 0) && (BB->nodesleft == 0))) {
676 if(BB->isfloor) {
677 if((BB->nodesleft == 0) && (lp->orig_lowbo[BB->varno] != 0))
678 return( OKstatus );
679 }
680 SOS_unmark(lp->SOS, 0, k);
681 }
682
683 /* Set new SOS marker */
684 if(BB->isfloor) {
685 SOS_set_marked(lp->SOS, 0, k, (MYBOOL) (BB->UPbound != 0));
686 /* Do case of high-order SOS where intervening variables need to be set */
687 if(BB->isSOS == AUTOMATIC) {
688
689 /* SOS_fix_list(lp->SOS, 0, k, BB->lowbo, NULL, AUTOMATIC, lp->bb_lowerchange); */
690 }
691 }
692 else {
693 SOS_set_marked(lp->SOS, 0, k, TRUE);
694 if(SOS_fix_unmarked(lp->SOS, 0, k, BB->upbo, 0, TRUE,
695 NULL, lp->bb_upperchange) < 0)
696 return( OKstatus );
697 }
698 }
699
700 /* Special GUB handling (three branches):
701 1) Undo and set new marker for k,
702 2) Restore upper bounds of the left/right/all non-k variables
703 set to 0 in the previous branch
704 3) Set new upper bounds for the non-k variables (k is set later) */
705 else if(BB->isGUB) {
706
707 /* First undo previous marker */
708 if(BB->nodessolved > 0)
709 SOS_unmark(lp->GUB, 0, k);
710
711 /* Make sure we take floor bound twice */
712 if((BB->nodesleft == 0) && !BB->isfloor)
713 BB->isfloor = !BB->isfloor;
714
715 /* Handle two floor instances;
716 (selected variable and left/right halves of non-selected variables at 0) */
717 SOS_set_marked(lp->GUB, 0, k, (MYBOOL) !BB->isfloor);
718 if(BB->isfloor) {
719 if(SOS_fix_list(lp->GUB, 0, k, BB->upbo,
720 BB->varmanaged, (MYBOOL) (BB->nodesleft > 0), lp->bb_upperchange) < 0)
721 return( OKstatus );
722 }
723 /* Handle one ceil instance;
724 (selected variable at 1, all other at 0) */
725 else {
726 if(SOS_fix_unmarked(lp->GUB, 0, k, BB->upbo, 0, TRUE,
727 NULL, lp->bb_upperchange) < 0)
728 return( OKstatus );
729 }
730 }
731
732 OKstatus = TRUE;
733
734 }
735 /* Initialize simplex status variables */
736 if(OKstatus) {
737 lp->bb_totalnodes++;
738 BB->nodestatus = NOTRUN;
739 BB->noderesult = lp->infinite;
740 }
741 return( OKstatus );
742 }
743
744
745 /* Cut generation and management routines */
746 STATIC MYBOOL initcuts_BB(lprec *lp)
747 {
748 return( TRUE );
749 }
750
751 STATIC int updatecuts_BB(lprec *lp)
752 {
753 return( 0 );
754 }
755
756 STATIC MYBOOL freecuts_BB(lprec *lp)
757 {
758 if(lp->bb_cuttype != NULL)
759 FREE(lp->bb_cuttype);
760 return( TRUE );
761 }
762
763 /* B&B solver routines */
764 STATIC int solve_LP(lprec *lp, BBrec *BB)
765 {
766 int tilted, restored, status;
767 REAL testOF, *upbo = BB->upbo, *lowbo = BB->lowbo;
768 BBrec *perturbed = NULL;
769
770 if(lp->bb_break)
771 return(PROCBREAK);
772
773 #ifdef Paranoia
774 debug_print(lp, "solve_LP: Starting solve for iter %.0f, B&B node level %d.\n",
775 (double) lp->total_iter, lp->bb_level);
776 if(lp->bb_trace &&
777 !validate_bounds(lp, upbo, lowbo))
778 report(lp, SEVERE, "solve_LP: Inconsistent bounds at iter %.0f, B&B node level %d.\n",
779 (double) lp->total_iter, lp->bb_level);
780 #endif
781
782 /* Copy user-specified entering bounds into lp_solve working bounds */
783 impose_bounds(lp, upbo, lowbo);
784
785 /* Restore previously pushed / saved basis for this level if we are in
786 the B&B mode and it is not the first call of the binary tree */
787 if(BB->nodessolved > 1)
788 restore_basis(lp);
789
790 /* Solve and possibly handle degeneracy cases via bound relaxations */
791 status = RUNNING;
792 tilted = 0;
793 restored = 0;
794
795 while(status == RUNNING) {
796
797 /* Copy user-specified entering bounds into lp_solve working bounds and run */
798 status = spx_run(lp, (MYBOOL) (tilted+restored > 0));
799 lp->bb_status = status;
800 lp->spx_perturbed = FALSE;
801
802 if(tilted < 0)
803 break;
804
805 else if((status == OPTIMAL) && (tilted > 0)) {
806 if(lp->spx_trace)
807 report(lp, DETAILED, "solve_LP: Restoring relaxed bounds at level %d.\n",
808 tilted);
809
810 /* Restore original pre-perturbed problem bounds, and solve again using the basis
811 found for the perturbed problem; also make sure we rebase and recompute. */
812 free_BB(&perturbed);
813 if((perturbed == NULL) || (perturbed == BB)) {
814 perturbed = NULL;
815 impose_bounds(lp, upbo, lowbo);
816 }
817 else
818 impose_bounds(lp, perturbed->upbo, perturbed->lowbo);
819 set_action(&lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE);
820 BB->UBzerobased = FALSE;
821 if(lp->bb_totalnodes == 0)
822 lp->real_solution = lp->infinite;
823 status = RUNNING;
824 tilted--;
825 restored++;
826 lp->spx_perturbed = TRUE;
827 }
828
829 else if(((lp->bb_level <= 1) || is_anti_degen(lp, ANTIDEGEN_DURINGBB)) &&
830 (((status == LOSTFEAS) && is_anti_degen(lp, ANTIDEGEN_LOSTFEAS)) ||
831 ((status == INFEASIBLE) && is_anti_degen(lp, ANTIDEGEN_INFEASIBLE)) ||
832 ((status == NUMFAILURE) && is_anti_degen(lp, ANTIDEGEN_NUMFAILURE)) ||
833 ((status == DEGENERATE) && is_anti_degen(lp, ANTIDEGEN_STALLING)))) {
834 /* Allow up to .. consecutive relaxations for non-B&B phases */
835 if((tilted <= DEF_MAXRELAX) && /* Conventional recovery case,... */
836 !((tilted == 0) && (restored > DEF_MAXRELAX))) { /* but not iterating infeasibility */
837
838 /* Create working copy of ingoing bounds if this is the first perturbation */
839 if(tilted == 0)
840 perturbed = BB;
841 perturbed = create_BB(lp, perturbed, TRUE);
842
843 /* Perturb/shift variable bounds; also make sure we rebase and recompute
844 (no refactorization is necessary, since the basis is unchanged) */
845 #if 1
846 perturb_bounds(lp, perturbed, TRUE, TRUE, TRUE);
847 #else
848 perturb_bounds(lp, perturbed, TRUE, TRUE, FALSE);
849 #endif
850 impose_bounds(lp, perturbed->upbo, perturbed->lowbo);
851 set_action(&lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE);
852 BB->UBzerobased = FALSE;
853 status = RUNNING;
854 tilted++;
855 lp->perturb_count++;
856 lp->spx_perturbed = TRUE;
857 if(lp->spx_trace)
858 report(lp, DETAILED, "solve_LP: Starting bound relaxation #%d ('%s')\n",
859 tilted, get_statustext(lp, status));
860 }
861 else {
862 if(lp->spx_trace)
863 report(lp, DETAILED, "solve_LP: Relaxation limit exceeded in resolving infeasibility\n");
864 while((perturbed != NULL) && (perturbed != BB))
865 free_BB(&perturbed);
866 perturbed = NULL;
867 }
868 }
869 }
870
871 /* Handle the different simplex outcomes */
872 if(status != OPTIMAL) {
873 if(lp->bb_level <= 1)
874 lp->bb_parentOF = lp->infinite;
875 if((status == USERABORT) || (status == TIMEOUT)) {
876 /* Construct the last feasible solution, if available */
877 if((lp->solutioncount == 0) &&
878 /*
879 30/01/08 <peno> added MIP_count test because in following situation thing were wrong:
880 - The model contains integers
881 - A break at first is set
882 - A timeout is set
883 - The timeout occurs before a first integer solution is found
884 - When the timeout occurs, the simplex algorithm is in phase 2 and has a feasible (but non-integer) solution, but not optimal yet.
885 If above situation occurs then a (sub-optimal) solution was returned while no integer
886 solution isn't found yet at this time
887 */
888 (MIP_count(lp) == 0) &&
889 ((lp->simplex_mode & (SIMPLEX_Phase2_PRIMAL | SIMPLEX_Phase2_DUAL)) > 0)) {
890 lp->solutioncount++;
891 construct_solution(lp, NULL);
892 transfer_solution(lp, TRUE);
893 }
894 /* Return messages */
895 report(lp, NORMAL, "\nlp_solve optimization was stopped %s.\n",
896 ((status == USERABORT) ? "by the user" : "due to time-out"));
897 }
898 else if(BB->varno == 0)
899 report(lp, NORMAL, "The model %s\n",
900 (status == UNBOUNDED) ? "is UNBOUNDED" :
901 ((status == INFEASIBLE) ? "is INFEASIBLE" : "FAILED"));
902 else {
903 #ifdef Paranoia
904 if((status != FATHOMED) && (status != INFEASIBLE))
905 report(lp, SEVERE, "spx_solve: Invalid return code %d during B&B\n", status);
906 #endif
907 /* If we fathomed a node due to an inferior OF having been detected, return infeasible */
908 if(status == FATHOMED)
909 lp->spx_status = INFEASIBLE;
910 }
911 }
912
913 else { /* ... there is a good solution */
914 construct_solution(lp, NULL);
915 if((lp->bb_level <= 1) && (restored > 0))
916 report(lp, DETAILED, "%s numerics encountered; validate accuracy\n",
917 (restored == 1) ? "Difficult" : "Severe");
918 /* Handle case where a user bound on the OF was found to
919 have been set too aggressively, giving an infeasible model */
920 if(lp->spx_status != OPTIMAL)
921 status = lp->spx_status;
922
923 else if((lp->bb_totalnodes == 0) && (MIP_count(lp) > 0)) {
924 if(lp->lag_status != RUNNING) {
925 report(lp, NORMAL, "\nRelaxed solution " RESULTVALUEMASK " after %10.0f iter is B&B base.\n",
926 lp->solution[0], (double) lp->total_iter);
927 report(lp, NORMAL, " \n");
928 }
929 if((lp->usermessage != NULL) && (lp->msgmask & MSG_LPOPTIMAL)) {
930 REAL *best_solution = lp->best_solution;
931
932 /* transfer_solution(lp, TRUE); */
933 lp->best_solution = lp->solution;
934 lp->usermessage(lp, lp->msghandle, MSG_LPOPTIMAL);
935 lp->best_solution = best_solution;
936 }
937 set_var_priority(lp);
938 }
939
940 /* Check if we have a numeric problem (an earlier version of this code used the
941 absolute difference, but it is not robust for large-valued OFs) */
942 testOF = my_chsign(is_maxim(lp), my_reldiff(lp->solution[0], lp->real_solution));
943 if(testOF < -lp->epsprimal) {
944 report(lp, DETAILED, "solve_LP: A MIP subproblem returned a value better than the base.\n");
945 status = INFEASIBLE;
946 lp->spx_status = status;
947 set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
948 }
949 else if(testOF < 0) /* Avoid problems later (could undo integer roundings, but usually Ok) */
950 lp->solution[0] = lp->real_solution;
951
952 }
953
954 /* status can have the following values:
955 OPTIMAL, SUBOPTIMAL, TIMEOUT, USERABORT, PROCFAIL, UNBOUNDED and INFEASIBLE. */
956
957 return( status );
958 } /* solve_LP */
959
960 STATIC BBrec *findself_BB(BBrec *BB)
961 {
962 int varno = BB->varno, vartype = BB->vartype;
963
964 BB = BB->parent;
965 while((BB != NULL) && (BB->vartype != vartype) && (BB->varno != varno))
966 BB = BB->parent;
967 return( BB );
968 }
969
970 /* Function to determine the opportunity for variable fixing and bound
971 tightening based on a previous best MILP solution and a variable's
972 reduced cost at the current relaxation - inspired by Wolsley */
973 STATIC int rcfbound_BB(BBrec *BB, int varno, MYBOOL isINT, REAL *newbound, MYBOOL *isfeasible)
974 {
975 int i = FR;
976 lprec *lp = BB->lp;
977 REAL deltaRC, rangeLU, deltaOF, lowbo, upbo;
978
979 /* Make sure we only accept non-basic variables */
980 if(lp->is_basic[varno])
981 return( i );
982
983 /* Make sure we only accept non-fixed variables */
984 lowbo = BB->lowbo[varno];
985 upbo = BB->upbo[varno];
986 rangeLU = upbo - lowbo;
987
988 if(rangeLU > lp->epsprimal) {
989 #if 1 /* v5.5 problematic - Gap between current node and the current best bound */
990 deltaOF = lp->rhs[0] - lp->bb_workOF;
991 #elif 0 /* v6 less aggressive - Gap between current best bound and the relaxed problem */
992 deltaOF = my_chsign(is_maxim(lp), lp->real_solution) - lp->bb_workOF;
993 #else /* v6 more aggressive - Gap between current node and the relaxed problem */
994 deltaOF = my_chsign(is_maxim(lp), lp->real_solution) - lp->rhs[0];
995 #endif
996
997 deltaRC = my_chsign(!lp->is_lower[varno], lp->drow[varno]);
998 /* Protect against divisions with tiny numbers and stray sign
999 reversals of the reduced cost */
1000 if(deltaRC < lp->epspivot)
1001 return( i );
1002 deltaRC = deltaOF / deltaRC; /* Should always be a positive number! */
1003 #ifdef Paranoia
1004 if(deltaRC <= 0)
1005 report(lp, SEVERE, "rcfbound_BB: A negative bound fixing level was encountered after node %.0f\n",
1006 (double) lp->bb_totalnodes);
1007 #endif
1008
1009 /* Check if bound implied by the reduced cost is less than existing range */
1010 if(deltaRC < rangeLU + lp->epsint) {
1011 if(lp->is_lower[varno]) {
1012 if(isINT)
1013 deltaRC = scaled_floor(lp, varno, unscaled_value(lp, deltaRC, varno)+lp->epsprimal, 1);
1014 upbo = lowbo + deltaRC;
1015 deltaRC = upbo;
1016 i = LE; /* Sets the upper bound */
1017 }
1018 else {
1019 if(isINT)
1020 deltaRC = scaled_ceil(lp, varno, unscaled_value(lp, deltaRC, varno)+lp->epsprimal, 1);
1021 lowbo = upbo - deltaRC;
1022 deltaRC = lowbo;
1023 i = GE; /* Sets the lower bound */
1024 }
1025
1026 /* Check and set feasibility status */
1027 if((isfeasible != NULL) && (upbo - lowbo < -lp->epsprimal))
1028 *isfeasible = FALSE;
1029
1030 /* Flag that we can fix the variable by returning the relation code negated */
1031 else if(fabs(upbo - lowbo) < lp->epsprimal)
1032 i = -i;
1033 if(newbound != NULL) {
1034 my_roundzero(deltaRC, lp->epsprimal);
1035 *newbound = deltaRC;
1036 }
1037 }
1038
1039 }
1040 return( i );
1041 }
1042
1043
1044 STATIC MYBOOL findnode_BB(BBrec *BB, int *varno, int *vartype, int *varcus)
1045 {
1046 int countsossc, countnint, k, reasonmsg = MSG_NONE;
1047 REAL varsol;
1048 MYBOOL is_better = FALSE, is_equal = FALSE, is_feasible = TRUE;
1049 lprec *lp = BB->lp;
1050
1051 /* Initialize result and return variables */
1052 *varno = 0;
1053 *vartype = BB_REAL;
1054 *varcus = 0;
1055 countnint = 0;
1056 BB->nodestatus = lp->spx_status;
1057 BB->noderesult = lp->solution[0];
1058
1059 /* If this solution is worse than the best so far, this branch dies.
1060 If we can only have integer OF values, and we only need the first solution
1061 then the OF must be at least (unscaled) 1 better than the best so far */
1062 if((lp->bb_limitlevel != 1) && (MIP_count(lp) > 0)) {
1063
1064 /* Check that we don't have a limit on the recursion level; two versions supported:
1065 1) Absolute B&B level (bb_limitlevel > 0), and
1066 2) B&B level relative to the "B&B order" (bb_limitlevel < 0). */
1067 countsossc = lp->sos_vars + lp->sc_vars;
1068 if((lp->bb_limitlevel > 0) && (lp->bb_level > lp->bb_limitlevel+countsossc))
1069 return( FALSE );
1070 else if((lp->bb_limitlevel < 0) &&
1071 (lp->bb_level > 2*(lp->int_vars+countsossc)*abs(lp->bb_limitlevel))) {
1072 if(lp->bb_limitlevel == DEF_BB_LIMITLEVEL)
1073 report(lp, IMPORTANT, "findnode_BB: Default B&B limit reached at %d; optionally change strategy or limit.\n\n",
1074 lp->bb_level);
1075 return( FALSE );
1076 }
1077
1078 /* First initialize or update pseudo-costs from previous optimal solution */
1079 if(BB->varno == 0) {
1080 varsol = lp->infinite;
1081 if((lp->int_vars+lp->sc_vars > 0) && (lp->bb_PseudoCost == NULL))
1082 lp->bb_PseudoCost = init_pseudocost(lp, get_bb_rule(lp));
1083 }
1084 else {
1085 varsol = lp->solution[BB->varno];
1086 if( ((lp->int_vars > 0) && (BB->vartype == BB_INT)) ||
1087 ((lp->sc_vars > 0) && (BB->vartype == BB_SC) && !is_int(lp, BB->varno-lp->rows)) )
1088 update_pseudocost(lp->bb_PseudoCost, BB->varno-lp->rows, BB->vartype, BB->isfloor, varsol);
1089 }
1090
1091 /* Make sure we don't have numeric problems (typically due to integer scaling) */
1092 if((lp->bb_totalnodes > 0) && !bb_better(lp, OF_RELAXED, OF_TEST_WE)) {
1093 if(lp->bb_trace)
1094 report(lp, IMPORTANT, "findnode_BB: Simplex failure due to loss of numeric accuracy\n");
1095 lp->spx_status = NUMFAILURE;
1096 return( FALSE );
1097 }
1098
1099 /* Abandon this branch if the solution is "worse" than a heuristically
1100 determined limit or the previous best MIP solution */
1101 if(((lp->solutioncount == 0) && !bb_better(lp, OF_HEURISTIC, OF_TEST_BE)) ||
1102 ((lp->solutioncount > 0) &&
1103 (!bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BE | OF_TEST_RELGAP) ||
1104 !bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BE)))) {
1105 return( FALSE );
1106 }
1107
1108 /* Collect violated SC variables (since they can also be real-valued); the
1109 approach is to get them out of the way, since a 0-value is assumed to be "cheap" */
1110 if(lp->sc_vars > 0) {
1111 *varno = find_sc_bbvar(lp, &countnint);
1112 if(*varno > 0)
1113 *vartype = BB_SC;
1114 }
1115
1116 /* Look among SOS variables if no SC candidate was found */
1117 if((SOS_count(lp) > 0) && (*varno == 0)) {
1118 *varno = find_sos_bbvar(lp, &countnint, FALSE);
1119 if(*varno < 0)
1120 *varno = 0;
1121 else if(*varno > 0)
1122 *vartype = BB_SOS;
1123 }
1124
1125 /* Then collect INTS that are not integer valued, and verify bounds */
1126 if((lp->int_vars > 0) && (*varno == 0)) {
1127 *varno = find_int_bbvar(lp, &countnint, BB, &is_feasible);
1128 if(*varno > 0) {
1129 *vartype = BB_INT;
1130 if((countnint == 1) && !is_feasible) {
1131 BB->lastrcf = 0;
1132 return( FALSE );
1133 }
1134 }
1135 }
1136
1137 #if 1 /* peno */
1138 /* Check if we have reached the depth limit for any individual variable
1139 (protects against infinite recursions of mainly integer variables) */
1140 k = *varno-lp->rows;
1141 if((*varno > 0) && (lp->bb_limitlevel != 0) && (lp->bb_varactive[k] >= abs(lp->bb_limitlevel) /* abs(DEF_BB_LIMITLEVEL) */)) {
1142 /* if(!is_action(lp->nomessage, NOMSG_BBLIMIT)) {*/
1143 /*
1144 report(lp, IMPORTANT, "findnode_BB: Reached B&B depth limit %d for variable %d; will not dive further.\n\n",
1145 lp->bb_varactive[k], k);
1146 */
1147 /* set_action(&lp->nomessage, NOMSG_BBLIMIT); */
1148 /* } */
1149 return( FALSE );
1150 }
1151 #endif
1152
1153 /* Check if the current MIP solution is optimal; equal or better */
1154 if(*varno == 0) {
1155 is_better = (MYBOOL) (lp->solutioncount == 0) || bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT);
1156 #if 1
1157 is_better &= bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT | OF_TEST_RELGAP);
1158 #else
1159 /* Check if we can determine clear improvement */
1160 is_better = (MYBOOL) (lp->solutioncount == 0) ||
1161 (MYBOOL) ((lp->bb_deltaOF > 0) &&
1162 (my_chsign(is_maxim(lp), lp->solution[0]-lp->best_solution[0]) < 0));
1163
1164 /* Apply gap-based improvement testing if the current solution is not clearly better */
1165
1166 if(!is_better) {
1167 is_better = bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT);
1168 is_better |= bb_better(lp, OF_INCUMBENT | OF_DELTA, OF_TEST_BT | OF_TEST_RELGAP);
1169 }
1170 #endif
1171 is_equal = !is_better;
1172
1173 if(is_equal) {
1174 if((lp->solutionlimit <= 0) || (lp->solutioncount < lp->solutionlimit)) {
1175 lp->solutioncount++;
1176 SETMIN(lp->bb_solutionlevel, lp->bb_level);
1177 reasonmsg = MSG_MILPEQUAL;
1178 }
1179 }
1180
1181 /* Current solution is better */
1182 else if(is_better) {
1183
1184 /* Update grand total solution count and check if we should go from
1185 depth-first to best-first variable selection mode */
1186 if(lp->bb_varactive != NULL) {
1187 lp->bb_varactive[0]++;
1188 if((lp->bb_varactive[0] == 1) &&
1189 is_bb_mode(lp, NODE_DEPTHFIRSTMODE) && is_bb_mode(lp, NODE_DYNAMICMODE))
1190 lp->bb_rule &= !NODE_DEPTHFIRSTMODE;
1191 }
1192
1193 if(lp->bb_trace ||
1194 ((lp->verbose >= NORMAL) && (lp->print_sol == FALSE) && (lp->lag_status != RUNNING))) {
1195 report(lp, IMPORTANT,
1196 "%s solution " RESULTVALUEMASK " after %10.0f iter, %9.0f nodes (gap %.1f%%)\n",
1197 (lp->bb_improvements == 0) ? "Feasible" : "Improved",
1198 lp->solution[0], (double) lp->total_iter, (double) lp->bb_totalnodes,
1199 100.0*fabs(my_reldiff(lp->solution[0], lp->bb_limitOF)));
1200 }
1201 if(MIP_count(lp) > 0) {
1202 if(lp->bb_improvements == 0)
1203 reasonmsg = MSG_MILPFEASIBLE;
1204 else
1205 reasonmsg = MSG_MILPBETTER;
1206 }
1207
1208 lp->bb_status = FEASFOUND;
1209 lp->bb_solutionlevel = lp->bb_level;
1210 lp->solutioncount = 1;
1211 lp->bb_improvements++;
1212 lp->bb_workOF = lp->rhs[0];
1213
1214 if(lp->bb_breakfirst ||
1215 (!is_infinite(lp, lp->bb_breakOF) && bb_better(lp, OF_USERBREAK, OF_TEST_BE)))
1216 lp->bb_break = TRUE;
1217 }
1218 }
1219 }
1220 else {
1221 is_better = TRUE;
1222 lp->solutioncount = 1;
1223 }
1224
1225 /* Transfer the successful solution vector */
1226 if(is_better || is_equal) {
1227 #ifdef ParanoiaMIP
1228 if((lp->bb_level > 0) &&
1229 (check_solution(lp, lp->columns, lp->solution,
1230 lp->orig_upbo, lp->orig_lowbo, lp->epssolution) != OPTIMAL)) {
1231 lp->solutioncount = 0;
1232 lp->spx_status = NUMFAILURE;
1233 lp->bb_status = lp->spx_status;
1234 lp->bb_break = TRUE;
1235 return( FALSE );
1236 }
1237 #endif
1238 transfer_solution(lp, (MYBOOL) ((lp->do_presolve & PRESOLVE_LASTMASKMODE) != PRESOLVE_NONE));
1239 if((MIP_count(lp) > 0) && (lp->bb_totalnodes > 0)) {
1240 if ((!construct_duals(lp)) ||
1241 (is_presolve(lp, PRESOLVE_SENSDUALS) &&
1242 (!construct_sensitivity_duals(lp) || !construct_sensitivity_obj(lp))
1243 )
1244 ) {
1245 }
1246 }
1247 if((reasonmsg != MSG_NONE) && (lp->msgmask & reasonmsg) && (lp->usermessage != NULL))
1248 lp->usermessage(lp, lp->msghandle, reasonmsg);
1249
1250 if(lp->print_sol != FALSE) {
1251 print_objective(lp);
1252 print_solution(lp, 1);
1253 }
1254 }
1255
1256 /* Do tracing and determine if we have arrived at the estimated lower MIP limit */
1257 *varcus = countnint;
1258 if(MIP_count(lp) > 0) {
1259 if((countnint == 0) && (lp->solutioncount == 1) && (lp->solutionlimit == 1) &&
1260 (bb_better(lp, OF_DUALLIMIT, OF_TEST_BE) || bb_better(lp, OF_USERBREAK, OF_TEST_BE | OF_TEST_RELGAP))) {
1261 lp->bb_break = (MYBOOL) (countnint == 0);
1262 return( FALSE );
1263 }
1264 else if(lp->bb_level > 0) {
1265 #ifdef MIPboundWithOF
1266 if((lp->constraintOF > 0) && (countnint == 0))
1267 set_rh(lp, lp->constraintOF, lp->solution[0] + my_chsign(!is_maxim(lp), lp->bb_deltaOF));
1268 #endif
1269 if(lp->spx_trace)
1270 report(lp, DETAILED, "B&B level %5d OPT %16s value " RESULTVALUEMASK "\n",
1271 lp->bb_level, (*varno) ? " " : "INT", lp->solution[0]);
1272 }
1273 return( (MYBOOL) (*varno > 0));
1274 }
1275 else
1276 return( FALSE );
1277
1278 }
1279
1280 STATIC int solve_BB(BBrec *BB)
1281 {
1282 int K, status;
1283 lprec *lp = BB->lp;
1284
1285 /* Protect against infinite recursions do to integer rounding effects */
1286 status = PROCFAIL;
1287
1288 /* Shortcut variables, set default bounds */
1289 K = BB->varno;
1290
1291 /* Load simple MIP bounds */
1292 if(K > 0) {
1293
1294 /* Update cuts, if specified */
1295 updatecuts_BB(lp);
1296
1297 /* BRANCH_FLOOR: Force the variable to be smaller than the B&B upper bound */
1298 if(BB->isfloor)
1299 modifyUndoLadder(lp->bb_upperchange, K, BB->upbo, BB->UPbound);
1300
1301 /* BRANCH_CEILING: Force the variable to be greater than the B&B lower bound */
1302 else
1303 modifyUndoLadder(lp->bb_lowerchange, K, BB->lowbo, BB->LObound);
1304
1305 /* Update MIP node count */
1306 BB->nodessolved++;
1307
1308 }
1309
1310 /* Solve! */
1311 status = solve_LP(lp, BB);
1312
1313 /* Do special feasibility assessment of high order SOS'es */
1314 #if 1
1315 if((status == OPTIMAL) && (BB->vartype == BB_SOS) && !SOS_is_feasible(lp->SOS, 0, lp->solution))
1316 status = INFEASIBLE;
1317 #endif
1318
1319 return( status );
1320 }
1321
1322 /* Routine to compute a "strong" pseudo-cost update for a node */
1323 STATIC MYBOOL strongbranch_BB(lprec *lp, BBrec *BB, int varno, int vartype, int varcus)
1324 {
1325 MYBOOL success = FALSE;
1326 int i;
1327 BBrec *strongBB;
1328
1329 /* Create new B&B level and solve each of the branches */
1330 lp->is_strongbranch = TRUE;
1331 push_basis(lp, lp->var_basic, lp->is_basic, lp->is_lower);
1332 strongBB = push_BB(lp, BB, lp->rows+varno, vartype, varcus);
1333 if(strongBB == BB)
1334 return( success );
1335
1336 do {
1337
1338 /* Solve incremental problem to local optimality */
1339 lp->bb_strongbranches++;
1340 /* set_action(&lp->spx_action, ACTION_REBASE | ACTION_RECOMPUTE); */
1341 if(solve_BB(strongBB) == OPTIMAL) {
1342
1343 /* Update result indicator*/
1344 success |= 1 << strongBB->isfloor;
1345
1346 /* Compute new count of non-ints */
1347 strongBB->lastvarcus = 0;
1348 for(i = 1; i <= lp->columns; i++) {
1349 if(is_int(lp, i) && !solution_is_int(lp, lp->rows+i, FALSE))
1350 strongBB->lastvarcus++;
1351 }
1352
1353 /* Perform the pseudo-cost update */
1354 update_pseudocost(lp->bb_PseudoCost, varno, strongBB->vartype, strongBB->isfloor,
1355 lp->solution[strongBB->varno]);
1356 }
1357 }
1358 while(nextbranch_BB(strongBB));
1359
1360 strongBB = pop_BB(strongBB);
1361 if(strongBB != BB)
1362 report(lp, SEVERE, "strongbranch_BB: Invalid bound settings restored for variable %d\n",
1363 varno);
1364 pop_basis(lp, TRUE);
1365 set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
1366
1367 lp->is_strongbranch = FALSE;
1368
1369 return( success );
1370 }
1371
1372 /* Future functions */
1373 STATIC MYBOOL pre_BB(lprec *lp)
1374 {
1375 return( TRUE );
1376 }
1377 STATIC MYBOOL post_BB(lprec *lp)
1378 {
1379 return( TRUE );
1380 }
1381
1382 /* This is the non-recursive B&B driver routine - beautifully simple, yet so subtle! */
1383 STATIC int run_BB(lprec *lp)
1384 {
1385 BBrec *currentBB;
1386 int varno, vartype, varcus, prevsolutions;
1387 int status = NOTRUN;
1388
1389 /* Initialize */
1390 pre_BB(lp);
1391 prevsolutions = lp->solutioncount;
1392 #ifdef UseMilpSlacksRCF /* Check if we should include ranged constraints */
1393 varno = lp->sum;
1394 #else
1395 varno = lp->columns;
1396 #endif
1397 lp->bb_upperchange = createUndoLadder(lp, varno, 2*MIP_count(lp));
1398 lp->bb_lowerchange = createUndoLadder(lp, varno, 2*MIP_count(lp));
1399 lp->rootbounds = currentBB = push_BB(lp, NULL, 0, BB_REAL, 0);
1400
1401 /* Perform the branch & bound loop */
1402 while(lp->bb_level > 0) {
1403 status = solve_BB(currentBB);
1404
1405 #if 0
1406 if((lp->bb_level == 1) && (MIP_count(lp) > 0)) {
1407 if(status == RUNNING)
1408 ;
1409
1410 /* Check if there was an integer solution of an aborted model */
1411 else if((status == SUBOPTIMAL) && (lp->solutioncount == 1) &&
1412 findnode_BB(currentBB, &varno, &vartype, &varcus))
1413 status = USERABORT;
1414 }
1415 #endif
1416
1417 if((status == OPTIMAL) && findnode_BB(currentBB, &varno, &vartype, &varcus))
1418 currentBB = push_BB(lp, currentBB, varno, vartype, varcus);
1419
1420 else while((lp->bb_level > 0) && !nextbranch_BB(currentBB))
1421 currentBB = pop_BB(currentBB);
1422
1423 }
1424
1425 /* Finalize */
1426 freeUndoLadder(&(lp->bb_upperchange));
1427 freeUndoLadder(&(lp->bb_lowerchange));
1428
1429 /* Check if we should adjust status */
1430 if(lp->solutioncount > prevsolutions) {
1431 if((status == PROCBREAK) || (status == USERABORT) || (status == TIMEOUT) || userabort(lp, -1))
1432 status = SUBOPTIMAL;
1433 else
1434 status = OPTIMAL;
1435 if(lp->bb_totalnodes > 0)
1436 lp->spx_status = OPTIMAL;
1437 }
1438 post_BB(lp);
1439 return( status );
1440 }
1441
1442