1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
/* gmicut.c (Gomory's mixed integer cut generator) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* Copyright (C) 2002-2016 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 "prob.h"
/***********************************************************************
* NAME
*
* glp_gmi_cut - generate Gomory's mixed integer cut (core routine)
*
* SYNOPSIS
*
* int glp_gmi_cut(glp_prob *P, int j, int ind[], double val[], double
* phi[]);
*
* DESCRIPTION
*
* This routine attempts to generate a Gomory's mixed integer cut for
* specified integer column (structural variable), whose primal value
* in current basic solution is integer infeasible (fractional).
*
* On entry to the routine the basic solution contained in the problem
* object P should be optimal, and the basis factorization should be
* valid. The parameter j should specify the ordinal number of column
* (structural variable x[j]), for which the cut should be generated,
* 1 <= j <= n, where n is the number of columns in the problem object.
* This column should be integer, non-fixed, and basic, and its primal
* value should be fractional.
*
* The cut generated by the routine is the following inequality:
*
* sum a[j] * x[j] >= b,
*
* which is expected to be violated at the current basic solution.
*
* If the cut has been successfully generated, the routine stores its
* non-zero coefficients a[j] and corresponding column indices j in the
* array locations val[1], ..., val[len] and ind[1], ..., ind[len],
* where 1 <= len <= n is the number of non-zero coefficients. The
* right-hand side value b is stored in val[0], and ind[0] is set to 0.
*
* The working array phi should have 1+m+n locations (location phi[0]
* is not used), where m and n is the number of rows and columns in the
* problem object, resp.
*
* RETURNS
*
* If the cut has been successfully generated, the routine returns
* len, the number of non-zero coefficients in the cut, 1 <= len <= n.
*
* Otherwise, the routine returns one of the following codes:
*
* -1 current basis factorization is not valid;
*
* -2 current basic solution is not optimal;
*
* -3 column ordinal number j is out of range;
*
* -4 variable x[j] is not of integral kind;
*
* -5 variable x[j] is either fixed or non-basic;
*
* -6 primal value of variable x[j] in basic solution is too close
* to nearest integer;
*
* -7 some coefficients in the simplex table row corresponding to
* variable x[j] are too large in magnitude;
*
* -8 some free (unbounded) variables have non-zero coefficients in
* the simplex table row corresponding to variable x[j].
*
* ALGORITHM
*
* See glpk/doc/notes/gomory (in Russian). */
#define f(x) ((x) - floor(x))
/* compute fractional part of x */
int glp_gmi_cut(glp_prob *P, int j,
int ind[/*1+n*/], double val[/*1+n*/], double phi[/*1+m+n*/])
{ int m = P->m;
int n = P->n;
GLPROW *row;
GLPCOL *col;
GLPAIJ *aij;
int i, k, len, kind, stat;
double lb, ub, alfa, beta, ksi, phi1, rhs;
/* sanity checks */
if (!(P->m == 0 || P->valid))
{ /* current basis factorization is not valid */
return -1;
}
if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
{ /* current basic solution is not optimal */
return -2;
}
if (!(1 <= j && j <= n))
{ /* column ordinal number is out of range */
return -3;
}
col = P->col[j];
if (col->kind != GLP_IV)
{ /* x[j] is not of integral kind */
return -4;
}
if (col->type == GLP_FX || col->stat != GLP_BS)
{ /* x[j] is either fixed or non-basic */
return -5;
}
if (fabs(col->prim - floor(col->prim + 0.5)) < 0.001)
{ /* primal value of x[j] is too close to nearest integer */
return -6;
}
/* compute row of the simplex tableau, which (row) corresponds
* to specified basic variable xB[i] = x[j]; see (23) */
len = glp_eval_tab_row(P, 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 = P->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 are 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 */
row = P->row[k];
kind = GLP_CV;
lb = row->lb;
ub = row->ub;
stat = row->stat;
}
else
{ /* structural variable */
col = P->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 magnitude, report failure */
if (fabs(ksi) > 1e+05)
return -7;
/* if ksi[i,j] is too small in 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];
* report failure */
return -8;
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++)
{ if (fabs(phi[i]) < 1e-10)
continue;
/* auxiliary variable x[i] has non-zero cut coefficient */
row = P->row[i];
/* x[i] cannot be fixed variable */
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++)
{ if (fabs(phi[m+j]) < 1e-10)
continue;
/* structural variable x[m+j] has non-zero cut coefficient */
col = P->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;
ind[0] = 0, val[0] = rhs;
/* the cut has been successfully generated */
return len;
}
/* eof */
|
