1 /*===========================================================================*/
2 /* */
3 /* This file is part of a demonstration application for use with the */
4 /* SYMPHONY Branch, Cut, and Price Library. This application is a solver for */
5 /* Capacitated Network Routing Problems. */
6 /* */
7 /* (c) Copyright 2000-2007 Ted Ralphs. All Rights Reserved. */
8 /* */
9 /* This application was developed by Ted Ralphs (ted@lehigh.edu) */
10 /* */
11 /* This software is licensed under the Eclipse Public License. Please see */
12 /* accompanying file for terms. */
13 /* */
14 /*===========================================================================*/
15
16 /* system include files */
17 #include <string.h>
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21
22 /* SYMPHONY include files */
23 #include "sym_macros.h"
24 #include "symphony.h"
25 #include "sym_master.h"
26
27 /* CNRP include files */
28 #include "cnrp_const.h"
29 #include "cnrp_types.h"
30 #include "cnrp_io.h"
31 #include "cnrp_master_functions.h"
32 #include "small_graph.h"
33
34 void cnrp_load_problem(sym_environment *env, cnrp_problem *cnrp);
35
36 int cnrp_test(sym_environment *env);
37
38 /*===========================================================================*\
39 * This is main() for the CNRP application
40 \*===========================================================================*/
41
main(int argc,char ** argv)42 int main(int argc, char **argv)
43 {
44 /* Allocate the data structure to store the problem-specific data */
45 cnrp_problem *cnrp = (cnrp_problem *) calloc(1, sizeof(cnrp_problem));
46 double ub;
47 char *pf;
48
49 /* Open the SYMPHONY environemt */
50 sym_environment *env = sym_open_environment();
51
52 /* Print version info */
53 sym_version();
54
55 /* Parse the command line arguments */
56 sym_parse_command_line(env, argc, argv);
57
58 /* Tell SYMPHONY where the user data is stored, so that this pointer can be
59 made available for user callbacks. */
60 sym_set_user_data(env, cnrp);
61
62 /* Read parameters from a file if there is one */
63 sym_get_str_param(env, "param_file", &pf);
64 cnrp_readparams(cnrp, pf, argc, argv);
65 cnrp->par.base_variable_selection = SOME_ARE_BASE;
66
67 sym_set_dbl_param(env, "granularity", .999999);
68
69 if(cnrp->par.test){
70
71 cnrp_test(env);
72
73 } else {
74
75 /* Read in the data */
76 cnrp_io(cnrp, cnrp->par.infile);
77
78 /* Read in a sparse graph representation of the network from a file */
79 if (cnrp->par.use_small_graph == LOAD_SMALL_GRAPH){
80 read_small_graph(cnrp);
81 }
82
83 if (!cnrp->numroutes && cnrp->par.prob_type == VRP){
84 printf("\nError: Number of trucks not specified or computed "
85 "for VRP\n\n");
86 exit(1);
87 }
88
89 if (cnrp->numroutes > 1){
90 printf("NUMBER OF TRUCKS: \t%i\n", cnrp->numroutes);
91 printf("TIGHTNESS: \t\t%.2f\n",
92 cnrp->demand[0]/(cnrp->capacity*(double)cnrp->numroutes));
93 }
94
95 /* Select the cheapest edges adjacent to each node for the base set */
96 if (cnrp->par.use_small_graph == SAVE_SMALL_GRAPH){
97 if (!cnrp->g) make_small_graph(cnrp, 0);
98 save_small_graph(cnrp);
99 }else if (!cnrp->g){
100 make_small_graph(cnrp, 0);
101 }
102
103 /* Check the a priori upper bound and set it if needed */
104 sym_get_primal_bound(env, &ub);
105
106 printf("ub is: %f", ub);
107 if (ub < SYM_INFINITY){
108 cnrp->cur_tour->cost = (int) ub;
109 }else{
110 cnrp->cur_tour->cost = MAXINT;
111 }
112 cnrp->cur_tour->numroutes = cnrp->numroutes;
113
114 if (cnrp->par.use_small_graph == LOAD_SMALL_GRAPH){
115 if (ub == SYM_INFINITY && cnrp->cur_tour->cost > 0){
116 ub = sym_set_primal_bound(env, cnrp->cur_tour->cost);
117 }
118 cnrp->numroutes = cnrp->cur_tour->numroutes;
119 }
120
121 #if 0
122 if(cnrp->par.prob_type == BPP)
123 sym_set_primal_bound(env, 1);
124 #endif
125
126 if (ub < SYM_INFINITY && !(cnrp->par.prob_type == BPP)){
127 printf("INITIAL UPPER BOUND: \t%i\n\n", (int)(ub));
128 }else if (!(cnrp->par.prob_type == BPP)){
129 printf("INITIAL UPPER BOUND: \tNone\n\n");
130 }else{
131 printf("\n\n");
132 }
133
134 /* Create the variable set and decide which variables are in the base set */
135 cnrp_create_variables(cnrp);
136
137 /* Construct the problem instance and load it into SYMPHONY */
138 cnrp_load_problem(env, cnrp);
139
140 /* Solve the problem */
141
142 #ifdef MULTI_CRITERIA
143 sym_mc_solve(env);
144 #else
145 sym_solve(env);
146 #endif
147 }
148
149 /* Close the SYMPHONY environment */
150 sym_close_environment(env);
151
152 return(0);
153 }
154
155 /*===========================================================================*\
156 * This is the function that creates the formulation
157 \*===========================================================================*/
158
cnrp_load_problem(sym_environment * env,cnrp_problem * cnrp)159 void cnrp_load_problem(sym_environment *env, cnrp_problem *cnrp)
160 {
161 int n, m, nz;
162 char *is_int, *sense;
163 int *matbeg, *matind;
164 double *matval, *obj, *obj2, *rhs, *rngval;
165 double *lb, *ub;
166 int *costs = cnrp->dist.cost;
167 int i, j;
168 char resize = FALSE;
169 int vertnum = cnrp->vertnum;
170 int total_edgenum = vertnum*(vertnum-1)/2;
171 char prob_type = cnrp->par.prob_type;
172 int v0, v1;
173 double flow_capacity;
174 int *edges;
175 int *vars, phase, varindex;
176 #ifdef DIRECTED_X_VARS
177 /*whether or not we will have the out-degree constraints*/
178 char od_const = (cnrp->par.prob_type == TSP || cnrp->par.prob_type == VRP ||
179 cnrp->par.prob_type == BPP);
180 char d_x_vars = TRUE;
181 #else
182 char od_const = FALSE;
183 char d_x_vars = FALSE;
184 #endif
185 int varnum = cnrp->basevarnum + cnrp->extravarnum;
186
187 #if 0
188 if (varnum != 4*total_edgenum){
189 printf("Error: Variables missing from problem description.\n");
190 printf("Exiting...\n");
191 exit(0);
192 }
193 #endif
194
195 #if defined(ADD_FLOW_VARS)
196 #ifdef DIRECTED_X_VARS
197 flow_capacity = cnrp->capacity;
198 #else
199 if (cnrp->par.prob_type == CSTP || cnrp->par.prob_type == CTP)
200 flow_capacity = cnrp->capacity;
201 else
202 flow_capacity = cnrp->capacity/2;
203 #endif
204 #endif
205
206 #ifdef MULTI_CRITERIA
207 cnrp->lp_par.tau = cnrp->lp_par.gamma = 1.0;
208 #endif
209
210 /* set up the inital LP data */
211
212 n = varnum;
213 m = cnrp->basecutnum;
214 edges = cnrp->edges;
215
216 /*Estimate the number of nonzeros*/
217 #ifdef ADD_CAP_CUTS
218 nz = 12*varnum;
219 #elif defined(ADD_FLOW_VARS)
220 nz = 8*varnum;
221 #else
222 nz = 3*varnum;
223 #endif
224 #ifdef ADD_X_CUTS
225 nz += 2*varnum;
226 #endif
227
228 /* Allocate the arrays. These are owned by SYMPHONY after returning. */
229 matbeg = (int *) malloc((n + 1) * ISIZE);
230 matind = (int *) malloc(nz * ISIZE);
231 matval = (double *) malloc(nz * DSIZE);
232 obj = (double *) calloc(n, DSIZE);
233 obj2 = (double *) calloc(n, DSIZE);
234 ub = (double *) malloc(n * DSIZE);
235 lb = (double *) calloc(n, DSIZE); /* zero lower bounds */
236 rhs = (double *) malloc(m * DSIZE);
237 sense = (char *) malloc(m * CSIZE);
238 rngval = (double *) calloc(m, DSIZE);
239 is_int = (char *) calloc(n, CSIZE);
240
241 for (i = 0, j = 0; i < varnum; i++){
242 if (i < total_edgenum){
243 is_int[i] = TRUE;
244 ub[i] = 1.0;
245 matbeg[i] = j;
246 obj[i] = (double) (cnrp->lp_par.gamma*costs[i]);
247 if (prob_type == CSTP || prob_type == CTP){
248 /*cardinality constraint*/
249 matind[j] = 0;
250 matval[j++] = 1.0;
251 }
252 /*in-degree constraint*/
253 matval[j] = 1.0;
254 matind[j++] = edges[2*i+1];
255 #ifdef DIRECTED_X_VARS
256 /*out-degree constraint*/
257 if (od_const){
258 matval[j] = 1.0;
259 matind[j++] = vertnum + edges[2*i];
260 }
261 #else
262 if (prob_type == VRP || prob_type == TSP ||
263 prob_type == BPP || edges[2*vars[varindex]]){
264 matval[j] = 1.0;
265 matind[j++] = edges[2*i];
266 }
267 #endif
268 #ifdef ADD_CAP_CUTS
269 v0 = edges[2*i];
270 matval[j] = -flow_capacity + (v0 ? cnrp->demand[v0] : 0);
271 matind[j++] = (2 + od_const)*vertnum - 1 + i;
272 #ifndef DIRECTED_X_VARS
273 matval[j] = -flow_capacity+cnrp->demand[edges[2*i+1]];
274 matind[j++] = 2*cnrp->vertnum - 1 + total_edgenum + i;
275 #endif
276 #endif
277 #ifdef ADD_X_CUTS
278 matval[j] = 1.0;
279 matind[j++]=(2+od_const)*vertnum-1+2*total_edgenum+i;
280 #endif
281 }
282 #ifdef DIRECTED_X_VARS
283 else if (i < 2*total_edgenum){
284 is_int[i] = TRUE;
285 ub[i] = 1.0;
286 matbeg[i] = j;
287 obj[i] = (double)(cnrp->lp_par.gamma*
288 costs[i-total_edgenum]);
289 if (prob_type == CSTP || prob_type == CTP){
290 /*cardinality constraint*/
291 matind[j] = 0;
292 matval[j++] = 1.0;
293 }
294 /*in-degree constraint*/
295 if (od_const || edges[2*(i - total_edgenum)]){
296 matval[j] = 1.0;
297 matind[j++] = edges[2*(i - total_edgenum)];
298 }
299 /*out-degree constraint*/
300 if (od_const){
301 matval[j] = 1.0;
302 matind[j++] = vertnum+edges[2*(i-total_edgenum)+1];
303 }
304 #ifdef ADD_CAP_CUTS
305 matval[j] = -flow_capacity +
306 cnrp->demand[edges[2*(i - total_edgenum) + 1]];
307 matind[j++] = (2 + od_const)*vertnum - 1 + i;
308 #endif
309 #ifdef ADD_X_CUTS
310 matval[j] = 1.0;
311 matind[j++] = (2 + od_const)* vertnum - 1 + 2 * total_edgenum +
312 i - total_edgenum;
313 #endif
314 }
315 #endif
316 else if (i < (2+d_x_vars)*total_edgenum){
317 is_int[i] = FALSE;
318 v0 = edges[2*(i-(1+d_x_vars)*total_edgenum)];
319 ub[i] = flow_capacity - (v0 ? cnrp->demand[v0] : 0);
320 matbeg[i] = j;
321 obj2[i] = (double)(cnrp->lp_par.tau*
322 costs[i-(1+d_x_vars)*total_edgenum]);
323 #ifdef ADD_CAP_CUTS
324 matval[j] = 1.0;
325 matval[j+1] = 1.0;
326 if (edges[2*(i-(1+d_x_vars)*total_edgenum)])
327 matval[j+2] = -1.0;
328 matind[j++] = (2 + od_const)*vertnum - 1 + i -
329 (1+d_x_vars)*total_edgenum;
330 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
331 (1+d_x_vars)*total_edgenum) + 1] - 1;
332 if (edges[2*(i - (1+d_x_vars)*total_edgenum)])
333 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
334 (1+d_x_vars)*total_edgenum)] - 1;
335 #else
336 matval[j] = 1.0;
337 if (edges[2*(i-(1+d_x_vars)*total_edgenum)])
338 matval[j+1] = -1.0;
339 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
340 (1+d_x_vars)*total_edgenum) + 1] - 1;
341 if (edges[2*(i - (1+d_x_vars)*total_edgenum)])
342 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
343 (1+d_x_vars)*total_edgenum)] - 1;
344 #endif
345 }
346 else{
347 is_int[i] = FALSE;
348 v1 = edges[2*(i - (2+d_x_vars)*total_edgenum) + 1];
349 ub[i] = flow_capacity - cnrp->demand[v1];
350 matbeg[i] = j;
351 obj2[i] = (double) (cnrp->lp_par.tau*costs[i-
352 (2+d_x_vars)*total_edgenum]);
353 #ifdef ADD_CAP_CUTS
354 matval[j] = 1.0;
355 matval[j+1] = -1.0;
356 if (edges[2*(i - (2+d_x_vars)*total_edgenum)])
357 matval[j+2] = 1.0;
358 matind[j++] = (2+od_const)*vertnum - 1 + i -
359 (1+d_x_vars)*total_edgenum;
360 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
361 (2+d_x_vars)*total_edgenum)+1] - 1;
362 if (edges[2*(i - (2+d_x_vars)*total_edgenum)])
363 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
364 (2+d_x_vars)*total_edgenum)] - 1;
365 #else
366 matval[j] = -1.0;
367 if (edges[2*(i - (2+d_x_vars)*total_edgenum)])
368 matval[j+1] = 1.0;
369 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
370 (2+d_x_vars)*total_edgenum)+1] - 1;
371 if (edges[2*(i - (2+d_x_vars)*total_edgenum)])
372 matind[j++] = (1+od_const)*vertnum + edges[2*(i -
373 (2+d_x_vars)*total_edgenum)] - 1;
374 #endif
375 }
376 }
377 matbeg[i] = j;
378
379 /* set the initial right hand side */
380 if (od_const){
381 /*degree constraints for the depot*/
382 #if 0
383 rhs[0] = cnrp->numroutes;
384 sense[0] = 'E';
385 rhs[vertnum] = cnrp->numroutes;
386 sense[vertnum] = 'E';
387 #else
388 rhs[0] = 1.0;
389 sense[0] = 'G';
390 rhs[vertnum] = 1.0;
391 sense[vertnum] = 'G';
392 #endif
393 }else if (prob_type == VRP || prob_type == TSP || prob_type == BPP){
394 (rhs[0]) = 2*cnrp->numroutes;
395 sense[0] = 'E';
396 }else{
397 /*cardinality constraint*/
398 rhs[0] = vertnum - 1;
399 sense[0] = 'E';
400 }
401 for (i = vertnum - 1; i > 0; i--){
402 switch (prob_type){
403 case VRP:
404 case TSP:
405 case BPP:
406 if (od_const){
407 rhs[i] = 1.0;
408 sense[i] = 'E';
409 rhs[i+vertnum] = 1.0;
410 sense[i+vertnum] = 'E';
411 }else{
412 rhs[i] = 2.0;
413 sense[i] = 'E';
414 }
415 break;
416 case CSTP:
417 case CTP:
418 rhs[i] = 1.0;
419 #ifdef DIRECTED_X_VARS
420 sense[i] = 'E';
421 #else
422 sense[i] = 'G';
423 #endif
424 break;
425 }
426 #ifdef ADD_FLOW_VARS
427 rhs[(1+od_const)*vertnum + i - 1] = cnrp->demand[i];
428 sense[(1+od_const)*vertnum + i - 1] = 'E';
429 #endif
430 }
431 #ifdef ADD_CAP_CUTS
432 for (i = (2+od_const)*vertnum - 1;
433 i < (2+od_const)*vertnum - 1 + 2*total_edgenum; i++){
434 rhs[i] = 0.0;
435 sense[i] = 'L';
436 }
437 #endif
438 #ifdef ADD_X_CUTS
439 for (i = (2+od_const)*vertnum-1+2*total_edgenum;
440 i < (2+od_const)*vertnum-1+3*total_edgenum; i++){
441 rhs[i] = 1;
442 sense[i] = 'L';
443 }
444 #endif
445
446 /* Load in the problem */
447 #ifdef MULTI_CRITERIA
448 sym_explicit_load_problem(env, n, m, matbeg, matind, matval, lb, ub, is_int,
449 obj, obj2, sense, rhs, NULL, TRUE);
450 #else
451 sym_explicit_load_problem(env, n, m, matbeg, matind, matval, lb, ub, is_int,
452 obj, obj2, sense, rhs, NULL, FALSE);
453 #endif
454 }
455
456 /*===========================================================================*\
457 \*===========================================================================*/
458
cnrp_test(sym_environment * env)459 int cnrp_test(sym_environment *env)
460 {
461
462 int termcode, i, file_num = 34;
463 char input_files[34][MAX_FILE_NAME_LENGTH +1] = {"A/A-n34-k5.vrp",
464 "A/A-n32-k5.vrp",
465 "A/A-n33-k5.vrp",
466 "E/E-n13-k4.vrp",
467 "E/E-n22-k4.vrp",
468 "E/E-n23-k3.vrp",
469 "E/E-n30-k3.vrp",
470 "E/E-n33-k4.vrp",
471 "V/att-n48-k4.vrp",
472 "E/E-n51-k5.vrp",
473 "A/A-n33-k6.vrp",
474 "A/A-n36-k5.vrp",
475 "A/A-n37-k5.vrp",
476 "A/A-n38-k5.vrp",
477 "A/A-n39-k5.vrp",
478 "A/A-n39-k6.vrp",
479 "A/A-n45-k6.vrp",
480 "A/A-n46-k7.vrp",
481 "B/B-n31-k5.vrp",
482 "B/B-n34-k5.vrp",
483 "B/B-n35-k5.vrp",
484 "B/B-n38-k6.vrp",
485 "B/B-n39-k5.vrp",
486 "B/B-n41-k6.vrp",
487 "B/B-n43-k6.vrp",
488 "B/B-n44-k7.vrp",
489 "B/B-n45-k5.vrp",
490 "B/B-n50-k7.vrp",
491 "B/B-n51-k7.vrp",
492 "B/B-n52-k7.vrp",
493 "B/B-n56-k7.vrp",
494 "B/B-n64-k9.vrp",
495 "A/A-n48-k7.vrp",
496 "A/A-n53-k7.vrp"};
497
498 double sol[34] = {403, 403, 403,403, 403, 403, 403, 403, 403, 403, 403,
499 403, 403, 403, 403, 403, 403, 403, 403, 403, 403,
500 403, 403, 403, 403, 403, 403, 403, 403, 403, 403,
501 403, 403, 403};
502
503 char *input_dir = (char*)malloc(CSIZE*(MAX_FILE_NAME_LENGTH+1));
504 char *infile = (char*)malloc(CSIZE*(MAX_FILE_NAME_LENGTH+1));
505 double *obj_val = (double *)calloc(DSIZE,file_num);
506 double tol = 1e-06, ub;
507 cnrp_problem *cnrp = (cnrp_problem *) env->user;
508
509 if (strcmp(cnrp->par.test_dir, "") == 0){
510 strcpy(input_dir, "../../../VRPLIB");
511 } else{
512 strcpy(input_dir, cnrp->par.test_dir);
513 }
514
515 sym_set_int_param(env, "verbosity", -10);
516
517 for(i = 0; i<file_num; i++){
518
519 strcpy(infile, "");
520 sprintf(infile, "%s%s%s", input_dir, "/", input_files[i]);
521
522 cnrp_io(cnrp, infile);
523
524 if (!cnrp->g){
525 make_small_graph(cnrp, 0);
526 }
527 sym_get_primal_bound(env, &ub);
528
529 if (ub < SYM_INFINITY){
530 cnrp->cur_tour->cost = (int) ub;
531 }else{
532 cnrp->cur_tour->cost = MAXINT;
533 }
534 cnrp->cur_tour->numroutes = cnrp->numroutes;
535
536 cnrp_create_variables(cnrp);
537
538 cnrp_load_problem(env, cnrp);
539
540 printf("Solving %s...\n", input_files[i]);
541
542 sym_solve(env);
543
544 sym_get_obj_val(env, &obj_val[i]);
545
546 if((obj_val[i] < sol[i] + tol) &&
547 (obj_val[i] > sol[i] - tol)){
548 printf("Success!\n");
549 } else {
550 printf("Failure!(%f, %f) \n", obj_val[i], sol[i]);
551 }
552
553 if(env->mip->n && i + 1 < file_num ){
554 free_master_u(env);
555 strcpy(env->par.infile, "");
556
557 env->mip = (MIPdesc *) calloc(1, sizeof(MIPdesc));
558 cnrp_problem *cnrp = (cnrp_problem *) calloc(1, sizeof(cnrp_problem));
559 env->user = (void *) cnrp;
560 cnrp_readparams(cnrp, NULL, 0, NULL);
561 cnrp->par.prob_type = cnrp->cg_par.prob_type =
562 cnrp->lp_par.prob_type = CTP;
563 cnrp->par.base_variable_selection = SOME_ARE_BASE;
564 }
565 }
566
567 FREE(input_dir);
568 FREE(infile);
569 FREE(obj_val);
570
571 return (0);
572 }
573