lemon-project-template-glpk
view deps/glpk/src/glpios07.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glpios07.c (mixed cover cut generator) */
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 ***********************************************************************/
25 #include "glpios.h"
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 */
135 #if 1
136 #define MAXTRY 1000
137 #else
138 #define MAXTRY 10000
139 #endif
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 }
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 }
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 }
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 }
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. */
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 }
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. */
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 }
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. */
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 }
554 /* eof */