You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
282 lines
9.9 KiB
282 lines
9.9 KiB
/* glpios05.c (Gomory's mixed integer cut generator) */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
|
* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* 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 "env.h"
|
|
#include "glpios.h"
|
|
|
|
/***********************************************************************
|
|
* NAME
|
|
*
|
|
* ios_gmi_gen - generate Gomory's mixed integer cuts.
|
|
*
|
|
* SYNOPSIS
|
|
*
|
|
* #include "glpios.h"
|
|
* void ios_gmi_gen(glp_tree *tree, IOSPOOL *pool);
|
|
*
|
|
* DESCRIPTION
|
|
*
|
|
* The routine ios_gmi_gen generates Gomory's mixed integer cuts for
|
|
* the current point and adds them to the cut pool. */
|
|
|
|
#define MAXCUTS 50
|
|
/* maximal number of cuts to be generated for one round */
|
|
|
|
struct worka
|
|
{ /* Gomory's cut generator working area */
|
|
int *ind; /* int ind[1+n]; */
|
|
double *val; /* double val[1+n]; */
|
|
double *phi; /* double phi[1+m+n]; */
|
|
};
|
|
|
|
#define f(x) ((x) - floor(x))
|
|
/* compute fractional part of x */
|
|
|
|
static void gen_cut(glp_tree *tree, struct worka *worka, int j)
|
|
{ /* this routine tries to generate Gomory's mixed integer cut for
|
|
specified structural variable x[m+j] of integer kind, which is
|
|
basic and has fractional value in optimal solution to current
|
|
LP relaxation */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mip->m;
|
|
int n = mip->n;
|
|
int *ind = worka->ind;
|
|
double *val = worka->val;
|
|
double *phi = worka->phi;
|
|
int i, k, len, kind, stat;
|
|
double lb, ub, alfa, beta, ksi, phi1, rhs;
|
|
/* compute row of the simplex tableau, which (row) corresponds
|
|
to specified basic variable xB[i] = x[m+j]; see (23) */
|
|
len = glp_eval_tab_row(mip, m+j, ind, val);
|
|
/* determine beta[i], which a value of xB[i] in optimal solution
|
|
to current LP relaxation; note that this value is the same as
|
|
if it would be computed with formula (27); it is assumed that
|
|
beta[i] is fractional enough */
|
|
beta = mip->col[j]->prim;
|
|
/* compute cut coefficients phi and right-hand side rho, which
|
|
correspond to formula (30); dense format is used, because rows
|
|
of the simplex tableau is usually dense */
|
|
for (k = 1; k <= m+n; k++) phi[k] = 0.0;
|
|
rhs = f(beta); /* initial value of rho; see (28), (32) */
|
|
for (j = 1; j <= len; j++)
|
|
{ /* determine original number of non-basic variable xN[j] */
|
|
k = ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
/* determine the kind, bounds and current status of xN[j] in
|
|
optimal solution to LP relaxation */
|
|
if (k <= m)
|
|
{ /* auxiliary variable */
|
|
GLPROW *row = mip->row[k];
|
|
kind = GLP_CV;
|
|
lb = row->lb;
|
|
ub = row->ub;
|
|
stat = row->stat;
|
|
}
|
|
else
|
|
{ /* structural variable */
|
|
GLPCOL *col = mip->col[k-m];
|
|
kind = col->kind;
|
|
lb = col->lb;
|
|
ub = col->ub;
|
|
stat = col->stat;
|
|
}
|
|
/* xN[j] cannot be basic */
|
|
xassert(stat != GLP_BS);
|
|
/* determine row coefficient ksi[i,j] at xN[j]; see (23) */
|
|
ksi = val[j];
|
|
/* if ksi[i,j] is too large in the magnitude, do not generate
|
|
the cut */
|
|
if (fabs(ksi) > 1e+05) goto fini;
|
|
/* if ksi[i,j] is too small in the magnitude, skip it */
|
|
if (fabs(ksi) < 1e-10) goto skip;
|
|
/* compute row coefficient alfa[i,j] at y[j]; see (26) */
|
|
switch (stat)
|
|
{ case GLP_NF:
|
|
/* xN[j] is free (unbounded) having non-zero ksi[i,j];
|
|
do not generate the cut */
|
|
goto fini;
|
|
case GLP_NL:
|
|
/* xN[j] has active lower bound */
|
|
alfa = - ksi;
|
|
break;
|
|
case GLP_NU:
|
|
/* xN[j] has active upper bound */
|
|
alfa = + ksi;
|
|
break;
|
|
case GLP_NS:
|
|
/* xN[j] is fixed; skip it */
|
|
goto skip;
|
|
default:
|
|
xassert(stat != stat);
|
|
}
|
|
/* compute cut coefficient phi'[j] at y[j]; see (21), (28) */
|
|
switch (kind)
|
|
{ case GLP_IV:
|
|
/* y[j] is integer */
|
|
if (fabs(alfa - floor(alfa + 0.5)) < 1e-10)
|
|
{ /* alfa[i,j] is close to nearest integer; skip it */
|
|
goto skip;
|
|
}
|
|
else if (f(alfa) <= f(beta))
|
|
phi1 = f(alfa);
|
|
else
|
|
phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa));
|
|
break;
|
|
case GLP_CV:
|
|
/* y[j] is continuous */
|
|
if (alfa >= 0.0)
|
|
phi1 = + alfa;
|
|
else
|
|
phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa);
|
|
break;
|
|
default:
|
|
xassert(kind != kind);
|
|
}
|
|
/* compute cut coefficient phi[j] at xN[j] and update right-
|
|
hand side rho; see (31), (32) */
|
|
switch (stat)
|
|
{ case GLP_NL:
|
|
/* xN[j] has active lower bound */
|
|
phi[k] = + phi1;
|
|
rhs += phi1 * lb;
|
|
break;
|
|
case GLP_NU:
|
|
/* xN[j] has active upper bound */
|
|
phi[k] = - phi1;
|
|
rhs -= phi1 * ub;
|
|
break;
|
|
default:
|
|
xassert(stat != stat);
|
|
}
|
|
skip: ;
|
|
}
|
|
/* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut
|
|
coefficients are stored in the array phi in dense format;
|
|
x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc-
|
|
tural variables; see (30) */
|
|
/* eliminate auxiliary variables in order to express the cut only
|
|
through structural variables; see (33) */
|
|
for (i = 1; i <= m; i++)
|
|
{ GLPROW *row;
|
|
GLPAIJ *aij;
|
|
if (fabs(phi[i]) < 1e-10) continue;
|
|
/* auxiliary variable x[i] has non-zero cut coefficient */
|
|
row = mip->row[i];
|
|
/* x[i] cannot be fixed */
|
|
xassert(row->type != GLP_FX);
|
|
/* substitute x[i] = sum_j a[i,j] * x[m+j] */
|
|
for (aij = row->ptr; aij != NULL; aij = aij->r_next)
|
|
phi[m+aij->col->j] += phi[i] * aij->val;
|
|
}
|
|
/* convert the final cut to sparse format and substitute fixed
|
|
(structural) variables */
|
|
len = 0;
|
|
for (j = 1; j <= n; j++)
|
|
{ GLPCOL *col;
|
|
if (fabs(phi[m+j]) < 1e-10) continue;
|
|
/* structural variable x[m+j] has non-zero cut coefficient */
|
|
col = mip->col[j];
|
|
if (col->type == GLP_FX)
|
|
{ /* eliminate x[m+j] */
|
|
rhs -= phi[m+j] * col->lb;
|
|
}
|
|
else
|
|
{ len++;
|
|
ind[len] = j;
|
|
val[len] = phi[m+j];
|
|
}
|
|
}
|
|
if (fabs(rhs) < 1e-12) rhs = 0.0;
|
|
/* if the cut inequality seems to be badly scaled, reject it to
|
|
avoid numeric difficulties */
|
|
for (k = 1; k <= len; k++)
|
|
{ if (fabs(val[k]) < 1e-03) goto fini;
|
|
if (fabs(val[k]) > 1e+03) goto fini;
|
|
}
|
|
/* add the cut to the cut pool for further consideration */
|
|
#if 0
|
|
ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO,
|
|
rhs);
|
|
#else
|
|
glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO,
|
|
rhs);
|
|
#endif
|
|
fini: return;
|
|
}
|
|
|
|
struct var { int j; double f; };
|
|
|
|
static int fcmp(const void *p1, const void *p2)
|
|
{ const struct var *v1 = p1, *v2 = p2;
|
|
if (v1->f > v2->f) return -1;
|
|
if (v1->f < v2->f) return +1;
|
|
return 0;
|
|
}
|
|
|
|
void ios_gmi_gen(glp_tree *tree)
|
|
{ /* main routine to generate Gomory's cuts */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mip->m;
|
|
int n = mip->n;
|
|
struct var *var;
|
|
int k, nv, j, size;
|
|
struct worka _worka, *worka = &_worka;
|
|
/* allocate working arrays */
|
|
var = xcalloc(1+n, sizeof(struct var));
|
|
worka->ind = xcalloc(1+n, sizeof(int));
|
|
worka->val = xcalloc(1+n, sizeof(double));
|
|
worka->phi = xcalloc(1+m+n, sizeof(double));
|
|
/* build the list of integer structural variables, which are
|
|
basic and have fractional value in optimal solution to current
|
|
LP relaxation */
|
|
nv = 0;
|
|
for (j = 1; j <= n; j++)
|
|
{ GLPCOL *col = mip->col[j];
|
|
double frac;
|
|
if (col->kind != GLP_IV) continue;
|
|
if (col->type == GLP_FX) continue;
|
|
if (col->stat != GLP_BS) continue;
|
|
frac = f(col->prim);
|
|
if (!(0.05 <= frac && frac <= 0.95)) continue;
|
|
/* add variable to the list */
|
|
nv++, var[nv].j = j, var[nv].f = frac;
|
|
}
|
|
/* order the list by descending fractionality */
|
|
qsort(&var[1], nv, sizeof(struct var), fcmp);
|
|
/* try to generate cuts by one for each variable in the list, but
|
|
not more than MAXCUTS cuts */
|
|
size = glp_ios_pool_size(tree);
|
|
for (k = 1; k <= nv; k++)
|
|
{ if (glp_ios_pool_size(tree) - size >= MAXCUTS) break;
|
|
gen_cut(tree, worka, var[k].j);
|
|
}
|
|
/* free working arrays */
|
|
xfree(var);
|
|
xfree(worka->ind);
|
|
xfree(worka->val);
|
|
xfree(worka->phi);
|
|
return;
|
|
}
|
|
|
|
/* eof */
|