COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios07.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 18.6 KB
Line 
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 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
141static 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      }
171done: return ret;
172}
173
174static 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      }
206done: return ret;
207}
208
209static 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      }
243done: return ret;
244}
245
246static 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
318static 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
476static 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
507void 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 */
Note: See TracBrowser for help on using the repository browser.