glpios03.c

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

*  This routine performs branching on j-th column (structural variable)
*  of the current subproblem. The specified column must be of integer
*  kind and must have a fractional value in optimal basic solution of
*  LP relaxation of the current subproblem (i.e. only columns for which
*  the flag non_int[j] is set are valid candidates to branch on).
*
*  Let x be j-th structural variable, and beta be its primal fractional
*  value in the current basic solution. Branching on j-th variable is
*  dividing the current subproblem into two new subproblems, which are
*  identical to the current subproblem with the following exception: in
*  the first subproblem that begins the down-branch x has a new upper
*  bound x <= floor(beta), and in the second subproblem that begins the
*  up-branch x has a new lower bound x >= ceil(beta).
*
*  Depending on estimation of the local bound for down- and up-branches
*  this routine returns one of the following:
*
*  0 - both branches have been created;
*  1 - one branch has no feasible solution and has been pruned; other
*      branch has became the current subproblem;
*  2 - both branches have no feasible solution and have been pruned;
*      subproblem selection is needed. */

static int branch_on(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int n = mip->n;
      int j = tree->br_var;
      int next = tree->br_sel;
      int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
      double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
      /* determine bounds and value of x[j] in optimal solution to LP
         relaxation of the current subproblem */
      xassert(1 <= j && j <= n);
      type = mip->col[j]->type;
      lb = mip->col[j]->lb;
      ub = mip->col[j]->ub;
      beta = mip->col[j]->prim;
      /* determine new bounds of x[j] for down- and up-branches */
      new_ub = floor(beta);
      new_lb = ceil(beta);
      switch (type)
      {  case GLP_FR:
            dn_type = GLP_UP;
            up_type = GLP_LO;
            break;
         case GLP_LO:
            xassert(lb <= new_ub);
            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
            xassert(lb + 1.0 <= new_lb);
            up_type = GLP_LO;
            break;
         case GLP_UP:
            xassert(new_ub <= ub - 1.0);
            dn_type = GLP_UP;
            xassert(new_lb <= ub);
            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
            break;
         case GLP_DB:
            xassert(lb <= new_ub && new_ub <= ub - 1.0);
            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
            xassert(lb + 1.0 <= new_lb && new_lb <= ub);
            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
            break;
         default:
            xassert(type != type);
      }
      /* compute local bounds to LP relaxation for both branches */
      ios_eval_degrad(tree, j, &dn_lp, &up_lp);
      /* and improve them by rounding */
      dn_bnd = ios_round_bound(tree, dn_lp);
      up_bnd = ios_round_bound(tree, up_lp);
      /* check local bounds for down- and up-branches */
      dn_bad = !ios_is_hopeful(tree, dn_bnd);
      up_bad = !ios_is_hopeful(tree, up_bnd);
      if (dn_bad && up_bad)
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Both down- and up-branches are hopeless\n");
         ret = 2;
         goto done;
      }
      else if (up_bad)
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Up-branch is hopeless\n");
         glp_set_col_bnds(mip, j, dn_type, lb, new_ub);
         tree->curr->lp_obj = dn_lp;
         if (mip->dir == GLP_MIN)
         {  if (tree->curr->bound < dn_bnd)
                tree->curr->bound = dn_bnd;
         }
         else if (mip->dir == GLP_MAX)
         {  if (tree->curr->bound > dn_bnd)
                tree->curr->bound = dn_bnd;
         }
         else
            xassert(mip != mip);
         ret = 1;
         goto done;
      }
      else if (dn_bad)
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Down-branch is hopeless\n");
         glp_set_col_bnds(mip, j, up_type, new_lb, ub);
         tree->curr->lp_obj = up_lp;
         if (mip->dir == GLP_MIN)
         {  if (tree->curr->bound < up_bnd)
                tree->curr->bound = up_bnd;
         }
         else if (mip->dir == GLP_MAX)
         {  if (tree->curr->bound > up_bnd)
                tree->curr->bound = up_bnd;
         }
         else
            xassert(mip != mip);
         ret = 1;
         goto done;
      }
      /* both down- and up-branches seem to be hopeful */
      /* branching */
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
         xprintf("Branching on column %d, primal value is %.9e\n",
            j, beta);
      /* determine the reference number of the current subproblem */
      xassert(tree->curr != NULL);
      p = tree->curr->p;
      tree->curr->br_var = j;
      tree->curr->br_val = beta;
      /* freeze the current subproblem */
      ios_freeze_node(tree);
      /* create two clones of the current subproblem; the first clone
         begins the down-branch, the second one begins the up-branch */
      ios_clone_node(tree, p, 2, clone);
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
         xprintf("Node %d begins down branch, node %d begins up branch "
            "\n", clone[1], clone[2]);
      /* set new upper bound of j-th column in the down-branch */
      ios_revive_node(tree, clone[1]);
      glp_set_col_bnds(mip, j, dn_type, lb, new_ub);
      tree->curr->lp_obj = dn_lp;
      if (mip->dir == GLP_MIN)
      {  if (tree->curr->bound < dn_bnd)
             tree->curr->bound = dn_bnd;
      }
      else if (mip->dir == GLP_MAX)
      {  if (tree->curr->bound > dn_bnd)
             tree->curr->bound = dn_bnd;
      }
      else
         xassert(mip != mip);
      ios_freeze_node(tree);
      /* set new lower bound of j-th column in the up-branch */
      ios_revive_node(tree, clone[2]);
      glp_set_col_bnds(mip, j, up_type, new_lb, ub);
      tree->curr->lp_obj = up_lp;
      if (mip->dir == GLP_MIN)
      {  if (tree->curr->bound < up_bnd)
             tree->curr->bound = up_bnd;
      }
      else if (mip->dir == GLP_MAX)
      {  if (tree->curr->bound > up_bnd)
             tree->curr->bound = up_bnd;
      }
      else
         xassert(mip != mip);
      ios_freeze_node(tree);
      /* revive the subproblem to be solved next */
      if (next == GLP_NO_BRNCH)
         ;
      else if (next == GLP_DN_BRNCH)
         ios_revive_node(tree, clone[1]);
      else if (next == GLP_UP_BRNCH)
         ios_revive_node(tree, clone[2]);
      else
         xassert(next != next);
      ret = 0;
done: return ret;
}

/*----------------------------------------------------------------------
-- cleanup_the_tree - prune hopeless branches from the tree.
--
-- This routine walks through the active list and checks the local
-- bound for every active subproblem. If the local bound indicates that
-- the subproblem cannot have integer optimal solution better than the
-- incumbent objective value, the routine deletes such subproblem that,
-- in turn, involves pruning the corresponding branch of the tree. */

static void cleanup_the_tree(glp_tree *tree)
{     IOSNPD *node, *next_node;
      int count = 0;
      /* the global bound must exist */
      xassert(tree->mip->mip_stat == GLP_FEAS);
      /* walk through the list of active subproblems */
      for (node = tree->head; node != NULL; node = next_node)
      {  /* deleting some active problem node may involve deleting its
            parents recursively; however, all its parents being created
            *before* it are always *precede* it in the node list, so
            the next problem node is never affected by such deletion */
         next_node = node->next;
         /* if the branch is hopeless, prune it */
         if (!is_branch_hopeful(tree, node->p))
            ios_delete_node(tree, node->p), count++;
      }
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
      {  if (count == 1)
            xprintf("One hopeless branch has been pruned\n");
         else if (count > 1)
            xprintf("%d hopeless branches have been pruned\n", count);
      }
      return;
}

/*----------------------------------------------------------------------
-- btrack_most_feas - select "most integer feasible" subproblem.
--
-- This routine selects from the active list a subproblem to be solved
-- next, whose parent has minimal sum of integer infeasibilities. */

static void btrack_most_feas(glp_tree *tree)
{     IOSNPD *node;
      int p;
      double best;
      p = 0, best = DBL_MAX;
      for (node = tree->head; node != NULL; node = node->next)
      {  xassert(node->up != NULL);
         if (best > node->up->ii_sum)
            p = node->p, best = node->up->ii_sum;
      }
      ios_revive_node(tree, p);
      return;
}

/*----------------------------------------------------------------------
-- btrack_best_proj - select subproblem with best projection heur.
--
-- This routine selects from the active list a subproblem to be solved
-- next using the best projection heuristic. */

static void btrack_best_proj(glp_tree *tree)
{     IOSNPD *root, *node;
      int p;
      double best, deg, obj;
      /* the global bound must exist */
      xassert(tree->mip->mip_stat == GLP_FEAS);
      /* obtain pointer to the root node, which must exist */
      root = tree->slot[1].node;
      xassert(root != NULL);
      /* deg estimates degradation of the objective function per unit
         of the sum of integer infeasibilities */
      xassert(root->ii_sum > 0.0);
      deg = (tree->mip->mip_obj - root->bound) / root->ii_sum;
      /* nothing has been selected so far */
      p = 0, best = DBL_MAX;
      /* walk through the list of active subproblems */
      for (node = tree->head; node != NULL; node = node->next)
      {  xassert(node->up != NULL);
         /* obj estimates optimal objective value if the sum of integer
            infeasibilities were zero */
         obj = node->up->bound + deg * node->up->ii_sum;
         if (tree->mip->dir == GLP_MAX) obj = - obj;
         /* select the subproblem which has the best estimated optimal
            objective value */
         if (best > obj) p = node->p, best = obj;
      }
      ios_revive_node(tree, p);
      return;
}

#undef int_col
#undef print

static void btrack_best_node(glp_tree *tree)
{     IOSNPD *node, *best = NULL;
      double bound, eps;
      switch (tree->mip->dir)
      {  case GLP_MIN:
            bound = +DBL_MAX;
            for (node = tree->head; node != NULL; node = node->next)
               if (bound > node->bound) bound = node->bound;
            xassert(bound != +DBL_MAX);
            eps = 0.001 * (1.0 + fabs(bound));
            for (node = tree->head; node != NULL; node = node->next)
            {  if (node->bound <= bound + eps)
               {  xassert(node->up != NULL);
                  if (best == NULL ||
#if 1
                  best->up->ii_sum > node->up->ii_sum) best = node;
#else
                  best->lp_obj > node->lp_obj) best = node;
#endif
               }
            }
            break;
         case GLP_MAX:
            bound = -DBL_MAX;
            for (node = tree->head; node != NULL; node = node->next)
               if (bound < node->bound) bound = node->bound;
            xassert(bound != -DBL_MAX);
            eps = 0.001 * (1.0 + fabs(bound));
            for (node = tree->head; node != NULL; node = node->next)
            {  if (node->bound >= bound - eps)
               {  xassert(node->up != NULL);
                  if (best == NULL ||
#if 1
                  best->up->ii_sum > node->up->ii_sum) best = node;
#else
                  best->lp_obj < node->lp_obj) best = node;
#endif
               }
            }
            break;
         default:
            xassert(tree != tree);
      }
      xassert(best != NULL);
      ios_revive_node(tree, best->p);
      return;
}

static void write_sol(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int term_out = lib_link_env()->term_out;
      if (tree->parm->msg_lev < GLP_MSG_ALL)
         glp_term_out(GLP_OFF);
      if (tree->parm->fn_sol != NULL)
      {  int m = mip->m;
         mip->m = tree->orig_m;
         glp_write_mip(mip, tree->parm->fn_sol);
         mip->m = m;
      }
      if (tree->parm->msg_lev < GLP_MSG_ALL)
         glp_term_out(term_out);
      return;
}

static void display_cut_info(glp_tree *tree)
{     /* display number of cuts added to current subproblem */
      glp_prob *mip = tree->mip;
      int i, gmi = 0, mir = 0, cov = 0, clq = 0, app = 0;
      for (i = mip->m; i > 0; i--)
      {  GLPROW *row;
         row = mip->row[i];
         /* if (row->level < tree->curr->level) break; */
         if (row->origin == GLP_RF_CUT)
         {  if (row->klass == GLP_RF_GMI)
               gmi++;
            else if (row->klass == GLP_RF_MIR)
               mir++;
            else if (row->klass == GLP_RF_COV)
               cov++;
            else if (row->klass == GLP_RF_CLQ)
               clq++;