3
/*----------------------------------------------------------------------
4
-- This code is part of GNU Linear Programming Kit (GLPK).
6
-- Copyright (C) 2000, 01, 02, 03, 04, 05, 06 Andrew Makhorin,
7
-- Department for Applied Informatics, Moscow Aviation Institute,
8
-- Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
10
-- GLPK is free software; you can redistribute it and/or modify it
11
-- under the terms of the GNU General Public License as published by
12
-- the Free Software Foundation; either version 2, or (at your option)
15
-- GLPK is distributed in the hope that it will be useful, but WITHOUT
16
-- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17
-- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
18
-- License for more details.
20
-- You should have received a copy of the GNU General Public License
21
-- along with GLPK; see the file COPYING. If not, write to the Free
22
-- Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
24
----------------------------------------------------------------------*/
32
/*----------------------------------------------------------------------
33
-- show_progress - display progress of the search.
35
-- This routine displays some information about progress of the search.
37
-- This information includes:
39
-- the current number of iterations performed by the simplex solver;
41
-- the objective value for the best known integer feasible solution,
42
-- which is upper (minimization) or lower (maximization) global bound
43
-- for optimal solution of the original mip problem;
45
-- the best local bound for active nodes, which is lower (minimization)
46
-- or upper (maximization) global bound for optimal solution of the
47
-- original mip problem;
49
-- the relative mip gap, in percents;
51
-- the number of open (active) subproblems;
53
-- the number of completely explored subproblems, i.e. whose nodes have
54
-- been already removed from the tree. */
56
static void show_progress(MIPTREE *tree)
59
char best_mip[50], best_bound[50], *rho, rel_gap[50];
60
/* format the best known integer feasible solution */
62
sprintf(best_mip, "%17.9e", tree->best);
64
sprintf(best_mip, "%17s", "not found yet");
65
/* determine reference number of an active subproblem whose local
67
p = mip_best_node(tree);
68
/* format the best bound */
70
sprintf(best_bound, "%17s", "tree is empty");
72
{ temp = tree->slot[p].node->bound;
74
sprintf(best_bound, "%17s", "-inf");
75
else if (temp == +DBL_MAX)
76
sprintf(best_bound, "%17s", "+inf");
78
sprintf(best_bound, "%17.9e", temp);
80
/* choose the relation sign between global bounds */
82
{ case LPX_MIN: rho = ">="; break;
83
case LPX_MAX: rho = "<="; break;
84
default: insist(tree != tree);
86
/* format the relative mip gap */
87
temp = mip_relative_gap(tree);
89
sprintf(rel_gap, " 0.0%%");
90
else if (temp < 0.001)
91
sprintf(rel_gap, "< 0.1%%");
92
else if (temp <= 9.999)
93
sprintf(rel_gap, "%5.1f%%", 100.0 * temp);
95
sprintf(rel_gap, "%6s", "");
96
/* display progress of the search */
97
print("+%6d: mip = %s %s %s %s (%d; %d)",
98
lpx_get_int_parm(tree->lp, LPX_K_ITCNT), best_mip, rho,
99
best_bound, rel_gap, tree->a_cnt, tree->t_cnt - tree->n_cnt);
100
tree->tm_lag = utime();
104
/*----------------------------------------------------------------------
105
-- set_local_bound - set local bound for current subproblem.
107
-- This routine sets the local bound for integer optimal solution of
108
-- the current subproblem, which is optimal solution of LP relaxation.
109
-- If the latter is fractional while the objective function is known to
110
-- be integral for any integer feasible point, the bound is strengthen
111
-- by rounding up (minimization) or down (maximization). */
113
static void set_local_bound(MIPTREE *tree)
115
bound = lpx_get_obj_val(tree->lp);
117
{ /* the objective function is known to be integral */
119
temp = gnm_floor(bound + 0.5);
120
if (temp - 1e-5 <= bound && bound <= temp + 1e-5)
125
bound = gnm_ceil(bound); break;
127
bound = gnm_floor(bound); break;
129
insist(tree != tree);
133
insist(tree->curr != NULL);
134
tree->curr->bound = bound;
135
if (tree->msg_lev >= 3)
136
print("Local bound is %.9e", bound);
140
/*----------------------------------------------------------------------
141
-- is_branch_hopeful - check if specified branch is hopeful.
143
-- This routine checks if the specified subproblem can have an integer
144
-- optimal solution which is better than the best known one.
146
-- The check is based on comparison of the local objective bound stored
147
-- in the subproblem descriptor and the incumbent objective value which
148
-- is the global objective bound.
150
-- If there is a chance that the specified subproblem can have a better
151
-- integer optimal solution, the routine returns non-zero. Otherwise, if
152
-- the corresponding branch can pruned, zero is returned. */
154
static int is_branch_hopeful(MIPTREE *tree, int p)
157
{ gnm_float bound, eps;
158
insist(1 <= p && p <= tree->nslots);
159
insist(tree->slot[p].node != NULL);
160
bound = tree->slot[p].node->bound;
161
eps = tree->tol_obj * (1.0 + gnm_abs(tree->best));
164
if (bound >= tree->best - eps) ret = 0;
167
if (bound <= tree->best + eps) ret = 0;
170
insist(tree != tree);
176
/*----------------------------------------------------------------------
177
-- check_integrality - check integrality of basic solution.
179
-- This routine checks if the basic solution of LP relaxation of the
180
-- current subproblem satisfies to integrality conditions, i.e. that all
181
-- variables of integer kind have integral primal values. (The solution
182
-- is assumed to be optimal.)
184
-- For each variable of integer kind the routine computes the following
187
-- ii(x[j]) = min(x[j] - gnm_floor(x[j]), gnm_ceil(x[j]) - x[j]), (1)
189
-- which is a measure of the integer infeasibility (non-integrality) of
190
-- x[j] (for example, ii(2.1) = 0.1, ii(3.7) = 0.3, ii(5.0) = 0). It is
191
-- understood that 0 <= ii(x[j]) <= 0.5, and variable x[j] is integer
192
-- feasible if ii(x[j]) = 0. However, due to floating-point arithmetic
193
-- the routine checks less restrictive condition:
195
-- ii(x[j]) <= tol_int, (2)
197
-- where tol_int is a given tolerance (small positive number) and marks
198
-- each variable which does not satisfy to (2) as integer infeasible by
199
-- setting its fractionality flag.
201
-- In order to characterize integer infeasibility of the basic solution
202
-- in the whole the routine computes two parameters: ii_cnt, which is
203
-- the number of variables with the fractionality flag set, and ii_sum,
204
-- which is the sum of integer infeasibilities (1). */
206
static void check_integrality(MIPTREE *tree)
207
{ LPX *lp = tree->lp;
208
int j, type, ii_cnt = 0;
209
gnm_float lb, ub, x, temp1, temp2, ii_sum = 0.0;
210
/* walk through the set of columns (structural variables) */
211
for (j = 1; j <= tree->n; j++)
212
{ tree->non_int[j] = 0;
213
/* if the column is continuous, skip it */
214
if (!tree->int_col[j]) continue;
215
/* if the column is non-basic, it is integer feasible */
216
if (lpx_get_col_stat(lp, j) != LPX_BS) continue;
217
/* obtain the type and bounds of the column */
218
type = lpx_get_col_type(lp, j);
219
lb = lpx_get_col_lb(lp, j);
220
ub = lpx_get_col_ub(lp, j);
221
/* obtain value of the column in optimal basic solution */
222
x = lpx_get_col_prim(lp, j);
223
/* if the column's primal value is close to the lower bound,
224
the column is integer feasible within given tolerance */
225
if (type == LPX_LO || type == LPX_DB || type == LPX_FX)
226
{ temp1 = lb - tree->tol_int;
227
temp2 = lb + tree->tol_int;
228
if (temp1 <= x && x <= temp2) continue;
229
/* the lower bound must not be violated */
232
/* if the column's primal value is close to the upper bound,
233
the column is integer feasible within given tolerance */
234
if (type == LPX_UP || type == LPX_DB || type == LPX_FX)
235
{ temp1 = ub - tree->tol_int;
236
temp2 = ub + tree->tol_int;
237
if (temp1 <= x && x <= temp2) continue;
238
/* the upper bound must not be violated */
241
/* if the column's primal value is close to nearest integer,
242
the column is integer feasible within given tolerance */
243
temp1 = gnm_floor(x + 0.5) - tree->tol_int;
244
temp2 = gnm_floor(x + 0.5) + tree->tol_int;
245
if (temp1 <= x && x <= temp2) continue;
246
/* otherwise the column is integer infeasible */
247
tree->non_int[j] = 1;
248
/* increase the number of fractional-valued columns */
250
/* compute the sum of integer infeasibilities */
251
temp1 = x - gnm_floor(x);
252
temp2 = gnm_ceil(x) - x;
253
insist(temp1 > 0.0 && temp2 > 0.0);
254
ii_sum += (temp1 <= temp2 ? temp1 : temp2);
256
/* store ii_cnt and ii_sum in the current problem descriptor */
257
insist(tree->curr != NULL);
258
tree->curr->ii_cnt = ii_cnt;
259
tree->curr->ii_sum = ii_sum;
260
/* and also display these parameters */
261
if (tree->msg_lev >= 3)
263
print("There are no fractional columns");
264
else if (ii_cnt == 1)
265
print("There is one fractional column, integer infeasibilit"
266
"y is %.3e", ii_sum);
268
print("There are %d fractional columns, integer infeasibili"
269
"ty is %.3e", ii_cnt, ii_sum);
274
/*----------------------------------------------------------------------
275
-- record_solution - record better integer feasible solution.
277
-- This routine records optimal basic solution of LP relaxation of the
278
-- current subproblem, which being integer feasible is better than the
279
-- best known integer feasible solution. */
281
static void record_solution(MIPTREE *tree)
288
tree->best = lpx_get_obj_val(lp);
289
for (i = 1; i <= m; i++)
290
{ temp = lpx_get_row_prim(lp, i);
291
tree->mipx[i] = temp;
293
for (j = 1; j <= n; j++)
294
{ temp = lpx_get_col_prim(lp, j);
295
/* value of the integer column must be integral */
296
if (tree->int_col[j]) temp = gnm_floor(temp + 0.5);
297
tree->mipx[m+j] = temp;
302
/*----------------------------------------------------------------------
303
-- fix_by_red_cost - fix non-basic integer columns by reduced costs.
305
-- This routine fixes some non-basic integer columns if their reduced
306
-- costs indicate that increasing (decreasing) the column at least by
307
-- one involves the objective value becoming worse than the incumbent
308
-- objective value (i.e. the global bound). */
310
static void fix_by_red_cost(MIPTREE *tree)
313
int j, stat, fixed = 0;
314
gnm_float obj, lb, ub, dj;
315
/* the global bound must exist */
317
/* basic solution of LP relaxation must be optimal */
318
insist(lpx_get_status(lp) == LPX_OPT);
319
/* determine the objective function value */
320
obj = lpx_get_obj_val(lp);
321
/* walk through the column list */
322
for (j = 1; j <= n; j++)
323
{ /* if j-th column is not of integer kind, skip it */
324
if (!tree->int_col[j]) continue;
325
/* obtain bounds of j-th column */
326
lb = lpx_get_col_lb(lp, j);
327
ub = lpx_get_col_ub(lp, j);
328
/* and determine its status and reduced cost */
329
stat = lpx_get_col_stat(lp, j);
330
dj = lpx_get_col_dual(lp, j);
331
/* analyze the reduced cost */
336
{ /* j-th column is non-basic on its lower bound */
337
if (dj < 0.0) dj = 0.0;
338
if (obj + dj >= tree->best)
339
lpx_set_col_bnds(lp, j, LPX_FX, lb, lb), fixed++;
341
else if (stat == LPX_NU)
342
{ /* j-th column is non-basic on its upper bound */
343
if (dj > 0.0) dj = 0.0;
344
if (obj - dj >= tree->best)
345
lpx_set_col_bnds(lp, j, LPX_FX, ub, ub), fixed++;
351
{ /* j-th column is non-basic on its lower bound */
352
if (dj > 0.0) dj = 0.0;
353
if (obj + dj <= tree->best)
354
lpx_set_col_bnds(lp, j, LPX_FX, lb, lb), fixed++;
356
else if (stat == LPX_NU)
357
{ /* j-th column is non-basic on its upper bound */
358
if (dj < 0.0) dj = 0.0;
359
if (obj - dj <= tree->best)
360
lpx_set_col_bnds(lp, j, LPX_FX, ub, ub), fixed++;
364
insist(tree != tree);
367
if (tree->msg_lev >= 3)
369
/* nothing to say */;
371
print("One column has been fixed by reduced cost");
373
print("%d columns have been fixed by reduced costs", fixed);
375
/* fixing non-basic columns on their current bounds does not
376
change the basic solution itself, however, lpx_set_col_bnds
377
resets the primal and dual statuses to undefined state, and
378
therefore we need to restore them */
379
lpx_put_solution(lp, LPX_P_FEAS, LPX_D_FEAS, NULL, NULL, NULL,
384
/*----------------------------------------------------------------------
385
-- branch_on - perform branching on specified column.
387
-- This routine performs branching on j-th column (structural variable)
388
-- of the current subproblem. The specified column must be of integer
389
-- kind and must have a fractional value in optimal basic solution of
390
-- LP relaxation of the current subproblem (i.e. only columns for which
391
-- the flag non_int[j] is set are valid candidates to branch on).
393
-- Let x be j-th structural variable, and beta be its primal fractional
394
-- value in the current basic solution. Branching on j-th variable is
395
-- dividing the current subproblem into two new subproblems, which are
396
-- identical to the current subproblem with the following exception: in
397
-- the first subproblem that begins the down-branch x has a new upper
398
-- bound x <= gnm_floor(beta), and in the second subproblem that begins the
399
-- up-branch x has a new lower bound x >= gnm_ceil(beta).
401
-- The parameter next specifies which subproblem should be solved next
402
-- to continue the search:
404
-- -1 means that one which begins the down-branch;
405
-- +1 means that one which begins the up-branch. */
407
static void branch_on(MIPTREE *tree, int j, int next)
408
{ LPX *lp = tree->lp;
409
int p, type, clone[1+2];
410
gnm_float beta, lb, ub, new_lb, new_ub;
411
/* check input parameters for correctness */
412
insist(1 <= j && j <= tree->n);
413
insist(tree->non_int[j]);
414
insist(next == -1 || next == +1);
415
/* obtain primal value of j-th column in basic solution */
416
beta = lpx_get_col_prim(lp, j);
417
if (tree->msg_lev >= 3)
418
print("Branching on column %d, primal value is %.9e", j, beta);
419
/* determine the reference number of the current subproblem */
420
insist(tree->curr != NULL);
422
/* freeze the current subproblem */
423
mip_freeze_node(tree);
424
/* create two clones of the current subproblem; the first clone
425
begins the down-branch, the second one begins the up-branch */
426
mip_clone_node(tree, p, 2, clone);
427
if (tree->msg_lev >= 3)
428
print("Node %d begins down branch, node %d begins up branch",
430
/* set new upper bound of j-th column in the first subproblem */
431
mip_revive_node(tree, clone[1]);
432
type = lpx_get_col_type(lp, j);
433
lb = lpx_get_col_lb(lp, j);
434
ub = lpx_get_col_ub(lp, j);
435
new_ub = gnm_floor(beta);
441
insist(lb <= new_ub);
442
type = (lb < new_ub ? LPX_DB : LPX_FX);
445
insist(new_ub <= ub - 1.0);
448
insist(lb <= new_ub && new_ub <= ub - 1.0);
449
type = (lb < new_ub ? LPX_DB : LPX_FX);
452
insist(type != type);
454
lpx_set_col_bnds(lp, j, type, lb, new_ub);
455
mip_freeze_node(tree);
456
/* set new lower bound of j-th column in the second subproblem */
457
mip_revive_node(tree, clone[2]);
458
type = lpx_get_col_type(lp, j);
459
lb = lpx_get_col_lb(lp, j);
460
ub = lpx_get_col_ub(lp, j);
461
new_lb = gnm_ceil(beta);
467
insist(lb + 1.0 <= new_lb);
470
insist(new_lb <= ub);
471
type = (new_lb < ub ? LPX_DB : LPX_FX);
474
insist(lb + 1.0 <= new_lb && new_lb <= ub);
475
type = (new_lb < ub ? LPX_DB : LPX_FX);
478
insist(type != type);
480
lpx_set_col_bnds(lp, j, type, new_lb, ub);
481
mip_freeze_node(tree);
482
/* revive the subproblem to be solved next */
483
mip_revive_node(tree, clone[next < 0 ? 1 : 2]);
487
/*----------------------------------------------------------------------
488
-- branch_first - choose first branching variable.
490
-- This routine looks up the list of structural variables and chooses
491
-- the first one, which is of integer kind and has fractional value in
492
-- optimal solution of the current LP relaxation.
494
-- This routine also selects the branch to be solved next where integer
495
-- infeasibility of the chosen variable is less than in other one. */
497
static void branch_first(MIPTREE *tree)
502
/* choose the column to branch on */
503
for (j = 1; j <= n; j++)
504
if (tree->non_int[j]) break;
505
insist(1 <= j && j <= n);
506
/* select the branch to be solved next */
507
beta = lpx_get_col_prim(lp, j);
508
if (beta - gnm_floor(beta) < gnm_ceil(beta) - beta)
509
next = -1; /* down branch */
511
next = +1; /* up branch */
512
/* perform branching */
513
branch_on(tree, j, next);
517
/*----------------------------------------------------------------------
518
-- branch_last - choose last branching variable.
520
-- This routine looks up the list of structural variables and chooses
521
-- the last one, which is of integer kind and has fractional value in
522
-- optimal solution of the current LP relaxation.
524
-- This routine also selects the branch to be solved next where integer
525
-- infeasibility of the chosen variable is less than in other one. */
527
static void branch_last(MIPTREE *tree)
532
/* choose the column to branch on */
533
for (j = n; j >= 1; j--)
534
if (tree->non_int[j]) break;
535
insist(1 <= j && j <= n);
536
/* select the branch to be solved next */
537
beta = lpx_get_col_prim(lp, j);
538
if (beta - gnm_floor(beta) < gnm_ceil(beta) - beta)
539
next = -1; /* down branch */
541
next = +1; /* up branch */
542
/* perform branching */
543
branch_on(tree, j, next);
547
/*----------------------------------------------------------------------
548
-- branch_drtom - choose branching variable with Driebeck-Tomlin heur.
550
-- This routine chooses a structural variable, which is required to be
551
-- integral and has fractional value in optimal solution of the current
552
-- LP relaxation, using a heuristic proposed by Driebeck and Tomlin.
554
-- The routine also selects the branch to be solved next, again due to
555
-- Driebeck and Tomlin.
557
-- This routine is based on the heuristic proposed in:
559
-- Driebeck N.J. An algorithm for the solution of mixed-integer
560
-- programming problems, Management Science, 12: 576-87 (1966)
564
-- Tomlin J.A. Branch and bound methods for integer and non-convex
565
-- programming, in J.Abadie (ed.), Integer and Nonlinear Programming,
566
-- North-Holland, Amsterdam, pp. 437-50 (1970).
568
-- Must note that this heuristic is time-expensive, because computing
569
-- one-step degradation (see the routine below) requires one BTRAN for
570
-- every fractional-valued structural variable. */
572
static void branch_drtom(MIPTREE *tree)
576
int j, jj, k, t, next, kase, len, stat, *ind;
577
gnm_float x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,
578
dd_dn, dd_up, degrad, *val;
579
/* basic solution of LP relaxation must be optimal */
580
insist(lpx_get_status(lp) == LPX_OPT);
581
/* allocate working arrays */
582
ind = ucalloc(1+n, sizeof(int));
583
val = ucalloc(1+n, sizeof(gnm_float));
584
/* nothing has been chosen so far */
585
jj = 0, degrad = -1.0;
586
/* walk through the list of columns (structural variables) */
587
for (j = 1; j <= n; j++)
588
{ /* if j-th column is not marked as fractional, skip it */
589
if (!tree->non_int[j]) continue;
590
/* obtain (fractional) value of j-th column in basic solution
592
x = lpx_get_col_prim(lp, j);
593
/* since the value of j-th column is fractional, the column is
594
basic; compute corresponding row of the simplex table */
595
len = lpx_eval_tab_row(lp, m+j, ind, val);
596
/* the following fragment computes a change in the objective
597
function: delta Z = new Z - old Z, where old Z is the
598
objective value in the current optimal basis, and new Z is
599
the objective value in the adjacent basis, for two cases:
600
1) if new upper bound ub' = gnm_floor(x[j]) is introduced for
601
j-th column (down branch);
602
2) if new lower bound lb' = gnm_ceil(x[j]) is introduced for
603
j-th column (up branch);
604
since in both cases the solution remaining dual feasible
605
becomes primal infeasible, one implicit simplex iteration
606
is performed to determine the change delta Z;
607
it is obvious that new Z, which is never better than old Z,
608
is a lower (minimization) or upper (maximization) bound of
609
the objective function for down- and up-branches. */
610
for (kase = -1; kase <= +1; kase += 2)
611
{ /* if kase < 0, the new upper bound of x[j] is introduced;
612
in this case x[j] should decrease in order to leave the
613
basis and go to its new upper bound */
614
/* if kase > 0, the new lower bound of x[j] is introduced;
615
in this case x[j] should increase in order to leave the
616
basis and go to its new lower bound */
617
/* apply the dual ratio test in order to determine which
618
auxiliary or structural variable should enter the basis
619
to keep dual feasibility */
620
k = lpx_dual_ratio_test(lp, len, ind, val, kase, 1e-8);
621
/* if no non-basic variable has been chosen, LP relaxation
622
of corresponding branch being primal infeasible and dual
623
unbounded has no primal feasible solution; in this case
624
the change delta Z is formally set to infinity */
626
{ delta_z = (tree->dir == LPX_MIN ? +DBL_MAX : -DBL_MAX);
629
/* row of the simplex table that corresponds to non-basic
630
variable x[k] choosen by the dual ratio test is:
631
x[j] = ... + alfa * x[k] + ...
632
where alfa is the influence coefficient (an element of
633
the simplex table row) */
634
/* determine the coefficient alfa */
635
for (t = 1; t <= len; t++) if (ind[t] == k) break;
636
insist(1 <= t && t <= len);
638
/* since in the adjacent basis the variable x[j] becomes
639
non-basic, knowing its value in the current basis we can
640
determine its change delta x[j] = new x[j] - old x[j] */
641
delta_j = (kase < 0 ? gnm_floor(x) : gnm_ceil(x)) - x;
642
/* and knowing the coefficient alfa we can determine the
643
corresponding change delta x[k] = new x[k] - old x[k],
644
where old x[k] is a value of x[k] in the current basis,
645
and new x[k] is a value of x[k] in the adjacent basis */
646
delta_k = delta_j / alfa;
647
/* Tomlin noticed that if the variable x[k] is of integer
648
kind, its change cannot be less (eventually) than one in
650
if (k > m && tree->int_col[k-m])
651
{ /* x[k] is structural integer variable */
652
if (gnm_abs(delta_k - gnm_floor(delta_k + 0.5)) > 1e-3)
654
delta_k = gnm_ceil(delta_k); /* +3.14 -> +4 */
656
delta_k = gnm_floor(delta_k); /* -3.14 -> -4 */
659
/* now determine the status and reduced cost of x[k] in the
662
{ stat = lpx_get_row_stat(lp, k);
663
dk = lpx_get_row_dual(lp, k);
666
{ stat = lpx_get_col_stat(lp, k-m);
667
dk = lpx_get_col_dual(lp, k-m);
669
/* if the current basis is dual degenerative, some reduced
670
costs which are close to zero may have wrong sign due to
671
round-off errors, so correct the sign of d[k] */
674
if (stat == LPX_NL && dk < 0.0 ||
675
stat == LPX_NU && dk > 0.0 ||
676
stat == LPX_NF) dk = 0.0;
679
if (stat == LPX_NL && dk > 0.0 ||
680
stat == LPX_NU && dk < 0.0 ||
681
stat == LPX_NF) dk = 0.0;
684
insist(tree != tree);
686
/* now knowing the change of x[k] and its reduced cost d[k]
687
we can compute the corresponding change in the objective
688
function delta Z = new Z - old Z = d[k] * delta x[k];
689
note that due to Tomlin's modification new Z can be even
690
worse than in the adjacent basis */
691
delta_z = dk * delta_k;
692
skip: /* new Z is never better than old Z, therefore the change
693
delta Z is always non-negative (in case of minimization)
694
or non-positive (in case of maximization) */
696
{ case LPX_MIN: insist(delta_z >= 0.0); break;
697
case LPX_MAX: insist(delta_z <= 0.0); break;
698
default: insist(tree != tree);
700
/* save the change in the objective fnction for down- and
701
up-branches, respectively */
702
if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;
704
/* thus, in down-branch no integer feasible solution can be
705
better than Z + dz_dn, and in up-branch no integer feasible
706
solution can be better than Z + dz_up, where Z is value of
707
the objective function in the current basis */
708
/* following the heuristic by Driebeck and Tomlin we choose a
709
column (i.e. structural variable) which provides largest
710
degradation of the objective function in some of branches;
711
besides, we select the branch with smaller degradation to
712
be solved next and keep other branch with larger degradation
713
in the active list hoping to minimize the number of further
715
if (degrad < gnm_abs(dz_dn) || degrad < gnm_abs(dz_up))
717
if (gnm_abs(dz_dn) < gnm_abs(dz_up))
718
{ /* select down branch to be solved next */
720
degrad = gnm_abs(dz_up);
723
{ /* select up branch to be solved next */
725
degrad = gnm_abs(dz_dn);
727
/* save the objective changes for printing */
728
dd_dn = dz_dn, dd_up = dz_up;
729
/* if down- or up-branch has no feasible solution, we does
730
not need to consider other candidates (in principle, the
731
corresponding branch could be pruned right now) */
732
if (degrad == DBL_MAX) break;
735
/* free working arrays */
738
/* something must be chosen */
739
insist(1 <= jj && jj <= n);
740
if (tree->msg_lev >= 3)
741
{ print("branch_drtom: column %d chosen to branch on", jj);
742
if (gnm_abs(dd_dn) == DBL_MAX)
743
print("branch_drtom: down-branch is infeasible");
745
print("branch_drtom: down-branch bound is %.9e",
746
lpx_get_obj_val(lp) + dd_dn);
747
if (gnm_abs(dd_up) == DBL_MAX)
748
print("branch_drtom: up-branch is infeasible");
750
print("branch_drtom: up-branch bound is %.9e",
751
lpx_get_obj_val(lp) + dd_up);
753
/* perform branching */
754
branch_on(tree, jj, next);
758
/*----------------------------------------------------------------------
759
-- branch_mostf - choose most fractional branching variable.
761
-- This routine looks up the list of structural variables and chooses
762
-- that one, which is of integer kind and has most fractional value in
763
-- optimal solution of the current LP relaxation.
765
-- This routine also selects the branch to be solved next where integer
766
-- infeasibility of the chosen variable is less than in other one. */
768
static void branch_mostf(MIPTREE *tree)
772
gnm_float beta, most, temp;
773
/* choose the column to branch on */
774
jj = 0, most = DBL_MAX;
775
for (j = 1; j <= n; j++)
776
{ if (tree->non_int[j])
777
{ beta = lpx_get_col_prim(lp, j);
778
temp = gnm_floor(beta) + 0.5;
779
if (most > gnm_abs(beta - temp))
780
{ jj = j, most = gnm_abs(beta - temp);
782
next = -1; /* down branch */
784
next = +1; /* up branch */
788
/* perform branching */
789
branch_on(tree, jj, next);
793
/*----------------------------------------------------------------------
794
-- cleanup_the_tree - prune hopeless branches from the tree.
796
-- This routine walks through the active list and checks the local
797
-- bound for every active subproblem. If the local bound indicates that
798
-- the subproblem cannot have integer optimal solution better than the
799
-- incumbent objective value, the routine deletes such subproblem that,
800
-- in turn, involves pruning the corresponding branch of the tree. */
802
static void cleanup_the_tree(MIPTREE *tree)
803
{ MIPNODE *node, *next_node;
805
/* the global bound must exist */
807
/* walk through the list of active subproblems */
808
for (node = tree->head; node != NULL; node = next_node)
809
{ /* deleting some active problem node may involve deleting its
810
parents recursively; however, all its parents being created
811
*before* it are always *precede* it in the node list, so
812
the next problem node is never affected by such deletion */
813
next_node = node->next;
814
/* if the branch is hopeless, prune it */
815
if (!is_branch_hopeful(tree, node->p))
816
mip_delete_node(tree, node->p), count++;
818
if (tree->msg_lev >= 3)
820
print("One hopeless branch has been pruned");
822
print("%d hopeless branches have been pruned", count);
827
/*----------------------------------------------------------------------
828
-- btrack_most_feas - select "most integer feasible" subproblem.
830
-- This routine selects from the active list a subproblem to be solved
831
-- next, whose parent has minimal sum of integer infeasibilities. */
833
static void btrack_most_feas(MIPTREE *tree)
837
p = 0, best = DBL_MAX;
838
for (node = tree->head; node != NULL; node = node->next)
839
{ insist(node->up != NULL);
840
if (best > node->up->ii_sum)
841
p = node->p, best = node->up->ii_sum;
843
mip_revive_node(tree, p);
847
/*----------------------------------------------------------------------
848
-- btrack_best_proj - select subproblem with best projection heur.
850
-- This routine selects from the active list a subproblem to be solved
851
-- next using the best projection heuristic. */
853
static void btrack_best_proj(MIPTREE *tree)
854
{ MIPNODE *root, *node;
856
gnm_float best, deg, obj;
857
/* the global bound must exist */
859
/* obtain pointer to the root node, which must exist */
860
root = tree->slot[1].node;
861
insist(root != NULL);
862
/* deg estimates degradation of the objective function per unit
863
of the sum of integer infeasibilities */
864
insist(root->ii_sum > 0.0);
865
deg = (tree->best - root->bound) / root->ii_sum;
866
/* nothing has been selected so far */
867
p = 0, best = DBL_MAX;
868
/* walk through the list of active subproblems */
869
for (node = tree->head; node != NULL; node = node->next)
870
{ insist(node->up != NULL);
871
/* obj estimates optimal objective value if the sum of integer
872
infeasibilities were zero */
873
obj = node->up->bound + deg * node->up->ii_sum;
874
if (tree->dir == LPX_MAX) obj = - obj;
875
/* select the subproblem which has the best estimated optimal
877
if (best > obj) p = node->p, best = obj;
879
mip_revive_node(tree, p);
883
/*----------------------------------------------------------------------
884
-- mip_driver - branch-and-bound driver.
888
-- #include "glpmip.h"
889
-- int mip_driver(MIPTREE *tree);
893
-- The routine mip_driver is a branch-and-bound driver that manages the
894
-- process of solving MIP problem instance.
898
-- The routine mip_driver returns one of the following exit codes:
900
-- MIP_E_OK the search is finished;
902
-- MIP_E_ITLIM iterations limit exceeded;
904
-- MIP_E_TMLIM time limit exceeded;
906
-- MIP_E_ERROR an error occurred on solving the LP relaxation of the
907
-- current subproblem.
909
-- Should note that additional exit codes may appear in future versions
910
-- of this routine. */
912
int mip_driver(MIPTREE *tree)
913
{ LPX *lp = tree->lp;
914
int p, p_stat, d_stat, ret;
915
/* revive the root subproblem */
916
mip_revive_node(tree, 1);
917
loop: /* main loop starts here; at this point some subproblem has been
918
just selected from the active list and made current */
919
insist(tree->curr != NULL);
920
/* determine the reference number of the current subproblem */
922
if (tree->msg_lev >= 3)
923
{ int level = tree->slot[p].node->level;
924
print("-------------------------------------------------------"
925
"-----------------");
926
print("Processing node %d at level %d", p, level);
928
/* check if the time limit has been exhausted */
929
if (tree->tm_lim >= 0.0 && tree->tm_lim <= utime() - tree->tm_beg)
930
{ if (tree->msg_lev >= 3)
931
print("Time limit exceeded; search terminated");
935
/* display current progress of the search */
936
if (tree->msg_lev >= 3 || tree->msg_lev >= 2 &&
937
utime() - tree->tm_lag >= tree->out_frq - 0.001)
939
/* solve LP relaxation of the current subproblem */
940
if (tree->msg_lev >= 3)
941
print("Solving LP relaxation...");
942
ret = mip_solve_node(tree);
949
if (tree->msg_lev >= 3)
950
print("Iterations limit exceeded; search terminated");
954
if (tree->msg_lev >= 1)
955
print("mip_driver: cannot solve current LP relaxation");
959
/* analyze status of the basic solution */
960
p_stat = lpx_get_prim_stat(lp);
961
d_stat = lpx_get_dual_stat(lp);
962
if (p_stat == LPX_P_FEAS && d_stat == LPX_D_FEAS)
963
{ /* LP relaxation has optimal solution */
964
if (tree->msg_lev >= 3)
965
print("Found optimal solution to LP relaxation");
967
else if (p_stat == LPX_P_FEAS && d_stat == LPX_D_NOFEAS)
968
{ /* LP relaxation has unbounded solution */
969
/* since the current subproblem cannot have a larger feasible
970
region than its parent, there is something wrong */
971
if (tree->msg_lev >= 1)
972
print("mip_driver: current LP relaxation has unbounded solu"
977
else if (p_stat == LPX_P_INFEAS && d_stat == LPX_D_FEAS)
978
{ /* LP relaxation has no primal solution which is better than
979
the incumbent objective value */
981
if (tree->msg_lev >= 3)
982
print("LP relaxation has no solution better than incumbent "
984
/* prune the branch */
987
else if (p_stat == LPX_P_NOFEAS)
988
{ /* LP relaxation has no primal feasible solution */
989
if (tree->msg_lev >= 3)
990
print("LP relaxation has no feasible solution");
991
/* prune the branch */
995
{ /* other cases cannot appear */
998
/* at this point basic solution of LP relaxation of the current
999
subproblem is optimal */
1000
insist(p_stat == LPX_P_FEAS && d_stat == LPX_D_FEAS);
1001
/* thus, it defines a local bound for integer optimal solution of
1002
the current subproblem */
1003
set_local_bound(tree);
1004
/* if the local bound indicates that integer optimal solution of
1005
the current subproblem cannot be better than the global bound,
1007
if (!is_branch_hopeful(tree, p))
1008
{ if (tree->msg_lev >= 3)
1009
print("Current branch is hopeless and can be pruned");
1012
/* check integrality of the basic solution */
1013
check_integrality(tree);
1014
/* if the basic solution satisfies to all integrality conditions,
1015
it is a new, better integer feasible solution */
1016
if (tree->curr->ii_cnt == 0)
1017
{ if (tree->msg_lev >= 3)
1018
print("New integer feasible solution found");
1019
record_solution(tree);
1020
if (tree->msg_lev >= 2) show_progress(tree);
1021
/* the current subproblem is fathomed; prune its branch */
1024
/* basic solution of LP relaxation of the current subproblem is
1025
integer infeasible */
1026
/* try to fix some non-basic structural variables of integer kind
1027
on their current bounds using reduced costs */
1028
if (tree->found) fix_by_red_cost(tree);
1029
/* perform branching */
1030
switch (tree->branch)
1032
/* branch on first appropriate variable */
1036
/* branch on last appropriate variable */
1040
/* branch using the heuristic by Dreebeck and Tomlin */
1044
/* branch on most fractional variable */
1048
insist(tree != tree);
1050
/* continue the search from the subproblem which begins down- or
1051
up-branch (it has been revived by branching routine) */
1053
fath: /* the current subproblem has been fathomed */
1054
if (tree->msg_lev >= 3)
1055
print("Node %d fathomed", p);
1056
/* freeze the current subproblem */
1057
mip_freeze_node(tree);
1058
/* and prune the corresponding branch of the tree */
1059
mip_delete_node(tree, p);
1060
/* if a new integer feasible solution has just been found, other
1061
branches may become hopeless and therefore should be pruned */
1062
if (tree->found) cleanup_the_tree(tree);
1063
/* if the active list is empty, the search is finished */
1064
if (tree->head == NULL)
1065
{ if (tree->msg_lev >= 3)
1066
print("Active list is empty!");
1067
insist(tree->node_pool->count == 0);
1068
insist(tree->bnds_pool->count == 0);
1069
insist(tree->stat_pool->count == 0);
1073
/* perform backtracking */
1074
switch (tree->btrack)
1076
/* depth first search */
1077
insist(tree->tail != NULL);
1078
mip_revive_node(tree, tree->tail->p);
1081
/* breadth first search */
1082
insist(tree->head != NULL);
1083
mip_revive_node(tree, tree->head->p);
1087
{ /* "most integer feasible" subproblem */
1088
btrack_most_feas(tree);
1091
{ /* best projection heuristic */
1092
btrack_best_proj(tree);
1096
/* select node with best local bound */
1097
mip_revive_node(tree, mip_best_node(tree));
1100
insist(tree != tree);
1102
/* continue the search from the subproblem selected */
1104
done: /* display status of the search on exit from the solver */
1105
if (tree->msg_lev >= 2) show_progress(tree);
1106
/* decrease the time limit by spent amount of time */
1107
if (tree->tm_lim >= 0.0)
1108
{ tree->tm_lim -= (utime() - tree->tm_beg);
1109
if (tree->tm_lim < 0.0) tree->tm_lim = 0.0;
1111
/* return to the calling program */
added
removed
src/tools/solver/glpk/source/glpmip2.c
