1 | /* glpios05.c (Gomory's mixed integer cut generator) */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glpios.h" |
---|
26 | |
---|
27 | /*********************************************************************** |
---|
28 | * NAME |
---|
29 | * |
---|
30 | * ios_gmi_gen - generate Gomory's mixed integer cuts. |
---|
31 | * |
---|
32 | * SYNOPSIS |
---|
33 | * |
---|
34 | * #include "glpios.h" |
---|
35 | * void ios_gmi_gen(glp_tree *tree, IOSPOOL *pool); |
---|
36 | * |
---|
37 | * DESCRIPTION |
---|
38 | * |
---|
39 | * The routine ios_gmi_gen generates Gomory's mixed integer cuts for |
---|
40 | * the current point and adds them to the cut pool. */ |
---|
41 | |
---|
42 | #define MAXCUTS 50 |
---|
43 | /* maximal number of cuts to be generated for one round */ |
---|
44 | |
---|
45 | struct worka |
---|
46 | { /* Gomory's cut generator working area */ |
---|
47 | int *ind; /* int ind[1+n]; */ |
---|
48 | double *val; /* double val[1+n]; */ |
---|
49 | double *phi; /* double phi[1+m+n]; */ |
---|
50 | }; |
---|
51 | |
---|
52 | #define f(x) ((x) - floor(x)) |
---|
53 | /* compute fractional part of x */ |
---|
54 | |
---|
55 | static void gen_cut(glp_tree *tree, struct worka *worka, int j) |
---|
56 | { /* this routine tries to generate Gomory's mixed integer cut for |
---|
57 | specified structural variable x[m+j] of integer kind, which is |
---|
58 | basic and has fractional value in optimal solution to current |
---|
59 | LP relaxation */ |
---|
60 | glp_prob *mip = tree->mip; |
---|
61 | int m = mip->m; |
---|
62 | int n = mip->n; |
---|
63 | int *ind = worka->ind; |
---|
64 | double *val = worka->val; |
---|
65 | double *phi = worka->phi; |
---|
66 | int i, k, len, kind, stat; |
---|
67 | double lb, ub, alfa, beta, ksi, phi1, rhs; |
---|
68 | /* compute row of the simplex tableau, which (row) corresponds |
---|
69 | to specified basic variable xB[i] = x[m+j]; see (23) */ |
---|
70 | len = glp_eval_tab_row(mip, m+j, ind, val); |
---|
71 | /* determine beta[i], which a value of xB[i] in optimal solution |
---|
72 | to current LP relaxation; note that this value is the same as |
---|
73 | if it would be computed with formula (27); it is assumed that |
---|
74 | beta[i] is fractional enough */ |
---|
75 | beta = mip->col[j]->prim; |
---|
76 | /* compute cut coefficients phi and right-hand side rho, which |
---|
77 | correspond to formula (30); dense format is used, because rows |
---|
78 | of the simplex tableau is usually dense */ |
---|
79 | for (k = 1; k <= m+n; k++) phi[k] = 0.0; |
---|
80 | rhs = f(beta); /* initial value of rho; see (28), (32) */ |
---|
81 | for (j = 1; j <= len; j++) |
---|
82 | { /* determine original number of non-basic variable xN[j] */ |
---|
83 | k = ind[j]; |
---|
84 | xassert(1 <= k && k <= m+n); |
---|
85 | /* determine the kind, bounds and current status of xN[j] in |
---|
86 | optimal solution to LP relaxation */ |
---|
87 | if (k <= m) |
---|
88 | { /* auxiliary variable */ |
---|
89 | GLPROW *row = mip->row[k]; |
---|
90 | kind = GLP_CV; |
---|
91 | lb = row->lb; |
---|
92 | ub = row->ub; |
---|
93 | stat = row->stat; |
---|
94 | } |
---|
95 | else |
---|
96 | { /* structural variable */ |
---|
97 | GLPCOL *col = mip->col[k-m]; |
---|
98 | kind = col->kind; |
---|
99 | lb = col->lb; |
---|
100 | ub = col->ub; |
---|
101 | stat = col->stat; |
---|
102 | } |
---|
103 | /* xN[j] cannot be basic */ |
---|
104 | xassert(stat != GLP_BS); |
---|
105 | /* determine row coefficient ksi[i,j] at xN[j]; see (23) */ |
---|
106 | ksi = val[j]; |
---|
107 | /* if ksi[i,j] is too large in the magnitude, do not generate |
---|
108 | the cut */ |
---|
109 | if (fabs(ksi) > 1e+05) goto fini; |
---|
110 | /* if ksi[i,j] is too small in the magnitude, skip it */ |
---|
111 | if (fabs(ksi) < 1e-10) goto skip; |
---|
112 | /* compute row coefficient alfa[i,j] at y[j]; see (26) */ |
---|
113 | switch (stat) |
---|
114 | { case GLP_NF: |
---|
115 | /* xN[j] is free (unbounded) having non-zero ksi[i,j]; |
---|
116 | do not generate the cut */ |
---|
117 | goto fini; |
---|
118 | case GLP_NL: |
---|
119 | /* xN[j] has active lower bound */ |
---|
120 | alfa = - ksi; |
---|
121 | break; |
---|
122 | case GLP_NU: |
---|
123 | /* xN[j] has active upper bound */ |
---|
124 | alfa = + ksi; |
---|
125 | break; |
---|
126 | case GLP_NS: |
---|
127 | /* xN[j] is fixed; skip it */ |
---|
128 | goto skip; |
---|
129 | default: |
---|
130 | xassert(stat != stat); |
---|
131 | } |
---|
132 | /* compute cut coefficient phi'[j] at y[j]; see (21), (28) */ |
---|
133 | switch (kind) |
---|
134 | { case GLP_IV: |
---|
135 | /* y[j] is integer */ |
---|
136 | if (fabs(alfa - floor(alfa + 0.5)) < 1e-10) |
---|
137 | { /* alfa[i,j] is close to nearest integer; skip it */ |
---|
138 | goto skip; |
---|
139 | } |
---|
140 | else if (f(alfa) <= f(beta)) |
---|
141 | phi1 = f(alfa); |
---|
142 | else |
---|
143 | phi1 = (f(beta) / (1.0 - f(beta))) * (1.0 - f(alfa)); |
---|
144 | break; |
---|
145 | case GLP_CV: |
---|
146 | /* y[j] is continuous */ |
---|
147 | if (alfa >= 0.0) |
---|
148 | phi1 = + alfa; |
---|
149 | else |
---|
150 | phi1 = (f(beta) / (1.0 - f(beta))) * (- alfa); |
---|
151 | break; |
---|
152 | default: |
---|
153 | xassert(kind != kind); |
---|
154 | } |
---|
155 | /* compute cut coefficient phi[j] at xN[j] and update right- |
---|
156 | hand side rho; see (31), (32) */ |
---|
157 | switch (stat) |
---|
158 | { case GLP_NL: |
---|
159 | /* xN[j] has active lower bound */ |
---|
160 | phi[k] = + phi1; |
---|
161 | rhs += phi1 * lb; |
---|
162 | break; |
---|
163 | case GLP_NU: |
---|
164 | /* xN[j] has active upper bound */ |
---|
165 | phi[k] = - phi1; |
---|
166 | rhs -= phi1 * ub; |
---|
167 | break; |
---|
168 | default: |
---|
169 | xassert(stat != stat); |
---|
170 | } |
---|
171 | skip: ; |
---|
172 | } |
---|
173 | /* now the cut has the form sum_k phi[k] * x[k] >= rho, where cut |
---|
174 | coefficients are stored in the array phi in dense format; |
---|
175 | x[1,...,m] are auxiliary variables, x[m+1,...,m+n] are struc- |
---|
176 | tural variables; see (30) */ |
---|
177 | /* eliminate auxiliary variables in order to express the cut only |
---|
178 | through structural variables; see (33) */ |
---|
179 | for (i = 1; i <= m; i++) |
---|
180 | { GLPROW *row; |
---|
181 | GLPAIJ *aij; |
---|
182 | if (fabs(phi[i]) < 1e-10) continue; |
---|
183 | /* auxiliary variable x[i] has non-zero cut coefficient */ |
---|
184 | row = mip->row[i]; |
---|
185 | /* x[i] cannot be fixed */ |
---|
186 | xassert(row->type != GLP_FX); |
---|
187 | /* substitute x[i] = sum_j a[i,j] * x[m+j] */ |
---|
188 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
189 | phi[m+aij->col->j] += phi[i] * aij->val; |
---|
190 | } |
---|
191 | /* convert the final cut to sparse format and substitute fixed |
---|
192 | (structural) variables */ |
---|
193 | len = 0; |
---|
194 | for (j = 1; j <= n; j++) |
---|
195 | { GLPCOL *col; |
---|
196 | if (fabs(phi[m+j]) < 1e-10) continue; |
---|
197 | /* structural variable x[m+j] has non-zero cut coefficient */ |
---|
198 | col = mip->col[j]; |
---|
199 | if (col->type == GLP_FX) |
---|
200 | { /* eliminate x[m+j] */ |
---|
201 | rhs -= phi[m+j] * col->lb; |
---|
202 | } |
---|
203 | else |
---|
204 | { len++; |
---|
205 | ind[len] = j; |
---|
206 | val[len] = phi[m+j]; |
---|
207 | } |
---|
208 | } |
---|
209 | if (fabs(rhs) < 1e-12) rhs = 0.0; |
---|
210 | /* if the cut inequality seems to be badly scaled, reject it to |
---|
211 | avoid numeric difficulties */ |
---|
212 | for (k = 1; k <= len; k++) |
---|
213 | { if (fabs(val[k]) < 1e-03) goto fini; |
---|
214 | if (fabs(val[k]) > 1e+03) goto fini; |
---|
215 | } |
---|
216 | /* add the cut to the cut pool for further consideration */ |
---|
217 | #if 0 |
---|
218 | ios_add_cut_row(tree, pool, GLP_RF_GMI, len, ind, val, GLP_LO, |
---|
219 | rhs); |
---|
220 | #else |
---|
221 | glp_ios_add_row(tree, NULL, GLP_RF_GMI, 0, len, ind, val, GLP_LO, |
---|
222 | rhs); |
---|
223 | #endif |
---|
224 | fini: return; |
---|
225 | } |
---|
226 | |
---|
227 | struct var { int j; double f; }; |
---|
228 | |
---|
229 | static int fcmp(const void *p1, const void *p2) |
---|
230 | { const struct var *v1 = p1, *v2 = p2; |
---|
231 | if (v1->f > v2->f) return -1; |
---|
232 | if (v1->f < v2->f) return +1; |
---|
233 | return 0; |
---|
234 | } |
---|
235 | |
---|
236 | void ios_gmi_gen(glp_tree *tree) |
---|
237 | { /* main routine to generate Gomory's cuts */ |
---|
238 | glp_prob *mip = tree->mip; |
---|
239 | int m = mip->m; |
---|
240 | int n = mip->n; |
---|
241 | struct var *var; |
---|
242 | int k, nv, j, size; |
---|
243 | struct worka _worka, *worka = &_worka; |
---|
244 | /* allocate working arrays */ |
---|
245 | var = xcalloc(1+n, sizeof(struct var)); |
---|
246 | worka->ind = xcalloc(1+n, sizeof(int)); |
---|
247 | worka->val = xcalloc(1+n, sizeof(double)); |
---|
248 | worka->phi = xcalloc(1+m+n, sizeof(double)); |
---|
249 | /* build the list of integer structural variables, which are |
---|
250 | basic and have fractional value in optimal solution to current |
---|
251 | LP relaxation */ |
---|
252 | nv = 0; |
---|
253 | for (j = 1; j <= n; j++) |
---|
254 | { GLPCOL *col = mip->col[j]; |
---|
255 | double frac; |
---|
256 | if (col->kind != GLP_IV) continue; |
---|
257 | if (col->type == GLP_FX) continue; |
---|
258 | if (col->stat != GLP_BS) continue; |
---|
259 | frac = f(col->prim); |
---|
260 | if (!(0.05 <= frac && frac <= 0.95)) continue; |
---|
261 | /* add variable to the list */ |
---|
262 | nv++, var[nv].j = j, var[nv].f = frac; |
---|
263 | } |
---|
264 | /* order the list by descending fractionality */ |
---|
265 | qsort(&var[1], nv, sizeof(struct var), fcmp); |
---|
266 | /* try to generate cuts by one for each variable in the list, but |
---|
267 | not more than MAXCUTS cuts */ |
---|
268 | size = glp_ios_pool_size(tree); |
---|
269 | for (k = 1; k <= nv; k++) |
---|
270 | { if (glp_ios_pool_size(tree) - size >= MAXCUTS) break; |
---|
271 | gen_cut(tree, worka, var[k].j); |
---|
272 | } |
---|
273 | /* free working arrays */ |
---|
274 | xfree(var); |
---|
275 | xfree(worka->ind); |
---|
276 | xfree(worka->val); |
---|
277 | xfree(worka->phi); |
---|
278 | return; |
---|
279 | } |
---|
280 | |
---|
281 | /* eof */ |
---|