glpios03.c

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

/* glpios03.c (under construction so far) */

/***********************************************************************
*  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/>.
***********************************************************************/

#define _GLPSTD_STDIO
#include "glpios.h"

/***********************************************************************
*  show_progress - display progress of the search
*
*  This routine displays some information about progress of the search.
*
*  The information includes:
*
*  the current number of iterations performed by the simplex solver;
*
*  the objective value for the best known integer feasible solution,
*  which is upper (minimization) or lower (maximization) global bound
*  for optimal solution of the original mip problem;
*
*  the best local bound for active nodes, which is lower (minimization)
*  or upper (maximization) global bound for optimal solution of the
*  original mip problem;
*
*  the relative mip gap, in percents;
*
*  the number of open (active) subproblems;
*
*  the number of completely explored subproblems, i.e. whose nodes have
*  been removed from the tree. */

static void show_progress(glp_tree *tree, int bingo)
{     int p;
      double temp;
      char best_mip[50], best_bound[50], *rho, rel_gap[50];
      /* format the best known integer feasible solution */
      if (tree->mip->mip_stat == GLP_FEAS)
         sprintf(best_mip, "%17.9e", tree->mip->mip_obj);
      else
         sprintf(best_mip, "%17s", "not found yet");
      /* determine reference number of an active subproblem whose local
         bound is best */
      p = ios_best_node(tree);
      /* format the best bound */
      if (p == 0)
         sprintf(best_bound, "%17s", "tree is empty");
      else
      {  temp = tree->slot[p].node->bound;
         if (temp == -DBL_MAX)
            sprintf(best_bound, "%17s", "-inf");
         else if (temp == +DBL_MAX)
            sprintf(best_bound, "%17s", "+inf");
         else
            sprintf(best_bound, "%17.9e", temp);
      }
      /* choose the relation sign between global bounds */
      switch (tree->mip->dir)
      {  case GLP_MIN: rho = ">="; break;
         case GLP_MAX: rho = "<="; break;
         default: xassert(tree != tree);
      }
      /* format the relative mip gap */
      temp = ios_relative_gap(tree);
      if (temp == 0.0)
         sprintf(rel_gap, "  0.0%%");
      else if (temp < 0.001)
         sprintf(rel_gap, "< 0.1%%");
      else if (temp <= 9.999)
         sprintf(rel_gap, "%5.1f%%", 100.0 * temp);
      else
         sprintf(rel_gap, "%6s", "");
      /* display progress of the search */
      xprintf("+%6d: %s %s %s %s %s (%d; %d)\n",
         tree->mip->it_cnt, bingo ? ">>>>>" : "mip =", best_mip, rho,
         best_bound, rel_gap, tree->a_cnt, tree->t_cnt - tree->n_cnt);
      tree->tm_lag = xtime();
      return;
}

/***********************************************************************
*  is_branch_hopeful - check if specified branch is hopeful
*
*  This routine checks if the specified subproblem can have an integer
*  optimal solution which is better than the best known one.
*
*  The check is based on comparison of the local objective bound stored
*  in the subproblem descriptor and the incumbent objective value which
*  is the global objective bound.
*
*  If there is a chance that the specified subproblem can have a better
*  integer optimal solution, the routine returns non-zero. Otherwise, if
*  the corresponding branch can pruned, zero is returned. */

static int is_branch_hopeful(glp_tree *tree, int p)
{     xassert(1 <= p && p <= tree->nslots);
      xassert(tree->slot[p].node != NULL);
      return ios_is_hopeful(tree, tree->slot[p].node->bound);
}

static int is_curr_node_hopeful(glp_tree *tree)
{     return
         ios_is_hopeful(tree, tree->curr->bound);
}

/***********************************************************************
*  check_integrality - check integrality of basic solution
*
*  This routine checks if the basic solution of LP relaxation of the
*  current subproblem satisfies to integrality conditions, i.e. that all
*  variables of integer kind have integral primal values. (The solution
*  is assumed to be optimal.)
*
*  For each variable of integer kind the routine computes the following
*  quantity:
*
*     ii(x[j]) = min(x[j] - floor(x[j]), ceil(x[j]) - x[j]),         (1)
*
*  which is a measure of the integer infeasibility (non-integrality) of
*  x[j] (for example, ii(2.1) = 0.1, ii(3.7) = 0.3, ii(5.0) = 0). It is
*  understood that 0 <= ii(x[j]) <= 0.5, and variable x[j] is integer
*  feasible if ii(x[j]) = 0. However, due to floating-point arithmetic
*  the routine checks less restrictive condition:
*
*     ii(x[j]) <= tol_int,                                           (2)
*
*  where tol_int is a given tolerance (small positive number) and marks
*  each variable which does not satisfy to (2) as integer infeasible by
*  setting its fractionality flag.
*
*  In order to characterize integer infeasibility of the basic solution
*  in the whole the routine computes two parameters: ii_cnt, which is
*  the number of variables with the fractionality flag set, and ii_sum,
*  which is the sum of integer infeasibilities (1). */

static void check_integrality(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int j, type, ii_cnt = 0;
      double lb, ub, x, temp1, temp2, ii_sum = 0.0;
      /* walk through the set of columns (structural variables) */
      for (j = 1; j <= mip->n; j++)
      {  GLPCOL *col = mip->col[j];
         tree->non_int[j] = 0;
         /* if the column is not integer, skip it */
         if (col->kind != GLP_IV) continue;
         /* if the column is non-basic, it is integer feasible */
         if (col->stat != GLP_BS) continue;
         /* obtain the type and bounds of the column */
         type = col->type, lb = col->lb, ub = col->ub;
         /* obtain value of the column in optimal basic solution */
         x = col->prim;
         /* if the column's primal value is close to the lower bound,
            the column is integer feasible within given tolerance */
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
         {  temp1 = lb - tree->parm->tol_int;
            temp2 = lb + tree->parm->tol_int;
            if (temp1 <= x && x <= temp2) continue;
#if 0
            /* the lower bound must not be violated */
            xassert(x >= lb);
#else
            if (x < lb) continue;
#endif
         }
         /* if the column's primal value is close to the upper bound,
            the column is integer feasible within given tolerance */
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
         {  temp1 = ub - tree->parm->tol_int;
            temp2 = ub + tree->parm->tol_int;
            if (temp1 <= x && x <= temp2) continue;
#if 0
            /* the upper bound must not be violated */
            xassert(x <= ub);
#else
            if (x > ub) continue;
#endif
         }
         /* if the column's primal value is close to nearest integer,
            the column is integer feasible within given tolerance */
         temp1 = floor(x + 0.5) - tree->parm->tol_int;
         temp2 = floor(x + 0.5) + tree->parm->tol_int;
         if (temp1 <= x && x <= temp2) continue;
         /* otherwise the column is integer infeasible */
         tree->non_int[j] = 1;
         /* increase the number of fractional-valued columns */
         ii_cnt++;
         /* compute the sum of integer infeasibilities */
         temp1 = x - floor(x);
         temp2 = ceil(x) - x;
         xassert(temp1 > 0.0 && temp2 > 0.0);
         ii_sum += (temp1 <= temp2 ? temp1 : temp2);
      }
      /* store ii_cnt and ii_sum in the current problem descriptor */
      xassert(tree->curr != NULL);
      tree->curr->ii_cnt = ii_cnt;
      tree->curr->ii_sum = ii_sum;
      /* and also display these parameters */
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
      {  if (ii_cnt == 0)
            xprintf("There are no fractional columns\n");
         else if (ii_cnt == 1)
            xprintf("There is one fractional column, integer infeasibil"
               "ity is %.3e\n", ii_sum);
         else
            xprintf("There are %d fractional columns, integer infeasibi"
               "lity is %.3e\n", ii_cnt, ii_sum);
      }
      return;
}

/***********************************************************************
*  record_solution - record better integer feasible solution
*
*  This routine records optimal basic solution of LP relaxation of the
*  current subproblem, which being integer feasible is better than the
*  best known integer feasible solution. */

static void record_solution(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int i, j;
      mip->mip_stat = GLP_FEAS;
      mip->mip_obj = mip->obj_val;
      for (i = 1; i <= mip->m; i++)
      {  GLPROW *row = mip->row[i];
         row->mipx = row->prim;
      }
      for (j = 1; j <= mip->n; j++)
      {  GLPCOL *col = mip->col[j];
         switch (col->kind)
         {  case GLP_CV:
               col->mipx = col->prim;
               break;
            case GLP_IV:
               /* value of the integer column must be integral */
               col->mipx = floor(col->prim + 0.5);
               break;
            default:
               xassert(col != col);
         }
      }
      tree->sol_cnt++;
      return;
}

/***********************************************************************
*  fix_by_red_cost - fix non-basic integer columns by reduced costs
*
*  This routine fixes some non-basic integer columns if their reduced
*  costs indicate that increasing (decreasing) the column at least by
*  one involves the objective value becoming worse than the incumbent
*  objective value (i.e. the global bound). */

static void fix_by_red_cost(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int j, stat, fixed = 0;
      double obj, lb, ub, dj;
      /* the global bound must exist */
      xassert(tree->mip->mip_stat == GLP_FEAS);
      /* basic solution of LP relaxation must be optimal */
      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
      /* determine the objective function value */
      obj = mip->obj_val;
      /* walk through the column list */
      for (j = 1; j <= mip->n; j++)
      {  GLPCOL *col = mip->col[j];
         /* if the column is not integer, skip it */
         if (col->kind != GLP_IV) continue;
         /* obtain bounds of j-th column */
         lb = col->lb, ub = col->ub;
         /* and determine its status and reduced cost */
         stat = col->stat, dj = col->dual;
         /* analyze the reduced cost */
         switch (mip->dir)
         {  case GLP_MIN:
               /* minimization */
               if (stat == GLP_NL)
               {  /* j-th column is non-basic on its lower bound */
                  if (dj < 0.0) dj = 0.0;
                  if (obj + dj >= mip->mip_obj)
                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
               }
               else if (stat == GLP_NU)
               {  /* j-th column is non-basic on its upper bound */
                  if (dj > 0.0) dj = 0.0;
                  if (obj - dj >= mip->mip_obj)
                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
               }
               break;
            case GLP_MAX:
               /* maximization */
               if (stat == GLP_NL)
               {  /* j-th column is non-basic on its lower bound */
                  if (dj > 0.0) dj = 0.0;
                  if (obj + dj <= mip->mip_obj)
                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
               }
               else if (stat == GLP_NU)
               {  /* j-th column is non-basic on its upper bound */
                  if (dj < 0.0) dj = 0.0;
                  if (obj - dj <= mip->mip_obj)
                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
               }
               break;
            default:
               xassert(tree != tree);
         }
      }
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
      {  if (fixed == 0)
            /* nothing to say */;
         else if (fixed == 1)
            xprintf("One column has been fixed by reduced cost\n");
         else
            xprintf("%d columns have been fixed by reduced costs\n",
               fixed);
      }
      /* fixing non-basic columns on their current bounds must not
         change the basic solution */
      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
      return;
}

/**********************************************************************/
/**********************************************************************/
/**********************************************************************/

/*----------------------------------------------------------------------
-- branch_first - choose first branching variable.
--
-- This routine looks up the list of structural variables and chooses
-- the first one, which is of integer kind and has fractional value in
-- optimal solution of the current LP relaxation.
--