|
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 */ |