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