glpios03.c

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

-- This routine also selects the branch to be solved next where integer
-- infeasibility of the chosen variable is less than in other one. */

static void branch_first(glp_tree *tree)
{     glp_prob *lp = tree->mip;
      int n = lp->n;
      int j, next;
      double beta;
      /* choose the column to branch on */
      for (j = 1; j <= n; j++)
         if (tree->non_int[j]) break;
      xassert(1 <= j && j <= n);
      /* select the branch to be solved next */
      beta = lpx_get_col_prim(lp, j);
      if (beta - floor(beta) < ceil(beta) - beta)
         next = GLP_DN_BRNCH;
      else
         next = GLP_UP_BRNCH;
      /* perform branching */
      glp_ios_branch_upon(tree, j, next);
      return;
}

/*----------------------------------------------------------------------
-- branch_last - choose last branching variable.
--
-- This routine looks up the list of structural variables and chooses
-- the last one, which is of integer kind and has fractional value in
-- optimal solution of the current LP relaxation.
--
-- This routine also selects the branch to be solved next where integer
-- infeasibility of the chosen variable is less than in other one. */

static void branch_last(glp_tree *tree)
{     glp_prob *lp = tree->mip;
      int n = lp->n;
      int j, next;
      double beta;
      /* choose the column to branch on */
      for (j = n; j >= 1; j--)
         if (tree->non_int[j]) break;
      xassert(1 <= j && j <= n);
      /* select the branch to be solved next */
      beta = lpx_get_col_prim(lp, j);
      if (beta - floor(beta) < ceil(beta) - beta)
         next = GLP_DN_BRNCH;
      else
         next = GLP_UP_BRNCH;
      /* perform branching */
      glp_ios_branch_upon(tree, j, next);
      return;
}

/*----------------------------------------------------------------------
-- branch_drtom - choose branching variable with Driebeck-Tomlin heur.
--
-- This routine chooses a structural variable, which is required to be
-- integral and has fractional value in optimal solution of the current
-- LP relaxation, using a heuristic proposed by Driebeck and Tomlin.
--
-- The routine also selects the branch to be solved next, again due to
-- Driebeck and Tomlin.
--
-- This routine is based on the heuristic proposed in:
--
-- Driebeck N.J. An algorithm for the solution of mixed-integer
-- programming problems, Management Science, 12: 576-87 (1966)
--
-- and improved in:
--
-- Tomlin J.A. Branch and bound methods for integer and non-convex
-- programming, in J.Abadie (ed.), Integer and Nonlinear Programming,
-- North-Holland, Amsterdam, pp. 437-50 (1970).
--
-- Must note that this heuristic is time-expensive, because computing
-- one-step degradation (see the routine below) requires one BTRAN for
-- each fractional-valued structural variable. */

#define int_col(j) (lpx_get_col_kind(tree->mip, j) == LPX_IV)

#if 0
struct col { int j; double f; };

static int fcmp(const void *x1, const void *x2)
{     const struct col *c1 = x1, *c2 = x2;
      if (c1->f > c2->f) return -1;
      if (c1->f < c2->f) return +1;
      return 0;
}
#endif

static void branch_drtom(glp_tree *tree)
{     glp_prob *lp = tree->mip;
      int m = lp->m;
      int n = lp->n;
      int *non_int = tree->non_int;
      int j, jj, k, t, next, kase, len, stat, *ind;
      double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,
         dd_dn, dd_up, degrad, *val;
#if 0
      struct col *list;
      int size;
#endif
      /* basic solution of LP relaxation must be optimal */
      xassert(lpx_get_status(lp) == LPX_OPT);
      /* allocate working arrays */
      ind = xcalloc(1+n, sizeof(int));
      val = xcalloc(1+n, sizeof(double));
#if 0
      list = xcalloc(1+n, sizeof(struct col));
      /* build the list of fractional-valued columns */
      size = 0;
      for (j = 1; j <= n; j++)
      {  if (non_int[j])
         {  double t1, t2;
            size++;
            list[size].j = j;
            x = lpx_get_col_prim(lp, j);
            t1 = x - floor(x);
            t2 = ceil(x) - x;
            list[size].f = (t1 <= t2 ? t1 : t2);
         }
      }
      /* sort the column list in descending fractionality */
      qsort(&list[1], size, sizeof(struct col), fcmp);
      /* limit the number of columns to be considered */
      if (size > 50) size = 50;
#endif
      /* nothing has been chosen so far */
      jj = 0, degrad = -1.0;
      /* walk through the list of columns (structural variables) */
#if 1
      for (j = 1; j <= n; j++)
      {  /* if j-th column is not marked as fractional, skip it */
         if (!non_int[j]) continue;
#else
      while (size > 0)
      {  j = list[size--].j;
#endif
         /* obtain (fractional) value of j-th column in basic solution
            of LP relaxation */
         x = lpx_get_col_prim(lp, j);
         /* since the value of j-th column is fractional, the column is
            basic; compute corresponding row of the simplex table */
         len = lpx_eval_tab_row(lp, m+j, ind, val);
         /* the following fragment computes a change in the objective
            function: delta Z = new Z - old Z, where old Z is the
            objective value in the current optimal basis, and new Z is
            the objective value in the adjacent basis, for two cases:
            1) if new upper bound ub' = floor(x[j]) is introduced for
               j-th column (down branch);
            2) if new lower bound lb' = ceil(x[j]) is introduced for
               j-th column (up branch);
            since in both cases the solution remaining dual feasible
            becomes primal infeasible, one implicit simplex iteration
            is performed to determine the change delta Z;
            it is obvious that new Z, which is never better than old Z,
            is a lower (minimization) or upper (maximization) bound of
            the objective function for down- and up-branches. */
         for (kase = -1; kase <= +1; kase += 2)
         {  /* if kase < 0, the new upper bound of x[j] is introduced;
               in this case x[j] should decrease in order to leave the
               basis and go to its new upper bound */
            /* if kase > 0, the new lower bound of x[j] is introduced;
               in this case x[j] should increase in order to leave the
               basis and go to its new lower bound */
            /* apply the dual ratio test in order to determine which
               auxiliary or structural variable should enter the basis
               to keep dual feasibility */
            k = lpx_dual_ratio_test(lp, len, ind, val, kase, 1e-8);
            /* if no non-basic variable has been chosen, LP relaxation
               of corresponding branch being primal infeasible and dual
               unbounded has no primal feasible solution; in this case
               the change delta Z is formally set to infinity */
            if (k == 0)
            {  delta_z =
                  (tree->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX);
               goto skip;
            }
            /* row of the simplex table that corresponds to non-basic
               variable x[k] choosen by the dual ratio test is:
                  x[j] = ... + alfa * x[k] + ...
               where alfa is the influence coefficient (an element of
               the simplex table row) */
            /* determine the coefficient alfa */
            for (t = 1; t <= len; t++) if (ind[t] == k) break;
            xassert(1 <= t && t <= len);
            alfa = val[t];
            /* since in the adjacent basis the variable x[j] becomes
               non-basic, knowing its value in the current basis we can
               determine its change delta x[j] = new x[j] - old x[j] */
            delta_j = (kase < 0 ? floor(x) : ceil(x)) - x;
            /* and knowing the coefficient alfa we can determine the
               corresponding change delta x[k] = new x[k] - old x[k],
               where old x[k] is a value of x[k] in the current basis,
               and new x[k] is a value of x[k] in the adjacent basis */
            delta_k = delta_j / alfa;
            /* Tomlin noticed that if the variable x[k] is of integer
               kind, its change cannot be less (eventually) than one in
               the magnitude */
            if (k > m && int_col(k-m))
            {  /* x[k] is structural integer variable */
               if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3)
               {  if (delta_k > 0.0)
                     delta_k = ceil(delta_k);  /* +3.14 -> +4 */
                  else
                     delta_k = floor(delta_k); /* -3.14 -> -4 */
               }
            }
            /* now determine the status and reduced cost of x[k] in the
               current basis */
            if (k <= m)
            {  stat = lpx_get_row_stat(lp, k);
               dk = lpx_get_row_dual(lp, k);
            }
            else
            {  stat = lpx_get_col_stat(lp, k-m);
               dk = lpx_get_col_dual(lp, k-m);
            }
            /* if the current basis is dual degenerative, some reduced
               costs which are close to zero may have wrong sign due to
               round-off errors, so correct the sign of d[k] */
            switch (tree->mip->dir)
            {  case GLP_MIN:
                  if (stat == LPX_NL && dk < 0.0 ||
                      stat == LPX_NU && dk > 0.0 ||
                      stat == LPX_NF) dk = 0.0;
                  break;
               case GLP_MAX:
                  if (stat == LPX_NL && dk > 0.0 ||
                      stat == LPX_NU && dk < 0.0 ||
                      stat == LPX_NF) dk = 0.0;
                  break;
               default:
                  xassert(tree != tree);
            }
            /* now knowing the change of x[k] and its reduced cost d[k]
               we can compute the corresponding change in the objective
               function delta Z = new Z - old Z = d[k] * delta x[k];
               note that due to Tomlin's modification new Z can be even
               worse than in the adjacent basis */
            delta_z = dk * delta_k;
skip:       /* new Z is never better than old Z, therefore the change
               delta Z is always non-negative (in case of minimization)
               or non-positive (in case of maximization) */
            switch (tree->mip->dir)
            {  case GLP_MIN: xassert(delta_z >= 0.0); break;
               case GLP_MAX: xassert(delta_z <= 0.0); break;
               default: xassert(tree != tree);
            }
            /* save the change in the objective fnction for down- and
               up-branches, respectively */
            if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;
         }
         /* thus, in down-branch no integer feasible solution can be
            better than Z + dz_dn, and in up-branch no integer feasible
            solution can be better than Z + dz_up, where Z is value of
            the objective function in the current basis */
         /* following the heuristic by Driebeck and Tomlin we choose a
            column (i.e. structural variable) which provides largest
            degradation of the objective function in some of branches;
            besides, we select the branch with smaller degradation to
            be solved next and keep other branch with larger degradation
            in the active list hoping to minimize the number of further
            backtrackings */
         if (degrad < fabs(dz_dn) || degrad < fabs(dz_up))
         {  jj = j;
            if (fabs(dz_dn) < fabs(dz_up))
            {  /* select down branch to be solved next */
               next = GLP_DN_BRNCH;
               degrad = fabs(dz_up);
            }
            else
            {  /* select up branch to be solved next */
               next = GLP_UP_BRNCH;
               degrad = fabs(dz_dn);
            }
            /* save the objective changes for printing */
            dd_dn = dz_dn, dd_up = dz_up;
            /* if down- or up-branch has no feasible solution, we does
               not need to consider other candidates (in principle, the
               corresponding branch could be pruned right now) */
            if (degrad == DBL_MAX) break;
         }
      }
      /* free working arrays */
      xfree(ind);
      xfree(val);
#if 0
      xfree(list);
#endif
      /* something must be chosen */
      xassert(1 <= jj && jj <= n);
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
      {  xprintf("branch_drtom: column %d chosen to branch on\n", jj);
         if (fabs(dd_dn) == DBL_MAX)
            xprintf("branch_drtom: down-branch is infeasible\n");
         else
            xprintf("branch_drtom: down-branch bound is %.9e\n",
               lpx_get_obj_val(lp) + dd_dn);
         if (fabs(dd_up) == DBL_MAX)
            xprintf("branch_drtom: up-branch   is infeasible\n");
         else
            xprintf("branch_drtom: up-branch   bound is %.9e\n",
               lpx_get_obj_val(lp) + dd_up);
      }
      /* perform branching */
      glp_ios_branch_upon(tree, jj, next);
      return;
}

/*----------------------------------------------------------------------
-- branch_mostf - choose most fractional branching variable.
--
-- This routine looks up the list of structural variables and chooses
-- that one, which is of integer kind and has most fractional value in
-- optimal solution of the current LP relaxation.
--
-- This routine also selects the branch to be solved next where integer
-- infeasibility of the chosen variable is less than in other one.
*
*  Martin notices that "...most infeasible is as good as random...". */

static void branch_mostf(glp_tree *tree)
{     glp_prob *lp = tree->mip;
      int n = lp->n;
      int j, jj, next;
      double beta, most, temp;
      /* choose the column to branch on */
      jj = 0, most = DBL_MAX;
      for (j = 1; j <= n; j++)
      {  if (tree->non_int[j])
         {  beta = lpx_get_col_prim(lp, j);
            temp = floor(beta) + 0.5;
            if (most > fabs(beta - temp))
            {  jj = j, most = fabs(beta - temp);
               if (beta < temp)
                  next = GLP_DN_BRNCH;
               else
                  next = GLP_UP_BRNCH;
            }
         }
      }
      /* perform branching */
      glp_ios_branch_upon(tree, jj, next);
      return;
}

/***********************************************************************
*  branch_on - perform branching on specified variable
*