glpios07.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* glpios07.c (mixed cover cut generator) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpios.h"

/*----------------------------------------------------------------------
-- COVER INEQUALITIES
--
-- Consider the set of feasible solutions to 0-1 knapsack problem:
--
--    sum a[j]*x[j] <= b,                                            (1)
--  j in J
--
--    x[j] is binary,                                                (2)
--
-- where, wlog, we assume that a[j] > 0 (since 0-1 variables can be
-- complemented) and a[j] <= b (since a[j] > b implies x[j] = 0).
--
-- A set C within J is called a cover if
--
--    sum a[j] > b.                                                  (3)
--  j in C
--
-- For any cover C the inequality
--
--    sum x[j] <= |C| - 1                                            (4)
--  j in C
--
-- is called a cover inequality and is valid for (1)-(2).
--
-- MIXED COVER INEQUALITIES
--
-- Consider the set of feasible solutions to mixed knapsack problem:
--
--    sum a[j]*x[j] + y <= b,                                        (5)
--  j in J
--
--    x[j] is binary,                                                (6)
--
--    0 <= y <= u is continuous,                                     (7)
--
-- where again we assume that a[j] > 0.
--
-- Let C within J be some set. From (1)-(4) it follows that
--
--    sum a[j] > b - y                                               (8)
--  j in C
--
-- implies
--
--    sum x[j] <= |C| - 1.                                           (9)
--  j in C
--
-- Thus, we need to modify the inequality (9) in such a way that it be
-- a constraint only if the condition (8) is satisfied.
--
-- Consider the following inequality:
--
--    sum x[j] <= |C| - t.                                          (10)
--  j in C
--
-- If 0 < t <= 1, then (10) is equivalent to (9), because all x[j] are
-- binary variables. On the other hand, if t <= 0, (10) being satisfied
-- for any values of x[j] is not a constraint.
--
-- Let
--
--    t' = sum a[j] + y - b.                                        (11)
--       j in C
--
-- It is understood that the condition t' > 0 is equivalent to (8).
-- Besides, from (6)-(7) it follows that t' has an implied upper bound:
--
--    t'max = sum a[j] + u - b.                                     (12)
--          j in C
--
-- This allows to express the parameter t having desired properties:
--
--    t = t' / t'max.                                               (13)
--
-- In fact, t <= 1 by definition, and t > 0 being equivalent to t' > 0
-- is equivalent to (8).
--
-- Thus, the inequality (10), where t is given by formula (13) is valid
-- for (5)-(7).
--
-- Note that if u = 0, then y = 0, so t = 1, and the conditions (8) and
-- (10) is transformed to the conditions (3) and (4).
--
-- GENERATING MIXED COVER CUTS
--
-- To generate a mixed cover cut in the form (10) we need to find such
-- set C which satisfies to the inequality (8) and for which, in turn,
-- the inequality (10) is violated in the current point.
--
-- Substituting t from (13) to (10) gives:
--
--                        1
--    sum x[j] <= |C| - -----  (sum a[j] + y - b),                  (14)
--  j in C              t'max j in C
--
-- and finally we have the cut inequality in the standard form:
--
--    sum x[j] + alfa * y <= beta,                                  (15)
--  j in C
--
-- where:
--
--    alfa = 1 / t'max,                                             (16)
--
--    beta = |C| - alfa *  (sum a[j] - b).                          (17)
--                        j in C                                      */

#if 1
#define MAXTRY 1000
#else
#define MAXTRY 10000
#endif

static int cover2(int n, double a[], double b, double u, double x[],
      double y, int cov[], double *_alfa, double *_beta)
{     /* try to generate mixed cover cut using two-element cover */
      int i, j, try = 0, ret = 0;
      double eps, alfa, beta, temp, rmax = 0.001;
      eps = 0.001 * (1.0 + fabs(b));
      for (i = 0+1; i <= n; i++)
      for (j = i+1; j <= n; j++)
      {  /* C = {i, j} */
         try++;
         if (try > MAXTRY) goto done;
         /* check if condition (8) is satisfied */
         if (a[i] + a[j] + y > b + eps)
         {  /* compute parameters for inequality (15) */
            temp = a[i] + a[j] - b;
            alfa = 1.0 / (temp + u);
            beta = 2.0 - alfa * temp;
            /* compute violation of inequality (15) */
            temp = x[i] + x[j] + alfa * y - beta;
            /* choose C providing maximum violation */
            if (rmax < temp)
            {  rmax = temp;
               cov[1] = i;
               cov[2] = j;
               *_alfa = alfa;
               *_beta = beta;
               ret = 1;
            }
         }
      }
done: return ret;
}

static int cover3(int n, double a[], double b, double u, double x[],
      double y, int cov[], double *_alfa, double *_beta)
{     /* try to generate mixed cover cut using three-element cover */
      int i, j, k, try = 0, ret = 0;
      double eps, alfa, beta, temp, rmax = 0.001;
      eps = 0.001 * (1.0 + fabs(b));
      for (i = 0+1; i <= n; i++)
      for (j = i+1; j <= n; j++)
      for (k = j+1; k <= n; k++)
      {  /* C = {i, j, k} */
         try++;
         if (try > MAXTRY) goto done;
         /* check if condition (8) is satisfied */
         if (a[i] + a[j] + a[k] + y > b + eps)
         {  /* compute parameters for inequality (15) */
            temp = a[i] + a[j] + a[k] - b;
            alfa = 1.0 / (temp + u);
            beta = 3.0 - alfa * temp;
            /* compute violation of inequality (15) */
            temp = x[i] + x[j] + x[k] + alfa * y - beta;
            /* choose C providing maximum violation */
            if (rmax < temp)
            {  rmax = temp;
               cov[1] = i;
               cov[2] = j;
               cov[3] = k;
               *_alfa = alfa;
               *_beta = beta;
               ret = 1;
            }
         }
      }
done: return ret;
}

static int cover4(int n, double a[], double b, double u, double x[],
      double y, int cov[], double *_alfa, double *_beta)
{     /* try to generate mixed cover cut using four-element cover */
      int i, j, k, l, try = 0, ret = 0;
      double eps, alfa, beta, temp, rmax = 0.001;
      eps = 0.001 * (1.0 + fabs(b));
      for (i = 0+1; i <= n; i++)
      for (j = i+1; j <= n; j++)
      for (k = j+1; k <= n; k++)
      for (l = k+1; l <= n; l++)
      {  /* C = {i, j, k, l} */
         try++;
         if (try > MAXTRY) goto done;
         /* check if condition (8) is satisfied */
         if (a[i] + a[j] + a[k] + a[l] + y > b + eps)
         {  /* compute parameters for inequality (15) */
            temp = a[i] + a[j] + a[k] + a[l] - b;
            alfa = 1.0 / (temp + u);
            beta = 4.0 - alfa * temp;
            /* compute violation of inequality (15) */
            temp = x[i] + x[j] + x[k] + x[l] + alfa * y - beta;
            /* choose C providing maximum violation */
            if (rmax < temp)
            {  rmax = temp;
               cov[1] = i;
               cov[2] = j;
               cov[3] = k;
               cov[4] = l;
               *_alfa = alfa;
               *_beta = beta;
               ret = 1;
            }
         }
      }
done: return ret;
}

static int cover(int n, double a[], double b, double u, double x[],
      double y, int cov[], double *alfa, double *beta)
{     /* try to generate mixed cover cut;
         input (see (5)):
         n        is the number of binary variables;
         a[1:n]   are coefficients at binary variables;
         b        is the right-hand side;
         u        is upper bound of continuous variable;
         x[1:n]   are values of binary variables at current point;
         y        is value of continuous variable at current point;
         output (see (15), (16), (17)):
         cov[1:r] are indices of binary variables included in cover C,
                  where r is the set cardinality returned on exit;
         alfa     coefficient at continuous variable;
         beta     is the right-hand side; */
      int j;
      /* perform some sanity checks */
      xassert(n >= 2);
      for (j = 1; j <= n; j++) xassert(a[j] > 0.0);
#if 1 /* ??? */
      xassert(b > -1e-5);
#else
      xassert(b > 0.0);
#endif
      xassert(u >= 0.0);
      for (j = 1; j <= n; j++) xassert(0.0 <= x[j] && x[j] <= 1.0);
      xassert(0.0 <= y && y <= u);
      /* try to generate mixed cover cut */
      if (cover2(n, a, b, u, x, y, cov, alfa, beta)) return 2;
      if (cover3(n, a, b, u, x, y, cov, alfa, beta)) return 3;
      if (cover4(n, a, b, u, x, y, cov, alfa, beta)) return 4;
      return 0;
}