1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file probdata_scflp.c
17 * @brief Problem data for Stochastic Capacitated Facility Location problem
18 * @author Stephen J. Maher
19 *
20 * This file handles the main problem data used in that project. For more details see \ref SCFLP_PROBLEMDATA page.
21 *
22 * @page SCFLP_SOLVEPROB Solving the deterministic equivalent
23 *
24 * The probdata_scflp.c is used to store the global problem data and build the monolithic MIP and decomposed problems.
25 * First, the structure of the problem data is describe. This is followed by a description of how to solve the problem
26 * directly using SCIP or using Benders' decomposition.
27 *
28 * @section SCFLP_PROBLEMDATA The global problem data
29 *
30 * The problem data is accessible in all plugins. The function SCIPgetProbData() returns the pointer to that structure.
31 * We use this data structure to store all the information of the SCFLP. Since this structure is not visible in the
32 * other plugins, we implemented setter and getter functions to access this data. The problem data structure
33 * SCIP_ProbData is shown below.
34 *
35 * \code
36 * ** @brief Problem data which is accessible in all places
37 * *
38 * * This problem data is used to store the input of the cap instance, all variables which are created, and all
39 * * constraints. In addition, the probdata stores the data structures for the decomposed problem. This permits the
40 * * use of Benders' decomposition to solve the stochastic program.
41 * *
42 * struct SCIP_ProbData
43 * {
44 * SCIP** subproblems; **< the Benders' decomposition subproblems * SCIP_VAR**
45 * SCIP_VAR** facilityvars; **< all variables representing facilities *
46 * SCIP_VAR*** subfacilityvars; **< duplicates of the facility variables in the subproblems *
47 * SCIP_VAR**** customervars; **< all variables representing the satisfaction of demand per scenario *
48 * SCIP_CONS*** capconss; **< capacity constraints per facility per scenario *
49 * SCIP_CONS*** demandconss; **< demand constraints per customer per scenario *
50 * SCIP_CONS* sufficientcap; **< ensuring sufficient capacity is provided to satisfy demand (relatively complete recourse) *
51 * SCIP_Real** costs; **< the transportation costs to a customer from a facility *
52 * SCIP_Real** demands; **< the customer demands per scenario *
53 * SCIP_Real* capacity; **< the capacity of each facility *
54 * SCIP_Real* fixedcost; **< the fixed cost of opening each facility *
55 * int ncustomers; **< the number of customers *
56 * int nfacilities; **< the number of facilities *
57 * int nscenarios; **< the number of scenarios *
58 * SCIP_Bool usebenders; **< whether Benders' decomposition is used *
59 * };
60 * \endcode
61 *
62 * The function SCIPprobdataCreate() manages the creation of the SCFLP instance in SCIP. There are two types of
63 * formulations that can be produced in this example. The first is the monolithic deterministic equivalent. The second
64 * is the reformulated problem that decomposes the stochastic problem by scenarios. This alternative formulations is
65 * solved using Benders' decomposition. Depending on the solution method, some members of SCIP_ProbData will be unused.
66 * For example, subproblems and subfacilityvars are only used when Benders' decomposition is applied to solve the SCFLP.
67 *
68 * The probdata_scflp.c also provide interface methods to the global problem data. A list of all interface methods can be
69 * found in probdata_scflp.h.
70 *
71 * @section SCFLP_DETEQUIV Directly solving the deterministic equivalent using SCIP
72 *
73 * Within probdata_scflp.c, both the monolithic determinstic equivalent or the decomposed problem can be built within
74 * SCIP. The monolithic deterministic equivalent involve a since SCIP instances that is solved directly as a MIP. The
75 * problem that is build in SCIP is given in \ref SCFLP_DETEQUIVMODEL.
76 *
77 * @section SCFLP_BENDERS Solving the SCFLP using Benders' decomposition
78 *
79 * The model that is used to build the decomposed problem is given in \ref SCFLP_BENDERSMODEL. In this example, the
80 * default Benders' decomposition plugin is used to employ the Benders' decomposition framework, see
81 * src/scip/benders_default.h. Before calling SCIPcreateBendersDefault() to invoke the Benders' decomposition framework,
82 * the SCIP instances for the master problem and the subproblems must be created.
83 *
84 * The SCIP instance for the master problem includes only the first stage variables (the facility variables \f$x_{i}\f$)
85 * and the first stage constraints. Note, the auxiliary variables are not added to the master problem by the user, nor
86 * are any Benders' decomposition cuts.
87 *
88 * For each subproblem \f$s\f$, the SCIP instance is formulated with the second stage variables (the customer variables
89 * \f$y^{s}_{ij}\f$) and the second stage constraints. Also, the first stage variables are created for each scenario.
90 * These variables are copies of the master variables from the master SCIP instances and must be created by calling
91 * SCIPcreateVarBasic() or SCIPcreateVar(). The master problem variable copies that are created in the subproblem SCIP
92 * instances must have an objective coefficient of 0.0. This is inline with the classical application of Benders'
93 * decomposition.
94 *
95 * IMPORTANT: the master variables that are created for the subproblem SCIP instances must have the same name as the
96 * corresponding master variables in the master problem SCIP instance. This is because the mapping between the master
97 * and subproblem variables relies on the variable names. This mapping is used for setting up the subproblems to
98 * evaluate solutions from the master problem and generating Benders' cuts.
99 *
100 * Once the master and subproblem SCIP instances are created, the Benders' decomposition is invoked by calling the
101 * interface function SCIPcreateBendersDefault(). The parameters for this function are a SCIP instance for the master
102 * problem, an array of SCIP instances for the subproblems and the number of subproblems.
103 *
104 * The Benders' decomposition framework involves the use of constraint handlers within SCIP, src/scip/cons_benders.h and
105 * src/scip/cons_benderslp.h. In order to solve the master problem by adding Benders' cuts, src/scip/cons_benders.h and
106 * src/scip/cons_benderslp.h must be activated. This is done by setting the parameter "constraints/benders/active" and
107 * "constraints/benderslp/active" to TRUE.
108 *
109 * NOTE: it is not necessary to activate src/scip/cons_benderslp.h. The purpose of this constraint handler is to
110 * generate Benders' decomposition cut from solutions to the LP relaxation in the root node. These solutions are
111 * fractional, since the enforcement priority of benderslp is higher than the integer constraint handler. The benderslp
112 * constraint handler allows the user to employ the multi-phase algorithm of McDaniel and Devine (1977).
113 *
114 * McDaniel D, Devine M. A modified Benders’ partitioning algorithm for mixed integer programming. Management Science
115 * 1977;24(2):312–9
116 */
117
118 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
119
120 #include <string.h>
121
122 #include "probdata_scflp.h"
123
124 #include "scip/scip.h"
125 #include "scip/scipdefplugins.h"
126
127 #define DEFAULT_SCALINGFACTOR 5000.0
128
129 /** @brief Problem data which is accessible in all places
130 *
131 * This problem data is used to store the input of the SCFLP, all variables which are created, and all constraints.
132 */
133 struct SCIP_ProbData
134 {
135 SCIP** subproblems; /**< the Benders' decomposition subproblems */
136 SCIP_VAR** facilityvars; /**< all variables representing facilities */
137 SCIP_VAR*** subfacilityvars; /**< duplicates of the facility variables in the subproblems */
138 SCIP_VAR**** customervars; /**< all variables representing the satisfaction of demand per scenario */
139 SCIP_CONS*** capconss; /**< capacity constraints per facility per scenario */
140 SCIP_CONS*** demandconss; /**< demand constraints per customer per scenario */
141 SCIP_CONS* sufficientcap; /**< ensuring sufficient capacity is provided to satisfy demand (relatively complete recourse) */
142 SCIP_Real** costs; /**< the transportation costs to a customer from a facility */
143 SCIP_Real** demands; /**< the customer demands per scenario */
144 SCIP_Real* capacity; /**< the capacity of each facility */
145 SCIP_Real* fixedcost; /**< the fixed cost of opening each facility */
146 int ncustomers; /**< the number of customers */
147 int nfacilities; /**< the number of facilities */
148
149 /* probdata parameters */
150 int nscenarios; /**< the number of scenarios */
151 SCIP_Bool usebenders; /**< whether Benders' decomposition is used */
152 SCIP_Bool quadcosts; /**< should the problem be formulated with quadratic costs */
153 };
154
155
156
157 /**@name Local methods
158 *
159 * @{
160 */
161
162 /** creates the original problem */
163 static
createOriginalproblem(SCIP * scip,SCIP_VAR ** facilityvars,SCIP_VAR **** customervars,SCIP_CONS *** capconss,SCIP_CONS *** demandconss,SCIP_CONS ** sufficientcap,SCIP_Real ** costs,SCIP_Real ** demands,SCIP_Real * capacity,SCIP_Real * fixedcost,int ncustomers,int nfacilities,int nscenarios,SCIP_Bool quadcosts)164 SCIP_RETCODE createOriginalproblem(
165 SCIP* scip, /**< SCIP data structure */
166 SCIP_VAR** facilityvars, /**< all variables representing facilities */
167 SCIP_VAR**** customervars, /**< all variables representing the satisfaction of demand */
168 SCIP_CONS*** capconss, /**< capacity constraints per facility */
169 SCIP_CONS*** demandconss, /**< demand constraints per customer */
170 SCIP_CONS** sufficientcap, /**< ensuring sufficient capacity is provided to satisfy demand */
171 SCIP_Real** costs, /**< the transportation costs from a facility to a customer */
172 SCIP_Real** demands, /**< the customer demands */
173 SCIP_Real* capacity, /**< the capacity of each facility */
174 SCIP_Real* fixedcost, /**< the fixed cost of opening a facility */
175 int ncustomers, /**< the number of customers */
176 int nfacilities, /**< the number of facilities */
177 int nscenarios, /**< the number of scenarios */
178 SCIP_Bool quadcosts /**< should the problem be formulated with quadratic costs */
179 )
180 {
181 SCIP_CONS* cons;
182 SCIP_VAR* var;
183 SCIP_VAR* sqrvar;
184 SCIP_Real maxdemand;
185 SCIP_Real coeff;
186 SCIP_Real custcoeff;
187 SCIP_Real one = 1.0;
188 SCIP_Real minusone = -1.0;
189 int i;
190 int j;
191 int k;
192 char name[SCIP_MAXSTRLEN];
193
194
195 assert(scip != NULL);
196
197 /* adding the sufficient capacity constraints */
198 maxdemand = 0;
199 for( i = 0; i < nscenarios; i++)
200 {
201 SCIP_Real sumdemand = 0;
202 for( j = 0; j < ncustomers; j++ )
203 sumdemand += demands[j][i];
204
205 if( sumdemand > maxdemand )
206 maxdemand = sumdemand;
207 }
208
209 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "sufficientcapacity");
210 SCIP_CALL( SCIPcreateConsBasicLinear(scip, sufficientcap, name, 0, NULL, NULL, maxdemand, SCIPinfinity(scip)) );
211
212 SCIP_CALL( SCIPaddCons(scip, (*sufficientcap)) );
213
214 /* adds the capacity constraints to the scenario */
215 for( i = 0; i < nfacilities; i++ )
216 {
217 for( j = 0; j < nscenarios; j++ )
218 {
219 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "capacity_%d_%d", i, j);
220 SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, -SCIPinfinity(scip), 0.0) );
221
222 SCIP_CALL( SCIPaddCons(scip, cons) );
223
224 capconss[i][j] = cons;
225 }
226 }
227
228 /* adds the demand constraints to the scenario */
229 for( i = 0; i < ncustomers; i++ )
230 {
231 for( j = 0; j < nscenarios; j++ )
232 {
233 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "demand_%d_%d", i, j);
234 SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, demands[i][j], SCIPinfinity(scip)) );
235
236 SCIP_CALL( SCIPaddCons(scip, cons) );
237
238 demandconss[i][j] = cons;
239 }
240 }
241
242 for( i = 0; i < nfacilities; i++ )
243 {
244 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "facility_%d", i);
245 SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, 1.0, fixedcost[i], SCIP_VARTYPE_BINARY) );
246
247 SCIP_CALL( SCIPaddVar(scip, var) );
248
249 /* storing the variable in the facility variable list */
250 facilityvars[i] = var;
251
252 /* adding the variable to the capacity constraints */
253 for( j = 0; j < nscenarios; j++ )
254 SCIP_CALL( SCIPaddCoefLinear(scip, capconss[i][j], var, -capacity[i]) );
255
256 /* adding the variable to the sufficient capacity constraints */
257 SCIP_CALL( SCIPaddCoefLinear(scip, (*sufficientcap), var, capacity[i]) );
258 }
259
260 /* adding the customer variables to the scenario */
261 for( i = 0; i < ncustomers; i++ )
262 {
263 for( j = 0; j < nfacilities; j++ )
264 {
265 for( k = 0; k < nscenarios; k++ )
266 {
267 custcoeff = costs[i][j]/(SCIP_Real)nscenarios;
268
269 /* if the quadratic costs are used, then the customer coefficient is zero and the quadratic costs are scaled
270 * by 10,000. */
271 if( quadcosts )
272 {
273 coeff = custcoeff / DEFAULT_SCALINGFACTOR;
274 custcoeff = 0.0;
275 }
276
277 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customer(%d,%d,%d)", i, j, k);
278 SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, SCIPinfinity(scip), custcoeff,
279 SCIP_VARTYPE_CONTINUOUS) );
280
281 SCIP_CALL( SCIPaddVar(scip, var) );
282
283 /* storing the customer variable in the list */
284 customervars[i][j][k] = var;
285
286 if( costs[i][j] > 0 )
287 {
288 /* adding the variable to the capacity constraints */
289 SCIP_CALL( SCIPaddCoefLinear(scip, capconss[j][k], customervars[i][j][k], 1.0) );
290
291 /* adding the variable to the demand constraints */
292 SCIP_CALL( SCIPaddCoefLinear(scip, demandconss[i][k], customervars[i][j][k], 1.0) );
293
294 /* if the quadratic costs are used, then variables representing the square of the customer supply
295 * must be added
296 */
297 if( quadcosts )
298 {
299 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customersqr(%d,%d,%d)", i, j, k);
300 SCIP_CALL( SCIPcreateVarBasic(scip, &sqrvar, name, 0.0, SCIPinfinity(scip), coeff, SCIP_VARTYPE_CONTINUOUS) );
301
302 SCIP_CALL( SCIPaddVar(scip, sqrvar) );
303
304 /* add constraint var^2 <= sqrvar */
305 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customersqrcons(%d,%d,%d)", i, j, k);
306 SCIP_CALL( SCIPcreateConsBasicQuadratic(scip, &cons, name, 1, &sqrvar, &minusone, 1, &var, &var,
307 &one, -SCIPinfinity(scip), 0.0) );
308
309 SCIP_CALL( SCIPaddCons(scip, cons) );
310
311 SCIP_CALL( SCIPreleaseCons(scip, &cons) );
312 SCIP_CALL( SCIPreleaseVar(scip, &sqrvar) );
313 }
314 }
315 }
316 }
317 }
318
319 return SCIP_OKAY;
320 }
321
322 /** creates the Benders' decomposition master problem */
323 static
createMasterproblem(SCIP * scip,SCIP_VAR ** facilityvars,SCIP_CONS ** sufficientcap,SCIP_Real * capacity,SCIP_Real * fixedcost,SCIP_Real ** demands,int ncustomers,int nfacilities,int nscenarios)324 SCIP_RETCODE createMasterproblem(
325 SCIP* scip, /**< SCIP data structure */
326 SCIP_VAR** facilityvars, /**< all variables representing facilities */
327 SCIP_CONS** sufficientcap, /**< ensuring sufficient capacity is provided to satisfy demand */
328 SCIP_Real* capacity, /**< the capacity of each facility */
329 SCIP_Real* fixedcost, /**< the fixed cost of opening a facility */
330 SCIP_Real** demands, /**< the customer demands */
331 int ncustomers, /**< the number of customers */
332 int nfacilities, /**< the number of facilities */
333 int nscenarios /**< the number of scenarios */
334 )
335 {
336 SCIP_VAR* var;
337 SCIP_Real maxdemand;
338 int i;
339 int j;
340 char name[SCIP_MAXSTRLEN];
341 assert(scip != NULL);
342
343 SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(scip), SCIPgetVerbLevel(scip), SCIP_VERBLEVEL_NORMAL,
344 "Creating the master problem\n============\n");
345
346 /* adding the sufficient capacity constraints */
347 maxdemand = 0;
348 for( i = 0; i < nscenarios; i++)
349 {
350 SCIP_Real sumdemand = 0;
351 for( j = 0; j < ncustomers; j++ )
352 sumdemand += demands[j][i];
353
354 if( sumdemand > maxdemand )
355 maxdemand = sumdemand;
356 }
357
358 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "sufficientcapacity");
359 SCIP_CALL( SCIPcreateConsBasicLinear(scip, sufficientcap, name, 0, NULL, NULL, maxdemand, SCIPinfinity(scip)) );
360
361 SCIP_CALL( SCIPaddCons(scip, (*sufficientcap)) );
362
363 /* adding the facility variables */
364 for( i = 0; i < nfacilities; i++ )
365 {
366 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "facility_%d", i);
367 SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, 1.0, fixedcost[i], SCIP_VARTYPE_BINARY) );
368
369 SCIP_CALL( SCIPaddVar(scip, var) );
370
371 /* storing the variable in the facility variable list */
372 facilityvars[i] = var;
373
374 /* adding the variable to the sufficient capacity constraints */
375 SCIP_CALL( SCIPaddCoefLinear(scip, (*sufficientcap), var, capacity[i]) );
376 }
377
378 SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(scip), SCIPgetVerbLevel(scip), SCIP_VERBLEVEL_NORMAL,
379 "master problem has %d binary variables and 1 constraint\n\n", nfacilities);
380
381 return SCIP_OKAY;
382 }
383
384 /** creates the scenario subproblems */
385 static
createSubproblems(SCIP * scip,SCIP ** subproblems,SCIP_VAR ** facilityvars,SCIP_VAR *** subfacilityvars,SCIP_VAR **** customervars,SCIP_CONS *** capconss,SCIP_CONS *** demandconss,SCIP_Real ** costs,SCIP_Real ** demands,SCIP_Real * capacity,SCIP_Real * fixedcost,int ncustomers,int nfacilities,int nscenarios,SCIP_Bool quadcosts)386 SCIP_RETCODE createSubproblems(
387 SCIP* scip, /**< SCIP data structure */
388 SCIP** subproblems, /**< the Benders' decomposition subproblems */
389 SCIP_VAR** facilityvars, /**< all variables representing facilities */
390 SCIP_VAR*** subfacilityvars, /**< the copies of the facility variables in the subproblems */
391 SCIP_VAR**** customervars, /**< all variables representing the satisfaction of demand */
392 SCIP_CONS*** capconss, /**< capacity constraints per facility */
393 SCIP_CONS*** demandconss, /**< demand constraints per customer */
394 SCIP_Real** costs, /**< the transportation costs from a facility to a customer */
395 SCIP_Real** demands, /**< the customer demands */
396 SCIP_Real* capacity, /**< the capacity of each facility */
397 SCIP_Real* fixedcost, /**< the fixed cost of opening a facility */
398 int ncustomers, /**< the number of customers */
399 int nfacilities, /**< the number of facilities */
400 int nscenarios, /**< the number of scenarios */
401 SCIP_Bool quadcosts /**< should the problem be formulated with quadratic costs */
402 )
403 {
404 SCIP_CONS* cons;
405 SCIP_VAR* var;
406 SCIP_VAR* sqrvar;
407 SCIP_Real coeff;
408 SCIP_Real custcoeff;
409 SCIP_Real one = 1.0;
410 SCIP_Real minusone = -1.0;
411 int i;
412 int j;
413 int k;
414 char name[SCIP_MAXSTRLEN];
415
416
417 assert(scip != NULL);
418
419 SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(scip), SCIPgetVerbLevel(scip), SCIP_VERBLEVEL_NORMAL,
420 "Creating the subproblems\n============\n");
421
422 /* adds the capacity constraints to the scenario */
423 for( i = 0; i < nfacilities; i++ )
424 {
425 for( j = 0; j < nscenarios; j++ )
426 {
427 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "capacity_%d_%d", i, j);
428 SCIP_CALL( SCIPcreateConsBasicLinear(subproblems[j], &cons, name, 0, NULL, NULL, -SCIPinfinity(subproblems[j]), 0.0) );
429
430 SCIP_CALL( SCIPaddCons(subproblems[j], cons) );
431
432 capconss[i][j] = cons;
433 }
434 }
435
436 /* adds the demand constraints to the scenario */
437 for( i = 0; i < ncustomers; i++ )
438 {
439 for( j = 0; j < nscenarios; j++ )
440 {
441 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "demand_%d_%d", i, j);
442 SCIP_CALL( SCIPcreateConsBasicLinear(subproblems[j], &cons, name, 0, NULL, NULL, demands[i][j], SCIPinfinity(subproblems[j])) );
443
444 SCIP_CALL( SCIPaddCons(subproblems[j], cons) );
445
446 demandconss[i][j] = cons;
447 }
448 }
449
450 for( i = 0; i < nfacilities; i++ )
451 {
452 for( j = 0; j < nscenarios; j++ )
453 {
454 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "facility_%d", i);
455 SCIP_CALL( SCIPcreateVarBasic(subproblems[j], &var, name, 0.0, 1.0, 0.0, SCIP_VARTYPE_CONTINUOUS) );
456
457 SCIP_CALL( SCIPaddVar(subproblems[j], var) );
458
459 /* storing the variable in the facility variable list */
460 subfacilityvars[i][j] = var;
461
462 /* adding the variable to the capacity constraints */
463 SCIP_CALL( SCIPaddCoefLinear(subproblems[j], capconss[i][j], subfacilityvars[i][j], -capacity[i]) );
464 }
465 }
466
467 /* adding the customer variables to the scenario */
468 for( i = 0; i < ncustomers; i++ )
469 {
470 for( j = 0; j < nfacilities; j++ )
471 {
472 for( k = 0; k < nscenarios; k++ )
473 {
474 custcoeff = costs[i][j]/(SCIP_Real)nscenarios;
475
476 /* if the quadratic costs are used, then the customer coefficient is zero and the quadratic costs are scaled
477 * by 10,000. */
478 if( quadcosts )
479 {
480 coeff = custcoeff / 10000.0;
481 custcoeff = 0.0;
482 }
483
484 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customer(%d,%d,%d)", i, j, k);
485 SCIP_CALL( SCIPcreateVarBasic(subproblems[k], &var, name, 0.0, SCIPinfinity(subproblems[k]), custcoeff,
486 SCIP_VARTYPE_CONTINUOUS) );
487
488 SCIP_CALL( SCIPaddVar(subproblems[k], var) );
489
490 /* storing the customer variable in the list */
491 customervars[i][j][k] = var;
492
493 if( costs[i][j] > 0 )
494 {
495 /* adding the variable to the capacity constraints */
496 SCIP_CALL( SCIPaddCoefLinear(subproblems[k], capconss[j][k], customervars[i][j][k], 1.0) );
497
498 /* adding the variable to the demand constraints */
499 SCIP_CALL( SCIPaddCoefLinear(subproblems[k], demandconss[i][k], customervars[i][j][k], 1.0) );
500
501
502 /* if the quadratic costs are used, then variables representing the square of the customer supply
503 * must be added
504 */
505 if( quadcosts )
506 {
507 coeff = costs[i][j]/(SCIP_Real)nscenarios / DEFAULT_SCALINGFACTOR;
508 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customersqr(%d,%d,%d)", i, j, k);
509 SCIP_CALL( SCIPcreateVarBasic(subproblems[k], &sqrvar, name, 0.0, SCIPinfinity(subproblems[k]), coeff, SCIP_VARTYPE_CONTINUOUS) );
510
511 SCIP_CALL( SCIPaddVar(subproblems[k], sqrvar) );
512
513 /* add constraint var^2 <= sqrvar */
514 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "customersqrcons(%d,%d,%d)", i, j, k);
515 SCIP_CALL( SCIPcreateConsBasicQuadratic(subproblems[k], &cons, name, 1, &sqrvar, &minusone, 1, &var, &var, &one, -SCIPinfinity(subproblems[k]), 0.0) );
516
517 SCIP_CALL( SCIPaddCons(subproblems[k], cons) );
518
519 SCIP_CALL( SCIPreleaseCons(subproblems[k], &cons) );
520 SCIP_CALL( SCIPreleaseVar(subproblems[k], &sqrvar) );
521 }
522 }
523 }
524 }
525 }
526
527 SCIPmessagePrintVerbInfo(SCIPgetMessagehdlr(scip), SCIPgetVerbLevel(scip), SCIP_VERBLEVEL_NORMAL,
528 "%d subproblems have been created.\neach subproblem has %d continuous variables and %d constraint\n\n",
529 nscenarios, ncustomers*nfacilities + nfacilities, nfacilities + ncustomers);
530
531 return SCIP_OKAY;
532 }
533
534
535 /** creates problem data */
536 static
probdataCreate(SCIP * scip,SCIP_PROBDATA ** probdata,SCIP ** subproblems,SCIP_VAR ** facilityvars,SCIP_VAR *** subfacilityvars,SCIP_VAR **** customervars,SCIP_CONS *** capconss,SCIP_CONS *** demandconss,SCIP_CONS * sufficientcap,SCIP_Real ** costs,SCIP_Real ** demands,SCIP_Real * capacity,SCIP_Real * fixedcost,int ncustomers,int nfacilities,int nscenarios,SCIP_Bool usebenders,SCIP_Bool quadcosts)537 SCIP_RETCODE probdataCreate(
538 SCIP* scip, /**< SCIP data structure */
539 SCIP_PROBDATA** probdata, /**< pointer to problem data */
540 SCIP** subproblems, /**< the Benders' decomposition subproblems */
541 SCIP_VAR** facilityvars, /**< all variables representing facilities */
542 SCIP_VAR*** subfacilityvars, /**< the copies of the facility variables in the subproblems */
543 SCIP_VAR**** customervars, /**< all variables representing the satisfaction of demand */
544 SCIP_CONS*** capconss, /**< capacity constraints per facility per scenario */
545 SCIP_CONS*** demandconss, /**< demand constraints per customer per scenario */
546 SCIP_CONS* sufficientcap, /**< ensuring sufficient capacity is provided to satisfy demand */
547 SCIP_Real** costs, /**< the transportation costs to a customer from a facility */
548 SCIP_Real** demands, /**< the customer demands per scenario */
549 SCIP_Real* capacity, /**< the capacity of each facility */
550 SCIP_Real* fixedcost, /**< the fixed cost of opening a facility */
551 int ncustomers, /**< the number of customers */
552 int nfacilities, /**< the number of facilities */
553 int nscenarios, /**< the number of scenarios */
554 SCIP_Bool usebenders, /**< whether Benders' decomposition is used */
555 SCIP_Bool quadcosts /**< should the problem be formulated with quadratic costs */
556 )
557 {
558 int i;
559 int j;
560
561 assert(scip != NULL);
562 assert(probdata != NULL);
563
564 /* allocate memory */
565 SCIP_CALL( SCIPallocBlockMemory(scip, probdata) );
566
567 /* copying the subproblem information */
568 if( usebenders )
569 {
570 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->subproblems, subproblems, nscenarios) );
571
572 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->subfacilityvars, nfacilities) );
573 for( i = 0; i < nfacilities; i++ )
574 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->subfacilityvars[i], subfacilityvars[i], nscenarios) );
575 }
576
577 /* copy variable arrays */
578 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->facilityvars, facilityvars, nfacilities) );
579 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->customervars, ncustomers) );
580 for( i = 0; i < ncustomers; i++ )
581 {
582 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->customervars[i], nfacilities) );
583 for( j = 0; j < nfacilities; j++ )
584 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->customervars[i][j], customervars[i][j],
585 nscenarios) );
586 }
587
588 /* duplicate the constraint arrays */
589 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->capconss, nfacilities) );
590 for( i = 0; i < nfacilities; i++ )
591 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->capconss[i], capconss[i], nscenarios) );
592
593 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->demandconss, ncustomers) );
594 for( i = 0; i < ncustomers; i++ )
595 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->demandconss[i], demandconss[i], nscenarios) );
596
597 /* duplicate the data arrays */
598 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->demands, ncustomers) );
599 for( i = 0; i < ncustomers; i++ )
600 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->demands[i], demands[i], nscenarios) );
601
602 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->costs, ncustomers) );
603 for( i = 0; i < ncustomers; i++ )
604 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->costs[i], costs[i], nfacilities) );
605
606 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->capacity, capacity, nfacilities) );
607 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->fixedcost, fixedcost, nfacilities) );
608
609 (*probdata)->sufficientcap = sufficientcap;
610 (*probdata)->ncustomers = ncustomers;
611 (*probdata)->nfacilities = nfacilities;
612 (*probdata)->nscenarios = nscenarios;
613 (*probdata)->usebenders = usebenders;
614 (*probdata)->quadcosts = quadcosts;
615
616 return SCIP_OKAY;
617 }
618
619 /** frees the memory of the given problem data */
620 static
probdataFree(SCIP * scip,SCIP_PROBDATA ** probdata)621 SCIP_RETCODE probdataFree(
622 SCIP* scip, /**< SCIP data structure */
623 SCIP_PROBDATA** probdata /**< pointer to problem data */
624 )
625 {
626 int i;
627 int j;
628 int k;
629
630 assert(scip != NULL);
631 assert(probdata != NULL);
632
633 #if 1
634 /* release all variables */
635 for( i = 0; i < (*probdata)->nfacilities; i++ )
636 SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->facilityvars[i]) );
637
638 for( i = 0; i < (*probdata)->nscenarios; i++ )
639 {
640 SCIP* varscip;
641 if( (*probdata)->usebenders )
642 varscip = (*probdata)->subproblems[i];
643 else
644 varscip = scip;
645
646 for( j = 0; j < (*probdata)->nfacilities; j++ )
647 {
648 for( k = 0; k < (*probdata)->ncustomers; k++ )
649 SCIP_CALL( SCIPreleaseVar(varscip, &(*probdata)->customervars[k][j][i]) );
650 }
651 }
652
653 /* release all constraints */
654 for( i = 0; i < (*probdata)->nscenarios; ++i )
655 {
656 SCIP* consscip;
657 if( (*probdata)->usebenders )
658 consscip = (*probdata)->subproblems[i];
659 else
660 consscip = scip;
661
662 for( j = 0; j < (*probdata)->ncustomers; j++ )
663 SCIP_CALL( SCIPreleaseCons(consscip, &(*probdata)->demandconss[j][i]) );
664 }
665
666 for( i = 0; i < (*probdata)->nscenarios; ++i )
667 {
668 SCIP* consscip;
669 if( (*probdata)->usebenders )
670 consscip = (*probdata)->subproblems[i];
671 else
672 consscip = scip;
673
674 for( j = 0; j < (*probdata)->nfacilities; ++j )
675 SCIP_CALL( SCIPreleaseCons(consscip, &(*probdata)->capconss[j][i]) );
676 }
677
678 SCIP_CALL( SCIPreleaseCons(scip, &(*probdata)->sufficientcap) );
679 #endif
680
681 /* free memory of arrays */
682 SCIPfreeBlockMemoryArray(scip, &(*probdata)->fixedcost, (*probdata)->nfacilities);
683 SCIPfreeBlockMemoryArray(scip, &(*probdata)->capacity, (*probdata)->nfacilities);
684
685 for( i = (*probdata)->ncustomers - 1; i >= 0; i-- )
686 SCIPfreeBlockMemoryArray(scip, &(*probdata)->costs[i], (*probdata)->nfacilities);
687 SCIPfreeBlockMemoryArray(scip, &(*probdata)->costs, (*probdata)->ncustomers);
688
689 for( i = (*probdata)->ncustomers - 1; i >= 0; i-- )
690 SCIPfreeBlockMemoryArray(scip, &(*probdata)->demands[i], (*probdata)->nscenarios);
691 SCIPfreeBlockMemoryArray(scip, &(*probdata)->demands, (*probdata)->ncustomers);
692
693 /* freeing the constraint memory arrays */
694 for( i = (*probdata)->ncustomers - 1; i >= 0; i-- )
695 SCIPfreeBlockMemoryArray(scip, &(*probdata)->demandconss[i], (*probdata)->nscenarios);
696 SCIPfreeBlockMemoryArray(scip, &(*probdata)->demandconss, (*probdata)->ncustomers);
697
698 for( i = (*probdata)->nfacilities - 1; i >= 0; i-- )
699 SCIPfreeBlockMemoryArray(scip, &(*probdata)->capconss[i], (*probdata)->nscenarios);
700 SCIPfreeBlockMemoryArray(scip, &(*probdata)->capconss, (*probdata)->nfacilities);
701
702 /* freeing the variable memory arrays */
703 for( i = (*probdata)->ncustomers - 1; i >= 0; i-- )
704 {
705 for( j = (*probdata)->nfacilities - 1; j >= 0; j-- )
706 SCIPfreeBlockMemoryArray(scip, &(*probdata)->customervars[i][j], (*probdata)->nscenarios);
707
708 SCIPfreeBlockMemoryArray(scip, &(*probdata)->customervars[i], (*probdata)->nfacilities);
709 }
710 SCIPfreeBlockMemoryArray(scip, &(*probdata)->customervars, (*probdata)->ncustomers);
711
712 SCIPfreeBlockMemoryArray(scip, &(*probdata)->facilityvars, (*probdata)->nfacilities);
713
714 /* freeing the subproblem information */
715 if( (*probdata)->usebenders )
716 {
717 /* freeing the sub facility variables */
718 for( i = 0; i < (*probdata)->nscenarios; i++ )
719 {
720 for( j = 0; j < (*probdata)->nfacilities; j++ )
721 SCIP_CALL( SCIPreleaseVar((*probdata)->subproblems[i], &(*probdata)->subfacilityvars[j][i]) );
722 }
723
724 for( i = (*probdata)->nfacilities - 1; i >= 0; i-- )
725 SCIPfreeBlockMemoryArray(scip, &(*probdata)->subfacilityvars[i], (*probdata)->nscenarios);
726
727 SCIPfreeBlockMemoryArray(scip, &(*probdata)->subfacilityvars, (*probdata)->nfacilities);
728
729 for( i = (*probdata)->nscenarios - 1; i >= 0 ; i-- )
730 SCIP_CALL( SCIPfree(&(*probdata)->subproblems[i]) );
731 SCIPfreeBlockMemoryArray(scip, &(*probdata)->subproblems, (*probdata)->nscenarios);
732 }
733
734 /* free probdata */
735 SCIPfreeBlockMemory(scip, probdata);
736
737 return SCIP_OKAY;
738 }
739
740 /**@} */
741
742 /**@name Callback methods of problem data
743 *
744 * @{
745 */
746
747 /** frees user data of original problem (called when the original problem is freed) */
748 static
SCIP_DECL_PROBDELORIG(probdelorigScflp)749 SCIP_DECL_PROBDELORIG(probdelorigScflp)
750 {
751 assert(scip != NULL);
752 assert(probdata != NULL);
753
754 SCIPdebugMsg(scip, "free original problem data\n");
755
756 SCIP_CALL( probdataFree(scip, probdata) );
757
758 return SCIP_OKAY;
759 }
760
761 /** creates user data of transformed problem by transforming the original user problem data
762 * (called after problem was transformed) */
763 static
SCIP_DECL_PROBTRANS(probtransScflp)764 SCIP_DECL_PROBTRANS(probtransScflp)
765 {
766 SCIPdebugMsg(scip, "transforming problem data\n");
767
768 return SCIP_OKAY;
769 }
770
771 /** frees user data of transformed problem (called when the transformed problem is freed) */
772 static
SCIP_DECL_PROBDELTRANS(probdeltransScflp)773 SCIP_DECL_PROBDELTRANS(probdeltransScflp)
774 {
775 SCIPdebugMsg(scip, "free transformed problem data\n");
776
777 return SCIP_OKAY;
778 }
779
780 /**@} */
781
782
783 /**@name Interface methods
784 *
785 * @{
786 */
787
788 /** sets up the problem data */
SCIPprobdataCreate(SCIP * scip,const char * probname,SCIP_Real ** costs,SCIP_Real ** demands,SCIP_Real * capacity,SCIP_Real * fixedcost,int ncustomers,int nfacilities,int nscenarios,SCIP_Bool usebenders,SCIP_Bool quadcosts)789 SCIP_RETCODE SCIPprobdataCreate(
790 SCIP* scip, /**< SCIP data structure */
791 const char* probname, /**< problem name */
792 SCIP_Real** costs, /**< the transportation costs from a facility to a customer */
793 SCIP_Real** demands, /**< the customer demands */
794 SCIP_Real* capacity, /**< the capacity of each facility */
795 SCIP_Real* fixedcost, /**< the fixed cost of opening a facility */
796 int ncustomers, /**< the number of customers */
797 int nfacilities, /**< the number of facilities */
798 int nscenarios, /**< the number of Benders' decomposition scenarios */
799 SCIP_Bool usebenders, /**< whether Benders' decomposition is used */
800 SCIP_Bool quadcosts /**< should the problem be formulated with quadratic costs */
801 )
802 {
803 SCIP** subproblems;
804 SCIP_PROBDATA* probdata;
805 SCIP_CONS*** demandconss;
806 SCIP_CONS*** capconss;
807 SCIP_CONS* sufficientcap;
808 SCIP_VAR** facilityvars;
809 SCIP_VAR*** subfacilityvars;
810 SCIP_VAR**** customervars;
811 int i;
812 int j;
813
814 assert(scip != NULL);
815
816 /* create problem in SCIP and add non-NULL callbacks via setter functions */
817 SCIP_CALL( SCIPcreateProbBasic(scip, probname) );
818
819 SCIP_CALL( SCIPsetProbDelorig(scip, probdelorigScflp) );
820 SCIP_CALL( SCIPsetProbTrans(scip, probtransScflp) );
821 SCIP_CALL( SCIPsetProbDeltrans(scip, probdeltransScflp) );
822
823 /* set objective sense */
824 SCIP_CALL( SCIPsetObjsense(scip, SCIP_OBJSENSE_MINIMIZE) );
825
826 SCIP_CALL( SCIPallocBufferArray(scip, &demandconss, ncustomers) );
827 for( i = 0; i < ncustomers; i++ )
828 SCIP_CALL( SCIPallocBufferArray(scip, &demandconss[i], nscenarios) );
829 SCIP_CALL( SCIPallocBufferArray(scip, &capconss, nfacilities) );
830 for( i = 0; i < nfacilities; i++ )
831 SCIP_CALL( SCIPallocBufferArray(scip, &capconss[i], nscenarios) );
832
833 SCIP_CALL( SCIPallocBufferArray(scip, &facilityvars, nfacilities) );
834 SCIP_CALL( SCIPallocBufferArray(scip, &customervars, ncustomers) );
835 for( i = 0; i < ncustomers; i++ )
836 {
837 SCIP_CALL( SCIPallocBufferArray(scip, &customervars[i], nfacilities) );
838 for( j = 0; j < nfacilities; j++ )
839 SCIP_CALL( SCIPallocBufferArray(scip, &customervars[i][j], nscenarios) );
840 }
841
842 sufficientcap = NULL;
843
844 subproblems = NULL;
845 subfacilityvars = NULL;
846
847 /* if quadratic costs are used, then the costs are scaled by 10,000. The user must be informed about this scaling */
848 if( quadcosts )
849 {
850 SCIPinfoMessage(scip, NULL, "The problem will be formulated with quadratic costs. "
851 "The input costs will be scaled by %g\n\n", DEFAULT_SCALINGFACTOR);
852 }
853
854 if( usebenders )
855 {
856 char subprobname[SCIP_MAXSTRLEN];
857
858 /* allocting the memory for the subproblem specific information */
859 SCIP_CALL( SCIPallocBufferArray(scip, &subproblems, nscenarios) );
860 SCIP_CALL( SCIPallocBufferArray(scip, &subfacilityvars, nfacilities) );
861 for( i = 0; i < nfacilities; i++ )
862 SCIP_CALL( SCIPallocBufferArray(scip, &subfacilityvars[i], nscenarios) );
863
864 /* creating the subproblems */
865 for( i = 0; i < nscenarios; i++ )
866 {
867 SCIP_CALL( SCIPcreate(&subproblems[i]) );
868
869 /* include default SCIP plugins */
870 SCIP_CALL( SCIPincludeDefaultPlugins(subproblems[i]) );
871
872 (void) SCIPsnprintf(subprobname, SCIP_MAXSTRLEN, "sub_%s_%d", probname, i);
873 SCIP_CALL( SCIPcreateProbBasic(subproblems[i], subprobname) );
874 }
875
876 /* creating the master problem */
877 SCIP_CALL( createMasterproblem(scip, facilityvars, &sufficientcap, capacity, fixedcost, demands, ncustomers,
878 nfacilities, nscenarios) );
879 SCIP_CALL( createSubproblems(scip, subproblems, facilityvars, subfacilityvars, customervars, capconss,
880 demandconss, costs, demands, capacity, fixedcost, ncustomers, nfacilities, nscenarios, quadcosts) );
881
882 /* including the Benders' decomposition plugin */
883 SCIP_CALL( SCIPcreateBendersDefault(scip, subproblems, nscenarios) );
884
885 /* activating the Benders' decomposition constraint handlers */
886 SCIP_CALL( SCIPsetBoolParam(scip, "constraints/benders/active", TRUE) );
887 SCIP_CALL( SCIPsetBoolParam(scip, "constraints/benderslp/active", TRUE) );
888
889 SCIP_CALL( SCIPsetIntParam(scip, "constraints/benders/maxprerounds", 1) );
890 SCIP_CALL( SCIPsetIntParam(scip, "presolving/maxrounds", 1) );
891 }
892 else
893 {
894 /* creating the original problem */
895 SCIP_CALL( createOriginalproblem(scip, facilityvars, customervars, capconss, demandconss, &sufficientcap, costs,
896 demands, capacity, fixedcost, ncustomers, nfacilities, nscenarios, quadcosts) );
897 }
898
899 /* create problem data */
900 SCIP_CALL( probdataCreate(scip, &probdata, subproblems, facilityvars, subfacilityvars, customervars, capconss,
901 demandconss, sufficientcap, costs, demands, capacity, fixedcost, ncustomers, nfacilities, nscenarios,
902 usebenders, quadcosts) );
903
904 /* set user problem data */
905 SCIP_CALL( SCIPsetProbData(scip, probdata) );
906
907 /* free local buffer arrays */
908 if( usebenders )
909 {
910 SCIPfreeBufferArray(scip, &subproblems);
911
912 for( i = nfacilities - 1; i >= 0; i-- )
913 SCIPfreeBufferArray(scip, &subfacilityvars[i]);
914 SCIPfreeBufferArray(scip, &subfacilityvars);
915 }
916
917 for( i = ncustomers - 1; i >= 0; i-- )
918 {
919 for( j = nfacilities - 1; j >= 0; j-- )
920 SCIPfreeBufferArray(scip, &customervars[i][j]);
921 SCIPfreeBufferArray(scip, &customervars[i]);
922 }
923 SCIPfreeBufferArray(scip, &customervars);
924 SCIPfreeBufferArray(scip, &facilityvars);
925
926 for( i = nfacilities - 1; i >= 0; i-- )
927 SCIPfreeBufferArray(scip, &capconss[i]);
928 SCIPfreeBufferArray(scip, &capconss);
929
930 for( i = ncustomers - 1; i >= 0; i-- )
931 SCIPfreeBufferArray(scip, &demandconss[i]);
932 SCIPfreeBufferArray(scip, &demandconss);
933
934 return SCIP_OKAY;
935 }
936
937 /** returns the number of facilities */
SCIPprobdataGetNFacilities(SCIP_PROBDATA * probdata)938 int SCIPprobdataGetNFacilities(
939 SCIP_PROBDATA* probdata /**< problem data */
940 )
941 {
942 assert(probdata != NULL);
943
944 return probdata->nfacilities;
945 }
946
947 /** returns the number of customers */
SCIPprobdataGetNCustomers(SCIP_PROBDATA * probdata)948 int SCIPprobdataGetNCustomers(
949 SCIP_PROBDATA* probdata /**< problem data */
950 )
951 {
952 assert(probdata != NULL);
953
954 return probdata->ncustomers;
955 }
956
957 /** returns the facility variables */
SCIPprobdataGetFacilityVars(SCIP_PROBDATA * probdata)958 SCIP_VAR** SCIPprobdataGetFacilityVars(
959 SCIP_PROBDATA* probdata /**< problem data */
960 )
961 {
962 assert(probdata != NULL);
963
964 return probdata->facilityvars;
965 }
966
967 /**@} */
968