1 | /* glpios10.c (feasibility pump heuristic) */ |
---|
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 | #include "glprng.h" |
---|
27 | |
---|
28 | /*********************************************************************** |
---|
29 | * NAME |
---|
30 | * |
---|
31 | * ios_feas_pump - feasibility pump heuristic |
---|
32 | * |
---|
33 | * SYNOPSIS |
---|
34 | * |
---|
35 | * #include "glpios.h" |
---|
36 | * void ios_feas_pump(glp_tree *T); |
---|
37 | * |
---|
38 | * DESCRIPTION |
---|
39 | * |
---|
40 | * The routine ios_feas_pump is a simple implementation of the Feasi- |
---|
41 | * bility Pump heuristic. |
---|
42 | * |
---|
43 | * REFERENCES |
---|
44 | * |
---|
45 | * M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math. |
---|
46 | * Program., Ser. A 104, pp. 91-104 (2005). */ |
---|
47 | |
---|
48 | struct VAR |
---|
49 | { /* binary variable */ |
---|
50 | int j; |
---|
51 | /* ordinal number */ |
---|
52 | int x; |
---|
53 | /* value in the rounded solution (0 or 1) */ |
---|
54 | double d; |
---|
55 | /* sorting key */ |
---|
56 | }; |
---|
57 | |
---|
58 | static int fcmp(const void *x, const void *y) |
---|
59 | { /* comparison routine */ |
---|
60 | const struct VAR *vx = x, *vy = y; |
---|
61 | if (vx->d > vy->d) |
---|
62 | return -1; |
---|
63 | else if (vx->d < vy->d) |
---|
64 | return +1; |
---|
65 | else |
---|
66 | return 0; |
---|
67 | } |
---|
68 | |
---|
69 | void ios_feas_pump(glp_tree *T) |
---|
70 | { glp_prob *P = T->mip; |
---|
71 | int n = P->n; |
---|
72 | glp_prob *lp = NULL; |
---|
73 | struct VAR *var = NULL; |
---|
74 | RNG *rand = NULL; |
---|
75 | GLPCOL *col; |
---|
76 | glp_smcp parm; |
---|
77 | int j, k, new_x, nfail, npass, nv, ret, stalling; |
---|
78 | double dist, tol; |
---|
79 | xassert(glp_get_status(P) == GLP_OPT); |
---|
80 | /* this heuristic is applied only once on the root level */ |
---|
81 | if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done; |
---|
82 | /* determine number of binary variables */ |
---|
83 | nv = 0; |
---|
84 | for (j = 1; j <= n; j++) |
---|
85 | { col = P->col[j]; |
---|
86 | /* if x[j] is continuous, skip it */ |
---|
87 | if (col->kind == GLP_CV) continue; |
---|
88 | /* if x[j] is fixed, skip it */ |
---|
89 | if (col->type == GLP_FX) continue; |
---|
90 | /* x[j] is non-fixed integer */ |
---|
91 | xassert(col->kind == GLP_IV); |
---|
92 | if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0) |
---|
93 | { /* x[j] is binary */ |
---|
94 | nv++; |
---|
95 | } |
---|
96 | else |
---|
97 | { /* x[j] is general integer */ |
---|
98 | if (T->parm->msg_lev >= GLP_MSG_ALL) |
---|
99 | xprintf("FPUMP heuristic cannot be applied due to genera" |
---|
100 | "l integer variables\n"); |
---|
101 | goto done; |
---|
102 | } |
---|
103 | } |
---|
104 | /* there must be at least one binary variable */ |
---|
105 | if (nv == 0) goto done; |
---|
106 | if (T->parm->msg_lev >= GLP_MSG_ALL) |
---|
107 | xprintf("Applying FPUMP heuristic...\n"); |
---|
108 | /* build the list of binary variables */ |
---|
109 | var = xcalloc(1+nv, sizeof(struct VAR)); |
---|
110 | k = 0; |
---|
111 | for (j = 1; j <= n; j++) |
---|
112 | { col = P->col[j]; |
---|
113 | if (col->kind == GLP_IV && col->type == GLP_DB) |
---|
114 | var[++k].j = j; |
---|
115 | } |
---|
116 | xassert(k == nv); |
---|
117 | /* create working problem object */ |
---|
118 | lp = glp_create_prob(); |
---|
119 | more: /* copy the original problem object to keep it intact */ |
---|
120 | glp_copy_prob(lp, P, GLP_OFF); |
---|
121 | /* we are interested to find an integer feasible solution, which |
---|
122 | is better than the best known one */ |
---|
123 | if (P->mip_stat == GLP_FEAS) |
---|
124 | { int *ind; |
---|
125 | double *val, bnd; |
---|
126 | /* add a row and make it identical to the objective row */ |
---|
127 | glp_add_rows(lp, 1); |
---|
128 | ind = xcalloc(1+n, sizeof(int)); |
---|
129 | val = xcalloc(1+n, sizeof(double)); |
---|
130 | for (j = 1; j <= n; j++) |
---|
131 | { ind[j] = j; |
---|
132 | val[j] = P->col[j]->coef; |
---|
133 | } |
---|
134 | glp_set_mat_row(lp, lp->m, n, ind, val); |
---|
135 | xfree(ind); |
---|
136 | xfree(val); |
---|
137 | /* introduce upper (minimization) or lower (maximization) |
---|
138 | bound to the original objective function; note that this |
---|
139 | additional constraint is not violated at the optimal point |
---|
140 | to LP relaxation */ |
---|
141 | #if 0 /* modified by xypron <xypron.glpk@gmx.de> */ |
---|
142 | if (P->dir == GLP_MIN) |
---|
143 | { bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj)); |
---|
144 | if (bnd < P->obj_val) bnd = P->obj_val; |
---|
145 | glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0); |
---|
146 | } |
---|
147 | else if (P->dir == GLP_MAX) |
---|
148 | { bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj)); |
---|
149 | if (bnd > P->obj_val) bnd = P->obj_val; |
---|
150 | glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0); |
---|
151 | } |
---|
152 | else |
---|
153 | xassert(P != P); |
---|
154 | #else |
---|
155 | bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj; |
---|
156 | /* xprintf("bnd = %f\n", bnd); */ |
---|
157 | if (P->dir == GLP_MIN) |
---|
158 | glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0); |
---|
159 | else if (P->dir == GLP_MAX) |
---|
160 | glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0); |
---|
161 | else |
---|
162 | xassert(P != P); |
---|
163 | #endif |
---|
164 | } |
---|
165 | /* reset pass count */ |
---|
166 | npass = 0; |
---|
167 | /* invalidate the rounded point */ |
---|
168 | for (k = 1; k <= nv; k++) |
---|
169 | var[k].x = -1; |
---|
170 | pass: /* next pass starts here */ |
---|
171 | npass++; |
---|
172 | if (T->parm->msg_lev >= GLP_MSG_ALL) |
---|
173 | xprintf("Pass %d\n", npass); |
---|
174 | /* initialize minimal distance between the basic point and the |
---|
175 | rounded one obtained during this pass */ |
---|
176 | dist = DBL_MAX; |
---|
177 | /* reset failure count (the number of succeeded iterations failed |
---|
178 | to improve the distance) */ |
---|
179 | nfail = 0; |
---|
180 | /* if it is not the first pass, perturb the last rounded point |
---|
181 | rather than construct it from the basic solution */ |
---|
182 | if (npass > 1) |
---|
183 | { double rho, temp; |
---|
184 | if (rand == NULL) |
---|
185 | rand = rng_create_rand(); |
---|
186 | for (k = 1; k <= nv; k++) |
---|
187 | { j = var[k].j; |
---|
188 | col = lp->col[j]; |
---|
189 | rho = rng_uniform(rand, -0.3, 0.7); |
---|
190 | if (rho < 0.0) rho = 0.0; |
---|
191 | temp = fabs((double)var[k].x - col->prim); |
---|
192 | if (temp + rho > 0.5) var[k].x = 1 - var[k].x; |
---|
193 | } |
---|
194 | goto skip; |
---|
195 | } |
---|
196 | loop: /* innermost loop begins here */ |
---|
197 | /* round basic solution (which is assumed primal feasible) */ |
---|
198 | stalling = 1; |
---|
199 | for (k = 1; k <= nv; k++) |
---|
200 | { col = lp->col[var[k].j]; |
---|
201 | if (col->prim < 0.5) |
---|
202 | { /* rounded value is 0 */ |
---|
203 | new_x = 0; |
---|
204 | } |
---|
205 | else |
---|
206 | { /* rounded value is 1 */ |
---|
207 | new_x = 1; |
---|
208 | } |
---|
209 | if (var[k].x != new_x) |
---|
210 | { stalling = 0; |
---|
211 | var[k].x = new_x; |
---|
212 | } |
---|
213 | } |
---|
214 | /* if the rounded point has not changed (stalling), choose and |
---|
215 | flip some its entries heuristically */ |
---|
216 | if (stalling) |
---|
217 | { /* compute d[j] = |x[j] - round(x[j])| */ |
---|
218 | for (k = 1; k <= nv; k++) |
---|
219 | { col = lp->col[var[k].j]; |
---|
220 | var[k].d = fabs(col->prim - (double)var[k].x); |
---|
221 | } |
---|
222 | /* sort the list of binary variables by descending d[j] */ |
---|
223 | qsort(&var[1], nv, sizeof(struct VAR), fcmp); |
---|
224 | /* choose and flip some rounded components */ |
---|
225 | for (k = 1; k <= nv; k++) |
---|
226 | { if (k >= 5 && var[k].d < 0.35 || k >= 10) break; |
---|
227 | var[k].x = 1 - var[k].x; |
---|
228 | } |
---|
229 | } |
---|
230 | skip: /* check if the time limit has been exhausted */ |
---|
231 | if (T->parm->tm_lim < INT_MAX && |
---|
232 | (double)(T->parm->tm_lim - 1) <= |
---|
233 | 1000.0 * xdifftime(xtime(), T->tm_beg)) goto done; |
---|
234 | /* build the objective, which is the distance between the current |
---|
235 | (basic) point and the rounded one */ |
---|
236 | lp->dir = GLP_MIN; |
---|
237 | lp->c0 = 0.0; |
---|
238 | for (j = 1; j <= n; j++) |
---|
239 | lp->col[j]->coef = 0.0; |
---|
240 | for (k = 1; k <= nv; k++) |
---|
241 | { j = var[k].j; |
---|
242 | if (var[k].x == 0) |
---|
243 | lp->col[j]->coef = +1.0; |
---|
244 | else |
---|
245 | { lp->col[j]->coef = -1.0; |
---|
246 | lp->c0 += 1.0; |
---|
247 | } |
---|
248 | } |
---|
249 | /* minimize the distance with the simplex method */ |
---|
250 | glp_init_smcp(&parm); |
---|
251 | if (T->parm->msg_lev <= GLP_MSG_ERR) |
---|
252 | parm.msg_lev = T->parm->msg_lev; |
---|
253 | else if (T->parm->msg_lev <= GLP_MSG_ALL) |
---|
254 | { parm.msg_lev = GLP_MSG_ON; |
---|
255 | parm.out_dly = 10000; |
---|
256 | } |
---|
257 | ret = glp_simplex(lp, &parm); |
---|
258 | if (ret != 0) |
---|
259 | { if (T->parm->msg_lev >= GLP_MSG_ERR) |
---|
260 | xprintf("Warning: glp_simplex returned %d\n", ret); |
---|
261 | goto done; |
---|
262 | } |
---|
263 | ret = glp_get_status(lp); |
---|
264 | if (ret != GLP_OPT) |
---|
265 | { if (T->parm->msg_lev >= GLP_MSG_ERR) |
---|
266 | xprintf("Warning: glp_get_status returned %d\n", ret); |
---|
267 | goto done; |
---|
268 | } |
---|
269 | if (T->parm->msg_lev >= GLP_MSG_DBG) |
---|
270 | xprintf("delta = %g\n", lp->obj_val); |
---|
271 | /* check if the basic solution is integer feasible; note that it |
---|
272 | may be so even if the minimial distance is positive */ |
---|
273 | tol = 0.3 * T->parm->tol_int; |
---|
274 | for (k = 1; k <= nv; k++) |
---|
275 | { col = lp->col[var[k].j]; |
---|
276 | if (tol < col->prim && col->prim < 1.0 - tol) break; |
---|
277 | } |
---|
278 | if (k > nv) |
---|
279 | { /* okay; the basic solution seems to be integer feasible */ |
---|
280 | double *x = xcalloc(1+n, sizeof(double)); |
---|
281 | for (j = 1; j <= n; j++) |
---|
282 | { x[j] = lp->col[j]->prim; |
---|
283 | if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5); |
---|
284 | } |
---|
285 | #if 1 /* modified by xypron <xypron.glpk@gmx.de> */ |
---|
286 | /* reset direction and right-hand side of objective */ |
---|
287 | lp->c0 = P->c0; |
---|
288 | lp->dir = P->dir; |
---|
289 | /* fix integer variables */ |
---|
290 | for (k = 1; k <= nv; k++) |
---|
291 | { lp->col[var[k].j]->lb = x[var[k].j]; |
---|
292 | lp->col[var[k].j]->ub = x[var[k].j]; |
---|
293 | lp->col[var[k].j]->type = GLP_FX; |
---|
294 | } |
---|
295 | /* copy original objective function */ |
---|
296 | for (j = 1; j <= n; j++) |
---|
297 | lp->col[j]->coef = P->col[j]->coef; |
---|
298 | /* solve original LP and copy result */ |
---|
299 | ret = glp_simplex(lp, &parm); |
---|
300 | if (ret != 0) |
---|
301 | { if (T->parm->msg_lev >= GLP_MSG_ERR) |
---|
302 | xprintf("Warning: glp_simplex returned %d\n", ret); |
---|
303 | goto done; |
---|
304 | } |
---|
305 | ret = glp_get_status(lp); |
---|
306 | if (ret != GLP_OPT) |
---|
307 | { if (T->parm->msg_lev >= GLP_MSG_ERR) |
---|
308 | xprintf("Warning: glp_get_status returned %d\n", ret); |
---|
309 | goto done; |
---|
310 | } |
---|
311 | for (j = 1; j <= n; j++) |
---|
312 | if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim; |
---|
313 | #endif |
---|
314 | ret = glp_ios_heur_sol(T, x); |
---|
315 | xfree(x); |
---|
316 | if (ret == 0) |
---|
317 | { /* the integer solution is accepted */ |
---|
318 | if (ios_is_hopeful(T, T->curr->bound)) |
---|
319 | { /* it is reasonable to apply the heuristic once again */ |
---|
320 | goto more; |
---|
321 | } |
---|
322 | else |
---|
323 | { /* the best known integer feasible solution just found |
---|
324 | is close to optimal solution to LP relaxation */ |
---|
325 | goto done; |
---|
326 | } |
---|
327 | } |
---|
328 | } |
---|
329 | /* the basic solution is fractional */ |
---|
330 | if (dist == DBL_MAX || |
---|
331 | lp->obj_val <= dist - 1e-6 * (1.0 + dist)) |
---|
332 | { /* the distance is reducing */ |
---|
333 | nfail = 0, dist = lp->obj_val; |
---|
334 | } |
---|
335 | else |
---|
336 | { /* improving the distance failed */ |
---|
337 | nfail++; |
---|
338 | } |
---|
339 | if (nfail < 3) goto loop; |
---|
340 | if (npass < 5) goto pass; |
---|
341 | done: /* delete working objects */ |
---|
342 | if (lp != NULL) glp_delete_prob(lp); |
---|
343 | if (var != NULL) xfree(var); |
---|
344 | if (rand != NULL) rng_delete_rand(rand); |
---|
345 | return; |
---|
346 | } |
---|
347 | |
---|
348 | /* eof */ |
---|