glpapi09.c

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

/* glpapi09.c (mixed integer programming routines) */

/***********************************************************************
*  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"
#include "glpipp.h"

/***********************************************************************
*  NAME
*
*  glp_set_col_kind - set (change) column kind
*
*  SYNOPSIS
*
*  void glp_set_col_kind(glp_prob *mip, int j, int kind);
*
*  DESCRIPTION
*
*  The routine glp_set_col_kind sets (changes) the kind of j-th column
*  (structural variable) as specified by the parameter kind:
*
*  GLP_CV - continuous variable;
*  GLP_IV - integer variable;
*  GLP_BV - binary variable. */

void glp_set_col_kind(glp_prob *mip, int j, int kind)
{     GLPCOL *col;
      if (!(1 <= j && j <= mip->n))
         xerror("glp_set_col_kind: j = %d; column number out of range\n"
            , j);
      col = mip->col[j];
      switch (kind)
      {  case GLP_CV:
            col->kind = GLP_CV;
            break;
         case GLP_IV:
            col->kind = GLP_IV;
            break;
         case GLP_BV:
            col->kind = GLP_IV;
            if (!(col->type == GLP_DB && col->lb == 0.0 && col->ub ==
               1.0)) glp_set_col_bnds(mip, j, GLP_DB, 0.0, 1.0);
            break;
         default:
            xerror("glp_set_col_kind: j = %d; kind = %d; invalid column"
               " kind\n", j, kind);
      }
      return;
}

/***********************************************************************
*  NAME
*
*  glp_get_col_kind - retrieve column kind
*
*  SYNOPSIS
*
*  int glp_get_col_kind(glp_prob *mip, int j);
*
*  RETURNS
*
*  The routine glp_get_col_kind returns the kind of j-th column, i.e.
*  the kind of corresponding structural variable, as follows:
*
*  GLP_CV - continuous variable;
*  GLP_IV - integer variable;
*  GLP_BV - binary variable */

int glp_get_col_kind(glp_prob *mip, int j)
{     GLPCOL *col;
      int kind;
      if (!(1 <= j && j <= mip->n))
         xerror("glp_get_col_kind: j = %d; column number out of range\n"
            , j);
      col = mip->col[j];
      kind = col->kind;
      switch (kind)
      {  case GLP_CV:
            break;
         case GLP_IV:
            if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
               kind = GLP_BV;
            break;
         default:
            xassert(kind != kind);
      }
      return kind;
}

/***********************************************************************
*  NAME
*
*  glp_get_num_int - retrieve number of integer columns
*
*  SYNOPSIS
*
*  int glp_get_num_int(glp_prob *mip);
*
*  RETURNS
*
*  The routine glp_get_num_int returns the current number of columns,
*  which are marked as integer. */

int glp_get_num_int(glp_prob *mip)
{     GLPCOL *col;
      int j, count = 0;
      for (j = 1; j <= mip->n; j++)
      {  col = mip->col[j];
         if (col->kind == GLP_IV) count++;
      }
      return count;
}

/***********************************************************************
*  NAME
*
*  glp_get_num_bin - retrieve number of binary columns
*
*  SYNOPSIS
*
*  int glp_get_num_bin(glp_prob *mip);
*
*  RETURNS
*
*  The routine glp_get_num_bin returns the current number of columns,
*  which are marked as binary. */

int glp_get_num_bin(glp_prob *mip)
{     GLPCOL *col;
      int j, count = 0;
      for (j = 1; j <= mip->n; j++)
      {  col = mip->col[j];
         if (col->kind == GLP_IV && col->type == GLP_DB && col->lb ==
            0.0 && col->ub == 1.0) count++;
      }
      return count;
}

/***********************************************************************
*  NAME
*
*  glp_intopt - solve MIP problem with the branch-and-bound method
*
*  SYNOPSIS
*
*  int glp_intopt(glp_prob *mip, const glp_iocp *parm);
*
*  DESCRIPTION
*
*  The routine glp_intopt is a driver to the MIP solver based on the
*  branch-and-bound method.
*
*  On entry the problem object should contain optimal solution to LP
*  relaxation (which can be obtained with the routine glp_simplex).
*
*  The MIP solver has a set of control parameters. Values of the control
*  parameters can be passed in a structure glp_iocp, which the parameter
*  parm points to.
*
*  The parameter parm can be specified as NULL, in which case the MIP
*  solver uses default settings.
*
*  RETURNS
*
*  0  The MIP problem instance has been successfully solved. This code
*     does not necessarily mean that the solver has found optimal
*     solution. It only means that the solution process was successful.
*
*  GLP_EBOUND
*     Unable to start the search, because some double-bounded variables
*     have incorrect bounds or some integer variables have non-integer
*     (fractional) bounds.
*
*  GLP_EROOT
*     Unable to start the search, because optimal basis for initial LP
*     relaxation is not provided.
*
*  GLP_EFAIL
*     The search was prematurely terminated due to the solver failure.
*
*  GLP_EMIPGAP
*     The search was prematurely terminated, because the relative mip
*     gap tolerance has been reached.
*
*  GLP_ETMLIM
*     The search was prematurely terminated, because the time limit has
*     been exceeded.
*
*  GLP_ENOPFS
*     The MIP problem instance has no primal feasible solution (only if
*     the MIP presolver is used).
*
*  GLP_ENODFS
*     LP relaxation of the MIP problem instance has no dual feasible
*     solution (only if the MIP presolver is used).
*
*  GLP_ESTOP
*     The search was prematurely terminated by application. */

static int driver1(glp_prob *mip, const glp_iocp *parm)
{     /* base driver which does not use MIP presolver */
      glp_tree *tree;
      int ret;
      /* optimal solution to LP relaxation must be known */
      if (glp_get_status(mip) != GLP_OPT)
      {  if (parm->msg_lev >= GLP_MSG_ERR)
            xprintf("glp_intopt: optimal basis to initial LP relaxation"
               " not provided\n");
         ret = GLP_EROOT;
         goto done;
      }
      /* it seems all is ok */
      if (parm->msg_lev >= GLP_MSG_ALL)
         xprintf("Integer optimization begins...\n");
      /* create the branch-and-bound tree */
#if 0
      ((glp_iocp *)parm)->msg_lev = GLP_MSG_DBG;
#endif
      tree = ios_create_tree(mip, parm);
      /* try to solve the problem */
      ret = ios_driver(tree);
      /* analyze exit code reported by the mip driver */
      switch (ret)
      {  case 0:
            if (tree->mip->mip_stat == GLP_FEAS)
            {  if (parm->msg_lev >= GLP_MSG_ALL)
                  xprintf("INTEGER OPTIMAL SOLUTION FOUND\n");
               tree->mip->mip_stat = GLP_OPT;
            }
            else
            {  if (parm->msg_lev >= GLP_MSG_ALL)
                  xprintf("PROBLEM HAS NO INTEGER FEASIBLE SOLUTION\n");
               tree->mip->mip_stat = GLP_NOFEAS;
            }
            break;
         case GLP_EMIPGAP:
            if (parm->msg_lev >= GLP_MSG_ALL)
               xprintf("RELATIVE MIP GAP TOLERANCE REACHED; SEARCH TERM"
                  "INATED\n");
            break;
         case GLP_ETMLIM:
            if (parm->msg_lev >= GLP_MSG_ALL)
               xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
            break;
         case GLP_EFAIL:
            if (parm->msg_lev >= GLP_MSG_ERR)
               xprintf("glp_intopt: cannot solve current LP relaxation "
                  "\n");
            break;
         case GLP_ESTOP:
            if (parm->msg_lev >= GLP_MSG_ALL)
               xprintf("SEARCH TERMINATED BY APPLICATION\n");
            break;
         default:
            xassert(ret != ret);
      }
      /* delete the branch-and-bound tree */
      ios_delete_tree(tree);
done: return ret;
}

static int driver2(glp_prob *orig, const glp_iocp *parm)
{     /* extended driver which uses MIP presolver */
      LIBENV *env = lib_link_env();
      int term_out = env->term_out;
      IPP *ipp;
      glp_prob *prob = NULL;
      int j, i_stat, ret;
      /* the problem must have at least one row and one column */
      xassert(orig->m > 0 && orig->n > 0);
      /* reset the status of MIP solution */
      orig->mip_stat = GLP_UNDEF;
      /* create MIP presolver workspace */
      ipp = ipp_create_wksp();
      /* load the original problem into the presolver workspace */
      ipp_load_orig(ipp, orig);
      /* perform basic MIP presolve analysis */
      if (!term_out || parm->msg_lev < GLP_MSG_ALL)
         env->term_out = GLP_OFF;
      else
         env->term_out = GLP_ON;
      ret = ipp_basic_tech(ipp);
      env->term_out = term_out;
      switch (ret)
      {  case 0:
            /* no infeasibility is detected */
            break;
         case 1:
nopfs:      /* primal infeasibility is detected */
            if (parm->msg_lev >= GLP_MSG_ALL)
               xprintf("PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION\n");
            ret = GLP_ENOPFS;
            goto done;
         case 2:
nodfs:      /* dual infeasibility is detected */
            if (parm->msg_lev >= GLP_MSG_ALL)
               xprintf("LP RELAXATION HAS NO DUAL FEASIBLE SOLUTION\n");
            ret = GLP_ENODFS;
            goto done;
         default:
            xassert(ipp != ipp);
      }
      /* reduce column bounds */
      if (!term_out || parm->msg_lev < GLP_MSG_ALL)
         env->term_out = GLP_OFF;
      else
         env->term_out = GLP_ON;
      ret = ipp_reduce_bnds(ipp);
      env->term_out = term_out;
      switch (ret)
      {  case 0:  break;
         case 1:  goto nopfs;
         default: xassert(ipp != ipp);
      }
      /* perform basic MIP presolve analysis */
      if (!term_out || parm->msg_lev < GLP_MSG_ALL)
         env->term_out = GLP_OFF;
      else
         env->term_out = GLP_ON;
      ret = ipp_basic_tech(ipp);
      env->term_out = term_out;
      switch (ret)
      {  case 0:  break;
         case 1:  goto nopfs;
         case 2:  goto nodfs;
         default: xassert(ipp != ipp);
      }
      /* replace general integer variables by sum of binary variables,
         if required */
      if (parm->binarize)
      {  if (!term_out || parm->msg_lev < GLP_MSG_ALL)
            env->term_out = GLP_OFF;
         else
            env->term_out = GLP_ON;
         ipp_binarize(ipp);
         env->term_out = term_out;
      }
      /* perform coefficient reduction */
      if (!term_out || parm->msg_lev < GLP_MSG_ALL)
         env->term_out = GLP_OFF;
      else
         env->term_out = GLP_ON;
      ipp_reduction(ipp);
      env->term_out = term_out;
      /* if the resultant problem is empty, it has an empty solution,
         which is optimal */
      if (ipp->row_ptr == NULL || ipp->col_ptr == NULL)
      {  xassert(ipp->row_ptr == NULL);
         xassert(ipp->col_ptr == NULL);
         if (parm->msg_lev >= GLP_MSG_ALL)
         {  xprintf("Objective value = %.10g\n",
               ipp->orig_dir == LPX_MIN ? +ipp->c0 : -ipp->c0);
            xprintf("INTEGER OPTIMAL SOLUTION FOUND BY MIP PRESOLVER\n")
               ;
         }
         /* allocate recovered solution segment */
         ipp->col_stat = xcalloc(1+ipp->ncols, sizeof(int));
         ipp->col_mipx = xcalloc(1+ipp->ncols, sizeof(double));
         for (j = 1; j <= ipp->ncols; j++) ipp->col_stat[j] = 0;
         /* perform MIP postsolve processing */
         ipp_postsolve(ipp);
         /* unload recovered MIP solution and store it in the original
            problem object */
         ipp_unload_sol(ipp, orig, LPX_I_OPT);
         ret = 0;