[1] | 1 | /* glpios09.c (branching heuristics) */ |
---|
| 2 | |
---|
| 3 | /*********************************************************************** |
---|
| 4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
| 5 | * |
---|
| 6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
| 7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
| 8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
| 9 | * E-mail: <mao@gnu.org>. |
---|
| 10 | * |
---|
| 11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
| 12 | * under the terms of the GNU General Public License as published by |
---|
| 13 | * the Free Software Foundation, either version 3 of the License, or |
---|
| 14 | * (at your option) any later version. |
---|
| 15 | * |
---|
| 16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
| 17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
| 18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
| 19 | * License for more details. |
---|
| 20 | * |
---|
| 21 | * You should have received a copy of the GNU General Public License |
---|
| 22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
| 23 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #include "glpios.h" |
---|
| 26 | |
---|
| 27 | /*********************************************************************** |
---|
| 28 | * NAME |
---|
| 29 | * |
---|
| 30 | * ios_choose_var - select variable to branch on |
---|
| 31 | * |
---|
| 32 | * SYNOPSIS |
---|
| 33 | * |
---|
| 34 | * #include "glpios.h" |
---|
| 35 | * int ios_choose_var(glp_tree *T, int *next); |
---|
| 36 | * |
---|
| 37 | * The routine ios_choose_var chooses a variable from the candidate |
---|
| 38 | * list to branch on. Additionally the routine provides a flag stored |
---|
| 39 | * in the location next to suggests which of the child subproblems |
---|
| 40 | * should be solved next. |
---|
| 41 | * |
---|
| 42 | * RETURNS |
---|
| 43 | * |
---|
| 44 | * The routine ios_choose_var returns the ordinal number of the column |
---|
| 45 | * choosen. */ |
---|
| 46 | |
---|
| 47 | static int branch_first(glp_tree *T, int *next); |
---|
| 48 | static int branch_last(glp_tree *T, int *next); |
---|
| 49 | static int branch_mostf(glp_tree *T, int *next); |
---|
| 50 | static int branch_drtom(glp_tree *T, int *next); |
---|
| 51 | |
---|
| 52 | int ios_choose_var(glp_tree *T, int *next) |
---|
| 53 | { int j; |
---|
| 54 | if (T->parm->br_tech == GLP_BR_FFV) |
---|
| 55 | { /* branch on first fractional variable */ |
---|
| 56 | j = branch_first(T, next); |
---|
| 57 | } |
---|
| 58 | else if (T->parm->br_tech == GLP_BR_LFV) |
---|
| 59 | { /* branch on last fractional variable */ |
---|
| 60 | j = branch_last(T, next); |
---|
| 61 | } |
---|
| 62 | else if (T->parm->br_tech == GLP_BR_MFV) |
---|
| 63 | { /* branch on most fractional variable */ |
---|
| 64 | j = branch_mostf(T, next); |
---|
| 65 | } |
---|
| 66 | else if (T->parm->br_tech == GLP_BR_DTH) |
---|
| 67 | { /* branch using the heuristic by Dreebeck and Tomlin */ |
---|
| 68 | j = branch_drtom(T, next); |
---|
| 69 | } |
---|
| 70 | else if (T->parm->br_tech == GLP_BR_PCH) |
---|
| 71 | { /* hybrid pseudocost heuristic */ |
---|
| 72 | j = ios_pcost_branch(T, next); |
---|
| 73 | } |
---|
| 74 | else |
---|
| 75 | xassert(T != T); |
---|
| 76 | return j; |
---|
| 77 | } |
---|
| 78 | |
---|
| 79 | /*********************************************************************** |
---|
| 80 | * branch_first - choose first branching variable |
---|
| 81 | * |
---|
| 82 | * This routine looks up the list of structural variables and chooses |
---|
| 83 | * the first one, which is of integer kind and has fractional value in |
---|
| 84 | * optimal solution to the current LP relaxation. |
---|
| 85 | * |
---|
| 86 | * This routine also selects the branch to be solved next where integer |
---|
| 87 | * infeasibility of the chosen variable is less than in other one. */ |
---|
| 88 | |
---|
| 89 | static int branch_first(glp_tree *T, int *_next) |
---|
| 90 | { int j, next; |
---|
| 91 | double beta; |
---|
| 92 | /* choose the column to branch on */ |
---|
| 93 | for (j = 1; j <= T->n; j++) |
---|
| 94 | if (T->non_int[j]) break; |
---|
| 95 | xassert(1 <= j && j <= T->n); |
---|
| 96 | /* select the branch to be solved next */ |
---|
| 97 | beta = glp_get_col_prim(T->mip, j); |
---|
| 98 | if (beta - floor(beta) < ceil(beta) - beta) |
---|
| 99 | next = GLP_DN_BRNCH; |
---|
| 100 | else |
---|
| 101 | next = GLP_UP_BRNCH; |
---|
| 102 | *_next = next; |
---|
| 103 | return j; |
---|
| 104 | } |
---|
| 105 | |
---|
| 106 | /*********************************************************************** |
---|
| 107 | * branch_last - choose last branching variable |
---|
| 108 | * |
---|
| 109 | * This routine looks up the list of structural variables and chooses |
---|
| 110 | * the last one, which is of integer kind and has fractional value in |
---|
| 111 | * optimal solution to the current LP relaxation. |
---|
| 112 | * |
---|
| 113 | * This routine also selects the branch to be solved next where integer |
---|
| 114 | * infeasibility of the chosen variable is less than in other one. */ |
---|
| 115 | |
---|
| 116 | static int branch_last(glp_tree *T, int *_next) |
---|
| 117 | { int j, next; |
---|
| 118 | double beta; |
---|
| 119 | /* choose the column to branch on */ |
---|
| 120 | for (j = T->n; j >= 1; j--) |
---|
| 121 | if (T->non_int[j]) break; |
---|
| 122 | xassert(1 <= j && j <= T->n); |
---|
| 123 | /* select the branch to be solved next */ |
---|
| 124 | beta = glp_get_col_prim(T->mip, j); |
---|
| 125 | if (beta - floor(beta) < ceil(beta) - beta) |
---|
| 126 | next = GLP_DN_BRNCH; |
---|
| 127 | else |
---|
| 128 | next = GLP_UP_BRNCH; |
---|
| 129 | *_next = next; |
---|
| 130 | return j; |
---|
| 131 | } |
---|
| 132 | |
---|
| 133 | /*********************************************************************** |
---|
| 134 | * branch_mostf - choose most fractional branching variable |
---|
| 135 | * |
---|
| 136 | * This routine looks up the list of structural variables and chooses |
---|
| 137 | * that one, which is of integer kind and has most fractional value in |
---|
| 138 | * optimal solution to the current LP relaxation. |
---|
| 139 | * |
---|
| 140 | * This routine also selects the branch to be solved next where integer |
---|
| 141 | * infeasibility of the chosen variable is less than in other one. |
---|
| 142 | * |
---|
| 143 | * (Alexander Martin notices that "...most infeasible is as good as |
---|
| 144 | * random...".) */ |
---|
| 145 | |
---|
| 146 | static int branch_mostf(glp_tree *T, int *_next) |
---|
| 147 | { int j, jj, next; |
---|
| 148 | double beta, most, temp; |
---|
| 149 | /* choose the column to branch on */ |
---|
| 150 | jj = 0, most = DBL_MAX; |
---|
| 151 | for (j = 1; j <= T->n; j++) |
---|
| 152 | { if (T->non_int[j]) |
---|
| 153 | { beta = glp_get_col_prim(T->mip, j); |
---|
| 154 | temp = floor(beta) + 0.5; |
---|
| 155 | if (most > fabs(beta - temp)) |
---|
| 156 | { jj = j, most = fabs(beta - temp); |
---|
| 157 | if (beta < temp) |
---|
| 158 | next = GLP_DN_BRNCH; |
---|
| 159 | else |
---|
| 160 | next = GLP_UP_BRNCH; |
---|
| 161 | } |
---|
| 162 | } |
---|
| 163 | } |
---|
| 164 | *_next = next; |
---|
| 165 | return jj; |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | /*********************************************************************** |
---|
| 169 | * branch_drtom - choose branching var using Driebeck-Tomlin heuristic |
---|
| 170 | * |
---|
| 171 | * This routine chooses a structural variable, which is required to be |
---|
| 172 | * integral and has fractional value in optimal solution of the current |
---|
| 173 | * LP relaxation, using a heuristic proposed by Driebeck and Tomlin. |
---|
| 174 | * |
---|
| 175 | * The routine also selects the branch to be solved next, again due to |
---|
| 176 | * Driebeck and Tomlin. |
---|
| 177 | * |
---|
| 178 | * This routine is based on the heuristic proposed in: |
---|
| 179 | * |
---|
| 180 | * Driebeck N.J. An algorithm for the solution of mixed-integer |
---|
| 181 | * programming problems, Management Science, 12: 576-87 (1966); |
---|
| 182 | * |
---|
| 183 | * and improved in: |
---|
| 184 | * |
---|
| 185 | * Tomlin J.A. Branch and bound methods for integer and non-convex |
---|
| 186 | * programming, in J.Abadie (ed.), Integer and Nonlinear Programming, |
---|
| 187 | * North-Holland, Amsterdam, pp. 437-50 (1970). |
---|
| 188 | * |
---|
| 189 | * Must note that this heuristic is time-expensive, because computing |
---|
| 190 | * one-step degradation (see the routine below) requires one BTRAN for |
---|
| 191 | * each fractional-valued structural variable. */ |
---|
| 192 | |
---|
| 193 | static int branch_drtom(glp_tree *T, int *_next) |
---|
| 194 | { glp_prob *mip = T->mip; |
---|
| 195 | int m = mip->m; |
---|
| 196 | int n = mip->n; |
---|
| 197 | char *non_int = T->non_int; |
---|
| 198 | int j, jj, k, t, next, kase, len, stat, *ind; |
---|
| 199 | double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up, |
---|
| 200 | dd_dn, dd_up, degrad, *val; |
---|
| 201 | /* basic solution of LP relaxation must be optimal */ |
---|
| 202 | xassert(glp_get_status(mip) == GLP_OPT); |
---|
| 203 | /* allocate working arrays */ |
---|
| 204 | ind = xcalloc(1+n, sizeof(int)); |
---|
| 205 | val = xcalloc(1+n, sizeof(double)); |
---|
| 206 | /* nothing has been chosen so far */ |
---|
| 207 | jj = 0, degrad = -1.0; |
---|
| 208 | /* walk through the list of columns (structural variables) */ |
---|
| 209 | for (j = 1; j <= n; j++) |
---|
| 210 | { /* if j-th column is not marked as fractional, skip it */ |
---|
| 211 | if (!non_int[j]) continue; |
---|
| 212 | /* obtain (fractional) value of j-th column in basic solution |
---|
| 213 | of LP relaxation */ |
---|
| 214 | x = glp_get_col_prim(mip, j); |
---|
| 215 | /* since the value of j-th column is fractional, the column is |
---|
| 216 | basic; compute corresponding row of the simplex table */ |
---|
| 217 | len = glp_eval_tab_row(mip, m+j, ind, val); |
---|
| 218 | /* the following fragment computes a change in the objective |
---|
| 219 | function: delta Z = new Z - old Z, where old Z is the |
---|
| 220 | objective value in the current optimal basis, and new Z is |
---|
| 221 | the objective value in the adjacent basis, for two cases: |
---|
| 222 | 1) if new upper bound ub' = floor(x[j]) is introduced for |
---|
| 223 | j-th column (down branch); |
---|
| 224 | 2) if new lower bound lb' = ceil(x[j]) is introduced for |
---|
| 225 | j-th column (up branch); |
---|
| 226 | since in both cases the solution remaining dual feasible |
---|
| 227 | becomes primal infeasible, one implicit simplex iteration |
---|
| 228 | is performed to determine the change delta Z; |
---|
| 229 | it is obvious that new Z, which is never better than old Z, |
---|
| 230 | is a lower (minimization) or upper (maximization) bound of |
---|
| 231 | the objective function for down- and up-branches. */ |
---|
| 232 | for (kase = -1; kase <= +1; kase += 2) |
---|
| 233 | { /* if kase < 0, the new upper bound of x[j] is introduced; |
---|
| 234 | in this case x[j] should decrease in order to leave the |
---|
| 235 | basis and go to its new upper bound */ |
---|
| 236 | /* if kase > 0, the new lower bound of x[j] is introduced; |
---|
| 237 | in this case x[j] should increase in order to leave the |
---|
| 238 | basis and go to its new lower bound */ |
---|
| 239 | /* apply the dual ratio test in order to determine which |
---|
| 240 | auxiliary or structural variable should enter the basis |
---|
| 241 | to keep dual feasibility */ |
---|
| 242 | k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9); |
---|
| 243 | if (k != 0) k = ind[k]; |
---|
| 244 | /* if no non-basic variable has been chosen, LP relaxation |
---|
| 245 | of corresponding branch being primal infeasible and dual |
---|
| 246 | unbounded has no primal feasible solution; in this case |
---|
| 247 | the change delta Z is formally set to infinity */ |
---|
| 248 | if (k == 0) |
---|
| 249 | { delta_z = |
---|
| 250 | (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX); |
---|
| 251 | goto skip; |
---|
| 252 | } |
---|
| 253 | /* row of the simplex table that corresponds to non-basic |
---|
| 254 | variable x[k] choosen by the dual ratio test is: |
---|
| 255 | x[j] = ... + alfa * x[k] + ... |
---|
| 256 | where alfa is the influence coefficient (an element of |
---|
| 257 | the simplex table row) */ |
---|
| 258 | /* determine the coefficient alfa */ |
---|
| 259 | for (t = 1; t <= len; t++) if (ind[t] == k) break; |
---|
| 260 | xassert(1 <= t && t <= len); |
---|
| 261 | alfa = val[t]; |
---|
| 262 | /* since in the adjacent basis the variable x[j] becomes |
---|
| 263 | non-basic, knowing its value in the current basis we can |
---|
| 264 | determine its change delta x[j] = new x[j] - old x[j] */ |
---|
| 265 | delta_j = (kase < 0 ? floor(x) : ceil(x)) - x; |
---|
| 266 | /* and knowing the coefficient alfa we can determine the |
---|
| 267 | corresponding change delta x[k] = new x[k] - old x[k], |
---|
| 268 | where old x[k] is a value of x[k] in the current basis, |
---|
| 269 | and new x[k] is a value of x[k] in the adjacent basis */ |
---|
| 270 | delta_k = delta_j / alfa; |
---|
| 271 | /* Tomlin noticed that if the variable x[k] is of integer |
---|
| 272 | kind, its change cannot be less (eventually) than one in |
---|
| 273 | the magnitude */ |
---|
| 274 | if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV) |
---|
| 275 | { /* x[k] is structural integer variable */ |
---|
| 276 | if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3) |
---|
| 277 | { if (delta_k > 0.0) |
---|
| 278 | delta_k = ceil(delta_k); /* +3.14 -> +4 */ |
---|
| 279 | else |
---|
| 280 | delta_k = floor(delta_k); /* -3.14 -> -4 */ |
---|
| 281 | } |
---|
| 282 | } |
---|
| 283 | /* now determine the status and reduced cost of x[k] in the |
---|
| 284 | current basis */ |
---|
| 285 | if (k <= m) |
---|
| 286 | { stat = glp_get_row_stat(mip, k); |
---|
| 287 | dk = glp_get_row_dual(mip, k); |
---|
| 288 | } |
---|
| 289 | else |
---|
| 290 | { stat = glp_get_col_stat(mip, k-m); |
---|
| 291 | dk = glp_get_col_dual(mip, k-m); |
---|
| 292 | } |
---|
| 293 | /* if the current basis is dual degenerate, some reduced |
---|
| 294 | costs which are close to zero may have wrong sign due to |
---|
| 295 | round-off errors, so correct the sign of d[k] */ |
---|
| 296 | switch (T->mip->dir) |
---|
| 297 | { case GLP_MIN: |
---|
| 298 | if (stat == GLP_NL && dk < 0.0 || |
---|
| 299 | stat == GLP_NU && dk > 0.0 || |
---|
| 300 | stat == GLP_NF) dk = 0.0; |
---|
| 301 | break; |
---|
| 302 | case GLP_MAX: |
---|
| 303 | if (stat == GLP_NL && dk > 0.0 || |
---|
| 304 | stat == GLP_NU && dk < 0.0 || |
---|
| 305 | stat == GLP_NF) dk = 0.0; |
---|
| 306 | break; |
---|
| 307 | default: |
---|
| 308 | xassert(T != T); |
---|
| 309 | } |
---|
| 310 | /* now knowing the change of x[k] and its reduced cost d[k] |
---|
| 311 | we can compute the corresponding change in the objective |
---|
| 312 | function delta Z = new Z - old Z = d[k] * delta x[k]; |
---|
| 313 | note that due to Tomlin's modification new Z can be even |
---|
| 314 | worse than in the adjacent basis */ |
---|
| 315 | delta_z = dk * delta_k; |
---|
| 316 | skip: /* new Z is never better than old Z, therefore the change |
---|
| 317 | delta Z is always non-negative (in case of minimization) |
---|
| 318 | or non-positive (in case of maximization) */ |
---|
| 319 | switch (T->mip->dir) |
---|
| 320 | { case GLP_MIN: xassert(delta_z >= 0.0); break; |
---|
| 321 | case GLP_MAX: xassert(delta_z <= 0.0); break; |
---|
| 322 | default: xassert(T != T); |
---|
| 323 | } |
---|
| 324 | /* save the change in the objective fnction for down- and |
---|
| 325 | up-branches, respectively */ |
---|
| 326 | if (kase < 0) dz_dn = delta_z; else dz_up = delta_z; |
---|
| 327 | } |
---|
| 328 | /* thus, in down-branch no integer feasible solution can be |
---|
| 329 | better than Z + dz_dn, and in up-branch no integer feasible |
---|
| 330 | solution can be better than Z + dz_up, where Z is value of |
---|
| 331 | the objective function in the current basis */ |
---|
| 332 | /* following the heuristic by Driebeck and Tomlin we choose a |
---|
| 333 | column (i.e. structural variable) which provides largest |
---|
| 334 | degradation of the objective function in some of branches; |
---|
| 335 | besides, we select the branch with smaller degradation to |
---|
| 336 | be solved next and keep other branch with larger degradation |
---|
| 337 | in the active list hoping to minimize the number of further |
---|
| 338 | backtrackings */ |
---|
| 339 | if (degrad < fabs(dz_dn) || degrad < fabs(dz_up)) |
---|
| 340 | { jj = j; |
---|
| 341 | if (fabs(dz_dn) < fabs(dz_up)) |
---|
| 342 | { /* select down branch to be solved next */ |
---|
| 343 | next = GLP_DN_BRNCH; |
---|
| 344 | degrad = fabs(dz_up); |
---|
| 345 | } |
---|
| 346 | else |
---|
| 347 | { /* select up branch to be solved next */ |
---|
| 348 | next = GLP_UP_BRNCH; |
---|
| 349 | degrad = fabs(dz_dn); |
---|
| 350 | } |
---|
| 351 | /* save the objective changes for printing */ |
---|
| 352 | dd_dn = dz_dn, dd_up = dz_up; |
---|
| 353 | /* if down- or up-branch has no feasible solution, we does |
---|
| 354 | not need to consider other candidates (in principle, the |
---|
| 355 | corresponding branch could be pruned right now) */ |
---|
| 356 | if (degrad == DBL_MAX) break; |
---|
| 357 | } |
---|
| 358 | } |
---|
| 359 | /* free working arrays */ |
---|
| 360 | xfree(ind); |
---|
| 361 | xfree(val); |
---|
| 362 | /* something must be chosen */ |
---|
| 363 | xassert(1 <= jj && jj <= n); |
---|
| 364 | #if 1 /* 02/XI-2009 */ |
---|
| 365 | if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val))) |
---|
| 366 | { jj = branch_mostf(T, &next); |
---|
| 367 | goto done; |
---|
| 368 | } |
---|
| 369 | #endif |
---|
| 370 | if (T->parm->msg_lev >= GLP_MSG_DBG) |
---|
| 371 | { xprintf("branch_drtom: column %d chosen to branch on\n", jj); |
---|
| 372 | if (fabs(dd_dn) == DBL_MAX) |
---|
| 373 | xprintf("branch_drtom: down-branch is infeasible\n"); |
---|
| 374 | else |
---|
| 375 | xprintf("branch_drtom: down-branch bound is %.9e\n", |
---|
| 376 | lpx_get_obj_val(mip) + dd_dn); |
---|
| 377 | if (fabs(dd_up) == DBL_MAX) |
---|
| 378 | xprintf("branch_drtom: up-branch is infeasible\n"); |
---|
| 379 | else |
---|
| 380 | xprintf("branch_drtom: up-branch bound is %.9e\n", |
---|
| 381 | lpx_get_obj_val(mip) + dd_up); |
---|
| 382 | } |
---|
| 383 | done: *_next = next; |
---|
| 384 | return jj; |
---|
| 385 | } |
---|
| 386 | |
---|
| 387 | /**********************************************************************/ |
---|
| 388 | |
---|
| 389 | struct csa |
---|
| 390 | { /* common storage area */ |
---|
| 391 | int *dn_cnt; /* int dn_cnt[1+n]; */ |
---|
| 392 | /* dn_cnt[j] is the number of subproblems, whose LP relaxations |
---|
| 393 | have been solved and which are down-branches for variable x[j]; |
---|
| 394 | dn_cnt[j] = 0 means the down pseudocost is uninitialized */ |
---|
| 395 | double *dn_sum; /* double dn_sum[1+n]; */ |
---|
| 396 | /* dn_sum[j] is the sum of per unit degradations of the objective |
---|
| 397 | over all dn_cnt[j] subproblems */ |
---|
| 398 | int *up_cnt; /* int up_cnt[1+n]; */ |
---|
| 399 | /* up_cnt[j] is the number of subproblems, whose LP relaxations |
---|
| 400 | have been solved and which are up-branches for variable x[j]; |
---|
| 401 | up_cnt[j] = 0 means the up pseudocost is uninitialized */ |
---|
| 402 | double *up_sum; /* double up_sum[1+n]; */ |
---|
| 403 | /* up_sum[j] is the sum of per unit degradations of the objective |
---|
| 404 | over all up_cnt[j] subproblems */ |
---|
| 405 | }; |
---|
| 406 | |
---|
| 407 | void *ios_pcost_init(glp_tree *tree) |
---|
| 408 | { /* initialize working data used on pseudocost branching */ |
---|
| 409 | struct csa *csa; |
---|
| 410 | int n = tree->n, j; |
---|
| 411 | csa = xmalloc(sizeof(struct csa)); |
---|
| 412 | csa->dn_cnt = xcalloc(1+n, sizeof(int)); |
---|
| 413 | csa->dn_sum = xcalloc(1+n, sizeof(double)); |
---|
| 414 | csa->up_cnt = xcalloc(1+n, sizeof(int)); |
---|
| 415 | csa->up_sum = xcalloc(1+n, sizeof(double)); |
---|
| 416 | for (j = 1; j <= n; j++) |
---|
| 417 | { csa->dn_cnt[j] = csa->up_cnt[j] = 0; |
---|
| 418 | csa->dn_sum[j] = csa->up_sum[j] = 0.0; |
---|
| 419 | } |
---|
| 420 | return csa; |
---|
| 421 | } |
---|
| 422 | |
---|
| 423 | static double eval_degrad(glp_prob *P, int j, double bnd) |
---|
| 424 | { /* compute degradation of the objective on fixing x[j] at given |
---|
| 425 | value with a limited number of dual simplex iterations */ |
---|
| 426 | /* this routine fixes column x[j] at specified value bnd, |
---|
| 427 | solves resulting LP, and returns a lower bound to degradation |
---|
| 428 | of the objective, degrad >= 0 */ |
---|
| 429 | glp_prob *lp; |
---|
| 430 | glp_smcp parm; |
---|
| 431 | int ret; |
---|
| 432 | double degrad; |
---|
| 433 | /* the current basis must be optimal */ |
---|
| 434 | xassert(glp_get_status(P) == GLP_OPT); |
---|
| 435 | /* create a copy of P */ |
---|
| 436 | lp = glp_create_prob(); |
---|
| 437 | glp_copy_prob(lp, P, 0); |
---|
| 438 | /* fix column x[j] at specified value */ |
---|
| 439 | glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd); |
---|
| 440 | /* try to solve resulting LP */ |
---|
| 441 | glp_init_smcp(&parm); |
---|
| 442 | parm.msg_lev = GLP_MSG_OFF; |
---|
| 443 | parm.meth = GLP_DUAL; |
---|
| 444 | parm.it_lim = 30; |
---|
| 445 | parm.out_dly = 1000; |
---|
| 446 | parm.meth = GLP_DUAL; |
---|
| 447 | ret = glp_simplex(lp, &parm); |
---|
| 448 | if (ret == 0 || ret == GLP_EITLIM) |
---|
| 449 | { if (glp_get_prim_stat(lp) == GLP_NOFEAS) |
---|
| 450 | { /* resulting LP has no primal feasible solution */ |
---|
| 451 | degrad = DBL_MAX; |
---|
| 452 | } |
---|
| 453 | else if (glp_get_dual_stat(lp) == GLP_FEAS) |
---|
| 454 | { /* resulting basis is optimal or at least dual feasible, |
---|
| 455 | so we have the correct lower bound to degradation */ |
---|
| 456 | if (P->dir == GLP_MIN) |
---|
| 457 | degrad = lp->obj_val - P->obj_val; |
---|
| 458 | else if (P->dir == GLP_MAX) |
---|
| 459 | degrad = P->obj_val - lp->obj_val; |
---|
| 460 | else |
---|
| 461 | xassert(P != P); |
---|
| 462 | /* degradation cannot be negative by definition */ |
---|
| 463 | /* note that the lower bound to degradation may be close |
---|
| 464 | to zero even if its exact value is zero due to round-off |
---|
| 465 | errors on computing the objective value */ |
---|
| 466 | if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val))) |
---|
| 467 | degrad = 0.0; |
---|
| 468 | } |
---|
| 469 | else |
---|
| 470 | { /* the final basis reported by the simplex solver is dual |
---|
| 471 | infeasible, so we cannot determine a non-trivial lower |
---|
| 472 | bound to degradation */ |
---|
| 473 | degrad = 0.0; |
---|
| 474 | } |
---|
| 475 | } |
---|
| 476 | else |
---|
| 477 | { /* the simplex solver failed */ |
---|
| 478 | degrad = 0.0; |
---|
| 479 | } |
---|
| 480 | /* delete the copy of P */ |
---|
| 481 | glp_delete_prob(lp); |
---|
| 482 | return degrad; |
---|
| 483 | } |
---|
| 484 | |
---|
| 485 | void ios_pcost_update(glp_tree *tree) |
---|
| 486 | { /* update history information for pseudocost branching */ |
---|
| 487 | /* this routine is called every time when LP relaxation of the |
---|
| 488 | current subproblem has been solved to optimality with all lazy |
---|
| 489 | and cutting plane constraints included */ |
---|
| 490 | int j; |
---|
| 491 | double dx, dz, psi; |
---|
| 492 | struct csa *csa = tree->pcost; |
---|
| 493 | xassert(csa != NULL); |
---|
| 494 | xassert(tree->curr != NULL); |
---|
| 495 | /* if the current subproblem is the root, skip updating */ |
---|
| 496 | if (tree->curr->up == NULL) goto skip; |
---|
| 497 | /* determine branching variable x[j], which was used in the |
---|
| 498 | parent subproblem to create the current subproblem */ |
---|
| 499 | j = tree->curr->up->br_var; |
---|
| 500 | xassert(1 <= j && j <= tree->n); |
---|
| 501 | /* determine the change dx[j] = new x[j] - old x[j], |
---|
| 502 | where new x[j] is a value of x[j] in optimal solution to LP |
---|
| 503 | relaxation of the current subproblem, old x[j] is a value of |
---|
| 504 | x[j] in optimal solution to LP relaxation of the parent |
---|
| 505 | subproblem */ |
---|
| 506 | dx = tree->mip->col[j]->prim - tree->curr->up->br_val; |
---|
| 507 | xassert(dx != 0.0); |
---|
| 508 | /* determine corresponding change dz = new dz - old dz in the |
---|
| 509 | objective function value */ |
---|
| 510 | dz = tree->mip->obj_val - tree->curr->up->lp_obj; |
---|
| 511 | /* determine per unit degradation of the objective function */ |
---|
| 512 | psi = fabs(dz / dx); |
---|
| 513 | /* update history information */ |
---|
| 514 | if (dx < 0.0) |
---|
| 515 | { /* the current subproblem is down-branch */ |
---|
| 516 | csa->dn_cnt[j]++; |
---|
| 517 | csa->dn_sum[j] += psi; |
---|
| 518 | } |
---|
| 519 | else /* dx > 0.0 */ |
---|
| 520 | { /* the current subproblem is up-branch */ |
---|
| 521 | csa->up_cnt[j]++; |
---|
| 522 | csa->up_sum[j] += psi; |
---|
| 523 | } |
---|
| 524 | skip: return; |
---|
| 525 | } |
---|
| 526 | |
---|
| 527 | void ios_pcost_free(glp_tree *tree) |
---|
| 528 | { /* free working area used on pseudocost branching */ |
---|
| 529 | struct csa *csa = tree->pcost; |
---|
| 530 | xassert(csa != NULL); |
---|
| 531 | xfree(csa->dn_cnt); |
---|
| 532 | xfree(csa->dn_sum); |
---|
| 533 | xfree(csa->up_cnt); |
---|
| 534 | xfree(csa->up_sum); |
---|
| 535 | xfree(csa); |
---|
| 536 | tree->pcost = NULL; |
---|
| 537 | return; |
---|
| 538 | } |
---|
| 539 | |
---|
| 540 | static double eval_psi(glp_tree *T, int j, int brnch) |
---|
| 541 | { /* compute estimation of pseudocost of variable x[j] for down- |
---|
| 542 | or up-branch */ |
---|
| 543 | struct csa *csa = T->pcost; |
---|
| 544 | double beta, degrad, psi; |
---|
| 545 | xassert(csa != NULL); |
---|
| 546 | xassert(1 <= j && j <= T->n); |
---|
| 547 | if (brnch == GLP_DN_BRNCH) |
---|
| 548 | { /* down-branch */ |
---|
| 549 | if (csa->dn_cnt[j] == 0) |
---|
| 550 | { /* initialize down pseudocost */ |
---|
| 551 | beta = T->mip->col[j]->prim; |
---|
| 552 | degrad = eval_degrad(T->mip, j, floor(beta)); |
---|
| 553 | if (degrad == DBL_MAX) |
---|
| 554 | { psi = DBL_MAX; |
---|
| 555 | goto done; |
---|
| 556 | } |
---|
| 557 | csa->dn_cnt[j] = 1; |
---|
| 558 | csa->dn_sum[j] = degrad / (beta - floor(beta)); |
---|
| 559 | } |
---|
| 560 | psi = csa->dn_sum[j] / (double)csa->dn_cnt[j]; |
---|
| 561 | } |
---|
| 562 | else if (brnch == GLP_UP_BRNCH) |
---|
| 563 | { /* up-branch */ |
---|
| 564 | if (csa->up_cnt[j] == 0) |
---|
| 565 | { /* initialize up pseudocost */ |
---|
| 566 | beta = T->mip->col[j]->prim; |
---|
| 567 | degrad = eval_degrad(T->mip, j, ceil(beta)); |
---|
| 568 | if (degrad == DBL_MAX) |
---|
| 569 | { psi = DBL_MAX; |
---|
| 570 | goto done; |
---|
| 571 | } |
---|
| 572 | csa->up_cnt[j] = 1; |
---|
| 573 | csa->up_sum[j] = degrad / (ceil(beta) - beta); |
---|
| 574 | } |
---|
| 575 | psi = csa->up_sum[j] / (double)csa->up_cnt[j]; |
---|
| 576 | } |
---|
| 577 | else |
---|
| 578 | xassert(brnch != brnch); |
---|
| 579 | done: return psi; |
---|
| 580 | } |
---|
| 581 | |
---|
| 582 | static void progress(glp_tree *T) |
---|
| 583 | { /* display progress of pseudocost initialization */ |
---|
| 584 | struct csa *csa = T->pcost; |
---|
| 585 | int j, nv = 0, ni = 0; |
---|
| 586 | for (j = 1; j <= T->n; j++) |
---|
| 587 | { if (glp_ios_can_branch(T, j)) |
---|
| 588 | { nv++; |
---|
| 589 | if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++; |
---|
| 590 | } |
---|
| 591 | } |
---|
| 592 | xprintf("Pseudocosts initialized for %d of %d variables\n", |
---|
| 593 | ni, nv); |
---|
| 594 | return; |
---|
| 595 | } |
---|
| 596 | |
---|
| 597 | int ios_pcost_branch(glp_tree *T, int *_next) |
---|
| 598 | { /* choose branching variable with pseudocost branching */ |
---|
| 599 | glp_long t = xtime(); |
---|
| 600 | int j, jjj, sel; |
---|
| 601 | double beta, psi, d1, d2, d, dmax; |
---|
| 602 | /* initialize the working arrays */ |
---|
| 603 | if (T->pcost == NULL) |
---|
| 604 | T->pcost = ios_pcost_init(T); |
---|
| 605 | /* nothing has been chosen so far */ |
---|
| 606 | jjj = 0, dmax = -1.0; |
---|
| 607 | /* go through the list of branching candidates */ |
---|
| 608 | for (j = 1; j <= T->n; j++) |
---|
| 609 | { if (!glp_ios_can_branch(T, j)) continue; |
---|
| 610 | /* determine primal value of x[j] in optimal solution to LP |
---|
| 611 | relaxation of the current subproblem */ |
---|
| 612 | beta = T->mip->col[j]->prim; |
---|
| 613 | /* estimate pseudocost of x[j] for down-branch */ |
---|
| 614 | psi = eval_psi(T, j, GLP_DN_BRNCH); |
---|
| 615 | if (psi == DBL_MAX) |
---|
| 616 | { /* down-branch has no primal feasible solution */ |
---|
| 617 | jjj = j, sel = GLP_DN_BRNCH; |
---|
| 618 | goto done; |
---|
| 619 | } |
---|
| 620 | /* estimate degradation of the objective for down-branch */ |
---|
| 621 | d1 = psi * (beta - floor(beta)); |
---|
| 622 | /* estimate pseudocost of x[j] for up-branch */ |
---|
| 623 | psi = eval_psi(T, j, GLP_UP_BRNCH); |
---|
| 624 | if (psi == DBL_MAX) |
---|
| 625 | { /* up-branch has no primal feasible solution */ |
---|
| 626 | jjj = j, sel = GLP_UP_BRNCH; |
---|
| 627 | goto done; |
---|
| 628 | } |
---|
| 629 | /* estimate degradation of the objective for up-branch */ |
---|
| 630 | d2 = psi * (ceil(beta) - beta); |
---|
| 631 | /* determine d = max(d1, d2) */ |
---|
| 632 | d = (d1 > d2 ? d1 : d2); |
---|
| 633 | /* choose x[j] which provides maximal estimated degradation of |
---|
| 634 | the objective either in down- or up-branch */ |
---|
| 635 | if (dmax < d) |
---|
| 636 | { dmax = d; |
---|
| 637 | jjj = j; |
---|
| 638 | /* continue the search from a subproblem, where degradation |
---|
| 639 | is less than in other one */ |
---|
| 640 | sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH); |
---|
| 641 | } |
---|
| 642 | /* display progress of pseudocost initialization */ |
---|
| 643 | if (T->parm->msg_lev >= GLP_ON) |
---|
| 644 | { if (xdifftime(xtime(), t) >= 10.0) |
---|
| 645 | { progress(T); |
---|
| 646 | t = xtime(); |
---|
| 647 | } |
---|
| 648 | } |
---|
| 649 | } |
---|
| 650 | if (dmax == 0.0) |
---|
| 651 | { /* no degradation is indicated; choose a variable having most |
---|
| 652 | fractional value */ |
---|
| 653 | jjj = branch_mostf(T, &sel); |
---|
| 654 | } |
---|
| 655 | done: *_next = sel; |
---|
| 656 | return jjj; |
---|
| 657 | } |
---|
| 658 | |
---|
| 659 | /* eof */ |
---|