lemon-project-template-glpk
comparison deps/glpk/src/glpios07.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1fb72ceeccd1 |
---|---|
1 /* glpios07.c (mixed cover cut generator) */ | |
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, 2011 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 -- COVER INEQUALITIES | |
29 -- | |
30 -- Consider the set of feasible solutions to 0-1 knapsack problem: | |
31 -- | |
32 -- sum a[j]*x[j] <= b, (1) | |
33 -- j in J | |
34 -- | |
35 -- x[j] is binary, (2) | |
36 -- | |
37 -- where, wlog, we assume that a[j] > 0 (since 0-1 variables can be | |
38 -- complemented) and a[j] <= b (since a[j] > b implies x[j] = 0). | |
39 -- | |
40 -- A set C within J is called a cover if | |
41 -- | |
42 -- sum a[j] > b. (3) | |
43 -- j in C | |
44 -- | |
45 -- For any cover C the inequality | |
46 -- | |
47 -- sum x[j] <= |C| - 1 (4) | |
48 -- j in C | |
49 -- | |
50 -- is called a cover inequality and is valid for (1)-(2). | |
51 -- | |
52 -- MIXED COVER INEQUALITIES | |
53 -- | |
54 -- Consider the set of feasible solutions to mixed knapsack problem: | |
55 -- | |
56 -- sum a[j]*x[j] + y <= b, (5) | |
57 -- j in J | |
58 -- | |
59 -- x[j] is binary, (6) | |
60 -- | |
61 -- 0 <= y <= u is continuous, (7) | |
62 -- | |
63 -- where again we assume that a[j] > 0. | |
64 -- | |
65 -- Let C within J be some set. From (1)-(4) it follows that | |
66 -- | |
67 -- sum a[j] > b - y (8) | |
68 -- j in C | |
69 -- | |
70 -- implies | |
71 -- | |
72 -- sum x[j] <= |C| - 1. (9) | |
73 -- j in C | |
74 -- | |
75 -- Thus, we need to modify the inequality (9) in such a way that it be | |
76 -- a constraint only if the condition (8) is satisfied. | |
77 -- | |
78 -- Consider the following inequality: | |
79 -- | |
80 -- sum x[j] <= |C| - t. (10) | |
81 -- j in C | |
82 -- | |
83 -- If 0 < t <= 1, then (10) is equivalent to (9), because all x[j] are | |
84 -- binary variables. On the other hand, if t <= 0, (10) being satisfied | |
85 -- for any values of x[j] is not a constraint. | |
86 -- | |
87 -- Let | |
88 -- | |
89 -- t' = sum a[j] + y - b. (11) | |
90 -- j in C | |
91 -- | |
92 -- It is understood that the condition t' > 0 is equivalent to (8). | |
93 -- Besides, from (6)-(7) it follows that t' has an implied upper bound: | |
94 -- | |
95 -- t'max = sum a[j] + u - b. (12) | |
96 -- j in C | |
97 -- | |
98 -- This allows to express the parameter t having desired properties: | |
99 -- | |
100 -- t = t' / t'max. (13) | |
101 -- | |
102 -- In fact, t <= 1 by definition, and t > 0 being equivalent to t' > 0 | |
103 -- is equivalent to (8). | |
104 -- | |
105 -- Thus, the inequality (10), where t is given by formula (13) is valid | |
106 -- for (5)-(7). | |
107 -- | |
108 -- Note that if u = 0, then y = 0, so t = 1, and the conditions (8) and | |
109 -- (10) is transformed to the conditions (3) and (4). | |
110 -- | |
111 -- GENERATING MIXED COVER CUTS | |
112 -- | |
113 -- To generate a mixed cover cut in the form (10) we need to find such | |
114 -- set C which satisfies to the inequality (8) and for which, in turn, | |
115 -- the inequality (10) is violated in the current point. | |
116 -- | |
117 -- Substituting t from (13) to (10) gives: | |
118 -- | |
119 -- 1 | |
120 -- sum x[j] <= |C| - ----- (sum a[j] + y - b), (14) | |
121 -- j in C t'max j in C | |
122 -- | |
123 -- and finally we have the cut inequality in the standard form: | |
124 -- | |
125 -- sum x[j] + alfa * y <= beta, (15) | |
126 -- j in C | |
127 -- | |
128 -- where: | |
129 -- | |
130 -- alfa = 1 / t'max, (16) | |
131 -- | |
132 -- beta = |C| - alfa * (sum a[j] - b). (17) | |
133 -- j in C */ | |
134 | |
135 #if 1 | |
136 #define MAXTRY 1000 | |
137 #else | |
138 #define MAXTRY 10000 | |
139 #endif | |
140 | |
141 static int cover2(int n, double a[], double b, double u, double x[], | |
142 double y, int cov[], double *_alfa, double *_beta) | |
143 { /* try to generate mixed cover cut using two-element cover */ | |
144 int i, j, try = 0, ret = 0; | |
145 double eps, alfa, beta, temp, rmax = 0.001; | |
146 eps = 0.001 * (1.0 + fabs(b)); | |
147 for (i = 0+1; i <= n; i++) | |
148 for (j = i+1; j <= n; j++) | |
149 { /* C = {i, j} */ | |
150 try++; | |
151 if (try > MAXTRY) goto done; | |
152 /* check if condition (8) is satisfied */ | |
153 if (a[i] + a[j] + y > b + eps) | |
154 { /* compute parameters for inequality (15) */ | |
155 temp = a[i] + a[j] - b; | |
156 alfa = 1.0 / (temp + u); | |
157 beta = 2.0 - alfa * temp; | |
158 /* compute violation of inequality (15) */ | |
159 temp = x[i] + x[j] + alfa * y - beta; | |
160 /* choose C providing maximum violation */ | |
161 if (rmax < temp) | |
162 { rmax = temp; | |
163 cov[1] = i; | |
164 cov[2] = j; | |
165 *_alfa = alfa; | |
166 *_beta = beta; | |
167 ret = 1; | |
168 } | |
169 } | |
170 } | |
171 done: return ret; | |
172 } | |
173 | |
174 static int cover3(int n, double a[], double b, double u, double x[], | |
175 double y, int cov[], double *_alfa, double *_beta) | |
176 { /* try to generate mixed cover cut using three-element cover */ | |
177 int i, j, k, try = 0, ret = 0; | |
178 double eps, alfa, beta, temp, rmax = 0.001; | |
179 eps = 0.001 * (1.0 + fabs(b)); | |
180 for (i = 0+1; i <= n; i++) | |
181 for (j = i+1; j <= n; j++) | |
182 for (k = j+1; k <= n; k++) | |
183 { /* C = {i, j, k} */ | |
184 try++; | |
185 if (try > MAXTRY) goto done; | |
186 /* check if condition (8) is satisfied */ | |
187 if (a[i] + a[j] + a[k] + y > b + eps) | |
188 { /* compute parameters for inequality (15) */ | |
189 temp = a[i] + a[j] + a[k] - b; | |
190 alfa = 1.0 / (temp + u); | |
191 beta = 3.0 - alfa * temp; | |
192 /* compute violation of inequality (15) */ | |
193 temp = x[i] + x[j] + x[k] + alfa * y - beta; | |
194 /* choose C providing maximum violation */ | |
195 if (rmax < temp) | |
196 { rmax = temp; | |
197 cov[1] = i; | |
198 cov[2] = j; | |
199 cov[3] = k; | |
200 *_alfa = alfa; | |
201 *_beta = beta; | |
202 ret = 1; | |
203 } | |
204 } | |
205 } | |
206 done: return ret; | |
207 } | |
208 | |
209 static int cover4(int n, double a[], double b, double u, double x[], | |
210 double y, int cov[], double *_alfa, double *_beta) | |
211 { /* try to generate mixed cover cut using four-element cover */ | |
212 int i, j, k, l, try = 0, ret = 0; | |
213 double eps, alfa, beta, temp, rmax = 0.001; | |
214 eps = 0.001 * (1.0 + fabs(b)); | |
215 for (i = 0+1; i <= n; i++) | |
216 for (j = i+1; j <= n; j++) | |
217 for (k = j+1; k <= n; k++) | |
218 for (l = k+1; l <= n; l++) | |
219 { /* C = {i, j, k, l} */ | |
220 try++; | |
221 if (try > MAXTRY) goto done; | |
222 /* check if condition (8) is satisfied */ | |
223 if (a[i] + a[j] + a[k] + a[l] + y > b + eps) | |
224 { /* compute parameters for inequality (15) */ | |
225 temp = a[i] + a[j] + a[k] + a[l] - b; | |
226 alfa = 1.0 / (temp + u); | |
227 beta = 4.0 - alfa * temp; | |
228 /* compute violation of inequality (15) */ | |
229 temp = x[i] + x[j] + x[k] + x[l] + alfa * y - beta; | |
230 /* choose C providing maximum violation */ | |
231 if (rmax < temp) | |
232 { rmax = temp; | |
233 cov[1] = i; | |
234 cov[2] = j; | |
235 cov[3] = k; | |
236 cov[4] = l; | |
237 *_alfa = alfa; | |
238 *_beta = beta; | |
239 ret = 1; | |
240 } | |
241 } | |
242 } | |
243 done: return ret; | |
244 } | |
245 | |
246 static int cover(int n, double a[], double b, double u, double x[], | |
247 double y, int cov[], double *alfa, double *beta) | |
248 { /* try to generate mixed cover cut; | |
249 input (see (5)): | |
250 n is the number of binary variables; | |
251 a[1:n] are coefficients at binary variables; | |
252 b is the right-hand side; | |
253 u is upper bound of continuous variable; | |
254 x[1:n] are values of binary variables at current point; | |
255 y is value of continuous variable at current point; | |
256 output (see (15), (16), (17)): | |
257 cov[1:r] are indices of binary variables included in cover C, | |
258 where r is the set cardinality returned on exit; | |
259 alfa coefficient at continuous variable; | |
260 beta is the right-hand side; */ | |
261 int j; | |
262 /* perform some sanity checks */ | |
263 xassert(n >= 2); | |
264 for (j = 1; j <= n; j++) xassert(a[j] > 0.0); | |
265 #if 1 /* ??? */ | |
266 xassert(b > -1e-5); | |
267 #else | |
268 xassert(b > 0.0); | |
269 #endif | |
270 xassert(u >= 0.0); | |
271 for (j = 1; j <= n; j++) xassert(0.0 <= x[j] && x[j] <= 1.0); | |
272 xassert(0.0 <= y && y <= u); | |
273 /* try to generate mixed cover cut */ | |
274 if (cover2(n, a, b, u, x, y, cov, alfa, beta)) return 2; | |
275 if (cover3(n, a, b, u, x, y, cov, alfa, beta)) return 3; | |
276 if (cover4(n, a, b, u, x, y, cov, alfa, beta)) return 4; | |
277 return 0; | |
278 } | |
279 | |
280 /*---------------------------------------------------------------------- | |
281 -- lpx_cover_cut - generate mixed cover cut. | |
282 -- | |
283 -- SYNOPSIS | |
284 -- | |
285 -- #include "glplpx.h" | |
286 -- int lpx_cover_cut(LPX *lp, int len, int ind[], double val[], | |
287 -- double work[]); | |
288 -- | |
289 -- DESCRIPTION | |
290 -- | |
291 -- The routine lpx_cover_cut generates a mixed cover cut for a given | |
292 -- row of the MIP problem. | |
293 -- | |
294 -- The given row of the MIP problem should be explicitly specified in | |
295 -- the form: | |
296 -- | |
297 -- sum{j in J} a[j]*x[j] <= b. (1) | |
298 -- | |
299 -- On entry indices (ordinal numbers) of structural variables, which | |
300 -- have non-zero constraint coefficients, should be placed in locations | |
301 -- ind[1], ..., ind[len], and corresponding constraint coefficients | |
302 -- should be placed in locations val[1], ..., val[len]. The right-hand | |
303 -- side b should be stored in location val[0]. | |
304 -- | |
305 -- The working array work should have at least nb locations, where nb | |
306 -- is the number of binary variables in (1). | |
307 -- | |
308 -- The routine generates a mixed cover cut in the same form as (1) and | |
309 -- stores the cut coefficients and right-hand side in the same way as | |
310 -- just described above. | |
311 -- | |
312 -- RETURNS | |
313 -- | |
314 -- If the cutting plane has been successfully generated, the routine | |
315 -- returns 1 <= len' <= n, which is the number of non-zero coefficients | |
316 -- in the inequality constraint. Otherwise, the routine returns zero. */ | |
317 | |
318 static int lpx_cover_cut(LPX *lp, int len, int ind[], double val[], | |
319 double work[]) | |
320 { int cov[1+4], j, k, nb, newlen, r; | |
321 double f_min, f_max, alfa, beta, u, *x = work, y; | |
322 /* substitute and remove fixed variables */ | |
323 newlen = 0; | |
324 for (k = 1; k <= len; k++) | |
325 { j = ind[k]; | |
326 if (lpx_get_col_type(lp, j) == LPX_FX) | |
327 val[0] -= val[k] * lpx_get_col_lb(lp, j); | |
328 else | |
329 { newlen++; | |
330 ind[newlen] = ind[k]; | |
331 val[newlen] = val[k]; | |
332 } | |
333 } | |
334 len = newlen; | |
335 /* move binary variables to the beginning of the list so that | |
336 elements 1, 2, ..., nb correspond to binary variables, and | |
337 elements nb+1, nb+2, ..., len correspond to rest variables */ | |
338 nb = 0; | |
339 for (k = 1; k <= len; k++) | |
340 { j = ind[k]; | |
341 if (lpx_get_col_kind(lp, j) == LPX_IV && | |
342 lpx_get_col_type(lp, j) == LPX_DB && | |
343 lpx_get_col_lb(lp, j) == 0.0 && | |
344 lpx_get_col_ub(lp, j) == 1.0) | |
345 { /* binary variable */ | |
346 int ind_k; | |
347 double val_k; | |
348 nb++; | |
349 ind_k = ind[nb], val_k = val[nb]; | |
350 ind[nb] = ind[k], val[nb] = val[k]; | |
351 ind[k] = ind_k, val[k] = val_k; | |
352 } | |
353 } | |
354 /* now the specified row has the form: | |
355 sum a[j]*x[j] + sum a[j]*y[j] <= b, | |
356 where x[j] are binary variables, y[j] are rest variables */ | |
357 /* at least two binary variables are needed */ | |
358 if (nb < 2) return 0; | |
359 /* compute implied lower and upper bounds for sum a[j]*y[j] */ | |
360 f_min = f_max = 0.0; | |
361 for (k = nb+1; k <= len; k++) | |
362 { j = ind[k]; | |
363 /* both bounds must be finite */ | |
364 if (lpx_get_col_type(lp, j) != LPX_DB) return 0; | |
365 if (val[k] > 0.0) | |
366 { f_min += val[k] * lpx_get_col_lb(lp, j); | |
367 f_max += val[k] * lpx_get_col_ub(lp, j); | |
368 } | |
369 else | |
370 { f_min += val[k] * lpx_get_col_ub(lp, j); | |
371 f_max += val[k] * lpx_get_col_lb(lp, j); | |
372 } | |
373 } | |
374 /* sum a[j]*x[j] + sum a[j]*y[j] <= b ===> | |
375 sum a[j]*x[j] + (sum a[j]*y[j] - f_min) <= b - f_min ===> | |
376 sum a[j]*x[j] + y <= b - f_min, | |
377 where y = sum a[j]*y[j] - f_min; | |
378 note that 0 <= y <= u, u = f_max - f_min */ | |
379 /* determine upper bound of y */ | |
380 u = f_max - f_min; | |
381 /* determine value of y at the current point */ | |
382 y = 0.0; | |
383 for (k = nb+1; k <= len; k++) | |
384 { j = ind[k]; | |
385 y += val[k] * lpx_get_col_prim(lp, j); | |
386 } | |
387 y -= f_min; | |
388 if (y < 0.0) y = 0.0; | |
389 if (y > u) y = u; | |
390 /* modify the right-hand side b */ | |
391 val[0] -= f_min; | |
392 /* now the transformed row has the form: | |
393 sum a[j]*x[j] + y <= b, where 0 <= y <= u */ | |
394 /* determine values of x[j] at the current point */ | |
395 for (k = 1; k <= nb; k++) | |
396 { j = ind[k]; | |
397 x[k] = lpx_get_col_prim(lp, j); | |
398 if (x[k] < 0.0) x[k] = 0.0; | |
399 if (x[k] > 1.0) x[k] = 1.0; | |
400 } | |
401 /* if a[j] < 0, replace x[j] by its complement 1 - x'[j] */ | |
402 for (k = 1; k <= nb; k++) | |
403 { if (val[k] < 0.0) | |
404 { ind[k] = - ind[k]; | |
405 val[k] = - val[k]; | |
406 val[0] += val[k]; | |
407 x[k] = 1.0 - x[k]; | |
408 } | |
409 } | |
410 /* try to generate a mixed cover cut for the transformed row */ | |
411 r = cover(nb, val, val[0], u, x, y, cov, &alfa, &beta); | |
412 if (r == 0) return 0; | |
413 xassert(2 <= r && r <= 4); | |
414 /* now the cut is in the form: | |
415 sum{j in C} x[j] + alfa * y <= beta */ | |
416 /* store the right-hand side beta */ | |
417 ind[0] = 0, val[0] = beta; | |
418 /* restore the original ordinal numbers of x[j] */ | |
419 for (j = 1; j <= r; j++) cov[j] = ind[cov[j]]; | |
420 /* store cut coefficients at binary variables complementing back | |
421 the variables having negative row coefficients */ | |
422 xassert(r <= nb); | |
423 for (k = 1; k <= r; k++) | |
424 { if (cov[k] > 0) | |
425 { ind[k] = +cov[k]; | |
426 val[k] = +1.0; | |
427 } | |
428 else | |
429 { ind[k] = -cov[k]; | |
430 val[k] = -1.0; | |
431 val[0] -= 1.0; | |
432 } | |
433 } | |
434 /* substitute y = sum a[j]*y[j] - f_min */ | |
435 for (k = nb+1; k <= len; k++) | |
436 { r++; | |
437 ind[r] = ind[k]; | |
438 val[r] = alfa * val[k]; | |
439 } | |
440 val[0] += alfa * f_min; | |
441 xassert(r <= len); | |
442 len = r; | |
443 return len; | |
444 } | |
445 | |
446 /*---------------------------------------------------------------------- | |
447 -- lpx_eval_row - compute explictily specified row. | |
448 -- | |
449 -- SYNOPSIS | |
450 -- | |
451 -- #include "glplpx.h" | |
452 -- double lpx_eval_row(LPX *lp, int len, int ind[], double val[]); | |
453 -- | |
454 -- DESCRIPTION | |
455 -- | |
456 -- The routine lpx_eval_row computes the primal value of an explicitly | |
457 -- specified row using current values of structural variables. | |
458 -- | |
459 -- The explicitly specified row may be thought as a linear form: | |
460 -- | |
461 -- y = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], | |
462 -- | |
463 -- where y is an auxiliary variable for this row, a[j] are coefficients | |
464 -- of the linear form, x[m+j] are structural variables. | |
465 -- | |
466 -- On entry column indices and numerical values of non-zero elements of | |
467 -- the row should be stored in locations ind[1], ..., ind[len] and | |
468 -- val[1], ..., val[len], where len is the number of non-zero elements. | |
469 -- The array ind and val are not changed on exit. | |
470 -- | |
471 -- RETURNS | |
472 -- | |
473 -- The routine returns a computed value of y, the auxiliary variable of | |
474 -- the specified row. */ | |
475 | |
476 static double lpx_eval_row(LPX *lp, int len, int ind[], double val[]) | |
477 { int n = lpx_get_num_cols(lp); | |
478 int j, k; | |
479 double sum = 0.0; | |
480 if (len < 0) | |
481 xerror("lpx_eval_row: len = %d; invalid row length\n", len); | |
482 for (k = 1; k <= len; k++) | |
483 { j = ind[k]; | |
484 if (!(1 <= j && j <= n)) | |
485 xerror("lpx_eval_row: j = %d; column number out of range\n", | |
486 j); | |
487 sum += val[k] * lpx_get_col_prim(lp, j); | |
488 } | |
489 return sum; | |
490 } | |
491 | |
492 /*********************************************************************** | |
493 * NAME | |
494 * | |
495 * ios_cov_gen - generate mixed cover cuts | |
496 * | |
497 * SYNOPSIS | |
498 * | |
499 * #include "glpios.h" | |
500 * void ios_cov_gen(glp_tree *tree); | |
501 * | |
502 * DESCRIPTION | |
503 * | |
504 * The routine ios_cov_gen generates mixed cover cuts for the current | |
505 * point and adds them to the cut pool. */ | |
506 | |
507 void ios_cov_gen(glp_tree *tree) | |
508 { glp_prob *prob = tree->mip; | |
509 int m = lpx_get_num_rows(prob); | |
510 int n = lpx_get_num_cols(prob); | |
511 int i, k, type, kase, len, *ind; | |
512 double r, *val, *work; | |
513 xassert(lpx_get_status(prob) == LPX_OPT); | |
514 /* allocate working arrays */ | |
515 ind = xcalloc(1+n, sizeof(int)); | |
516 val = xcalloc(1+n, sizeof(double)); | |
517 work = xcalloc(1+n, sizeof(double)); | |
518 /* look through all rows */ | |
519 for (i = 1; i <= m; i++) | |
520 for (kase = 1; kase <= 2; kase++) | |
521 { type = lpx_get_row_type(prob, i); | |
522 if (kase == 1) | |
523 { /* consider rows of '<=' type */ | |
524 if (!(type == LPX_UP || type == LPX_DB)) continue; | |
525 len = lpx_get_mat_row(prob, i, ind, val); | |
526 val[0] = lpx_get_row_ub(prob, i); | |
527 } | |
528 else | |
529 { /* consider rows of '>=' type */ | |
530 if (!(type == LPX_LO || type == LPX_DB)) continue; | |
531 len = lpx_get_mat_row(prob, i, ind, val); | |
532 for (k = 1; k <= len; k++) val[k] = - val[k]; | |
533 val[0] = - lpx_get_row_lb(prob, i); | |
534 } | |
535 /* generate mixed cover cut: | |
536 sum{j in J} a[j] * x[j] <= b */ | |
537 len = lpx_cover_cut(prob, len, ind, val, work); | |
538 if (len == 0) continue; | |
539 /* at the current point the cut inequality is violated, i.e. | |
540 sum{j in J} a[j] * x[j] - b > 0 */ | |
541 r = lpx_eval_row(prob, len, ind, val) - val[0]; | |
542 if (r < 1e-3) continue; | |
543 /* add the cut to the cut pool */ | |
544 glp_ios_add_row(tree, NULL, GLP_RF_COV, 0, len, ind, val, | |
545 GLP_UP, val[0]); | |
546 } | |
547 /* free working arrays */ | |
548 xfree(ind); | |
549 xfree(val); | |
550 xfree(work); | |
551 return; | |
552 } | |
553 | |
554 /* eof */ |