Scippy

SCIP

Solving Constraint Integer Programs

main_vrp.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 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 email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file
17  * @brief main file for VRP pricer example
18  * @author Andreas Bley
19  * @author Marc Pfetsch
20  *
21  * We want to solve the vehicle routing problem on a graph G = (V,E) with
22  * V = J cup {d}, where d is the depot and the distances are given by the
23  * length function l_e: E -> R_{<= 0}.
24  *
25  * Consider the MIP formulation
26  *
27  * min sum_{e in E} l_e y_e
28  * s.t. -y_e + sum_{t in T_k} a^t_e x_t <= 0, for all e in E
29  * sum_{t in T_k} a^t_j x_t == 1, for all j in J
30  * y(delta(j)) == 2, for all j in J
31  * y_e in {0,1,2}, for all e in E
32  * x_t in [0,1], for all t in T_k
33  *
34  * where T_k is the set of tours visiting at most k customers
35  * with repetitions of customers allowed and a^t_e (a^t_j) counts how often
36  * edge e (node j) is traversed in t in T_k.
37  */
38 
39 /* standard library includes */
40 #include <stdio.h>
41 #include <iostream>
42 #include <fstream>
43 #include <vector>
44 #include <string>
45 
46 /* scip includes */
47 #include "objscip/objscip.h"
48 #include "objscip/objscipdefplugins.h"
49 
50 /* user defined includes */
51 #include "pricer_vrp.h"
52 
53 
54 /* namespace usage */
55 using namespace std;
56 using namespace scip;
57 
58 
59 /** read VRP problem */
60 int read_problem(
61  const char* filename, /**< filename */
62  int& num_nodes, /**< number of nodes in instance */
63  int& capacity, /**< capacity in instance */
64  vector<int>& demand, /**< array of demands of instance */
65  vector<vector<int> >& distance /**< distances between nodes */
66  )
67 {
68  static const string DIMENSION = "DIMENSION";
69  static const string DEMAND_SECTION = "DEMAND_SECTION";
70  static const string DEPOT_SECTION = "DEPOT_SECTION";
71  static const string EDGE_WEIGHT_TYPE = "EDGE_WEIGHT_TYPE";
72  static const string EUC_2D = "EUC_2D";
73  static const string EXPLICIT = "EXPLICIT";
74  static const string LOWER_DIAG_ROW = "LOWER_DIAG_ROW";
75  static const string EDGE_WEIGHT_FORMAT = "EDGE_WEIGHT_FORMAT";
76  static const string EDGE_WEIGHT_SECTION = "EDGE_WEIGHT_SECTION";
77  static const string NODE_COORD_SECTION = "NODE_COORD_SECTION";
78  static const string CAPACITY = "CAPACITY";
79 
80  ifstream file(filename);
81 
82  if ( ! file )
83  {
84  cerr << "Cannot open file " << filename << endl;
85  return 1;
86  }
87 
88  string edge_weight_type = "";
89  string edge_weight_format = "";
90  vector<int> x;
91  vector<int> y;
92 
93  while ( file )
94  {
95  //--------------------
96  // Read keyword.
97  //--------------------
98  string key;
99  string dummy;
100  file >> key;
101 
102  if ( key == DIMENSION )
103  {
104  file >> dummy;
105  file >> num_nodes;
106 
107  demand.resize(num_nodes, 0);
108  distance.resize(num_nodes);
109  for (int i = 0; i < num_nodes; ++i)
110  distance[i].resize(i, 0);
111  }
112 
113  if ( key == CAPACITY )
114  {
115  file >> dummy;
116  file >> capacity;
117  }
118  else if ( key == EDGE_WEIGHT_TYPE )
119  {
120  file >> dummy;
121  file >> edge_weight_type;
122  if ( edge_weight_type != EUC_2D && edge_weight_type != EXPLICIT )
123  {
124  cerr << "Wrong " << EDGE_WEIGHT_TYPE << " " << edge_weight_type << endl;
125  return 1;
126  }
127  if ( edge_weight_type == EUC_2D )
128  {
129  x.resize(num_nodes, 0);
130  y.resize(num_nodes, 0);
131  }
132  }
133  else if ( key == EDGE_WEIGHT_FORMAT )
134  {
135  file >> dummy;
136  file >> edge_weight_format;
137  }
138  else if ( key == EDGE_WEIGHT_FORMAT + ":" )
139  {
140  file >> edge_weight_format;
141  }
142  else if ( key == EDGE_WEIGHT_SECTION )
143  {
144  if ( edge_weight_type != EXPLICIT || edge_weight_format != LOWER_DIAG_ROW )
145  {
146  cerr << "Error. Unsupported edge length type." << endl;
147  return 1;
148  }
149  for (int i = 0; i < num_nodes; ++i)
150  {
151  for (int j = 0; j < i; ++j)
152  {
153  int l;
154  file >> l;
155  distance[i][j] = l;
156  }
157  }
158  }
159  else if ( key == NODE_COORD_SECTION )
160  {
161  if ( edge_weight_type != EUC_2D )
162  {
163  cerr << "Error. Data file contains " << EDGE_WEIGHT_TYPE << " " << edge_weight_type << " and " << NODE_COORD_SECTION << endl;
164  return 1;
165  }
166  for (int i = 0; i < num_nodes; ++i)
167  {
168  int j, xi, yi;
169  file >> j;
170  file >> xi;
171  file >> yi;
172  if ( j != i+1 )
173  {
174  cerr << "Error reading " << NODE_COORD_SECTION << endl;
175  return 1;
176  }
177  x[i] = xi;
178  y[i] = yi;
179  }
180  for (int i = 0; i < num_nodes; ++i)
181  {
182  for (int j = 0; j < i; ++j)
183  {
184  int dx = x[i] - x[j];
185  int dy = y[i] - y[j];
186  distance[i][j] = int( sqrt((double)dx*dx + dy*dy) + 0.5 );
187  }
188  }
189  }
190  else if ( key == DEMAND_SECTION )
191  {
192  for (int i = 0; i < num_nodes; ++i)
193  {
194  int j, d;
195  file >> j;
196  file >> d;
197  if ( j != i+1 )
198  {
199  cerr << "Error reading " << DEMAND_SECTION << endl;
200  return 1;
201  }
202  demand[i] = d;
203  }
204  }
205  else if ( key == DEPOT_SECTION )
206  {
207  for (int i = 0; i != -1 ;)
208  {
209  file >> i;
210  if ( i != -1 && i != 1 )
211  {
212  cerr << "Error: This file specifies other depots than 1." << endl;
213  return 1;
214  }
215  }
216  }
217  else
218  {
219  getline(file, dummy);
220  }
221  }
222 
223  return 0;
224 }
225 
226 
227 
228 //------------------------------------------------------------
229 int main(int argc, char** argv)
230 {
231  SCIP* scip = NULL;
232 
233  cout << "Solving the vehicle routing problem using SCIP." << endl;
234  cout << "Implemented by Andreas Bley." << endl << endl;
235 
236  if ( argc != 2 && argc != 3 )
237  {
238  cerr << "Usage: vrp [-h] datafile" << endl;
239  cerr << "Options:" << endl;
240  cerr << " -h Uses hop limit instead of capacity limit for tours."<< endl;
241  return 1;
242  }
243 
244 
245  /**********************
246  * Setup problem data *
247  **********************/
248 
249  static const char* VRP_PRICER_NAME = "VRP_Pricer";
250 
251  vector<vector<int> > distance;
252  vector<int> demand;
253  int capacity;
254  int num_nodes;
255 
256  if ( read_problem(argv[argc-1], num_nodes, capacity, demand, distance) )
257  {
258  cerr << "Error reading data file " << argv[argc-1] << endl;
259  return 1;
260  }
261 
262  cout << "Number of nodes: " << num_nodes << endl;
263 
264  if ( argc == 3 )
265  {
266  if ( string("-h") != argv[1] )
267  {
268  cerr << "Unknow option " << argv[2] << endl;
269  return 1;
270  }
271 
272  int total_demand = 0;
273  for (int i = 1; i< num_nodes; ++i)
274  total_demand += demand[i];
275  capacity = (num_nodes - 1) * capacity / total_demand;
276  demand.assign(num_nodes, 1);
277  demand[0] = 0;
278  cout << "Max customers per tour: " << capacity << endl << endl;
279  }
280  else
281  cout << "Max demand per tour: " << capacity << endl << endl;
282 
283  /**************
284  * Setup SCIP *
285  **************/
286 
287  /* initialize SCIP environment */
288  SCIP_CALL( SCIPcreate(&scip) );
289 
290  /***********************
291  * Version information *
292  ***********************/
293 
294  SCIPprintVersion(scip, NULL);
295  SCIPinfoMessage(scip, NULL, "\n");
296 
297  /* include default plugins */
298  SCIP_CALL( SCIPincludeDefaultPlugins(scip) );
299 
300  /* set verbosity parameter */
301  SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 5) );
302  /* SCIP_CALL( SCIPsetBoolParam(scip, "display/lpinfo", TRUE) ); */
303 
304  /* create empty problem */
305  SCIP_CALL( SCIPcreateProb(scip, "VRP", 0, 0, 0, 0, 0, 0, 0) );
306 
307  /* add arc-routing variables */
308  char var_name[255];
309  vector< vector<SCIP_VAR*> > arc_var( num_nodes );
310  for (int i = 0; i < num_nodes; ++i)
311  {
312  arc_var[i].resize(i, (SCIP_VAR*) NULL);
313  for (int j = 0; j < i; ++j)
314  {
315  SCIP_VAR* var;
316  SCIPsnprintf(var_name, 255, "E%d_%d", i, j );
317 
318  SCIP_CALL( SCIPcreateVar(scip,
319  &var, // returns new index
320  var_name, // name
321  0.0, // lower bound
322  2.0, // upper bound
323  distance[i][j], // objective
324  SCIP_VARTYPE_INTEGER, // variable type
325  true, // initial
326  false, // forget the rest ...
327  0, 0, 0, 0, 0) );
328  SCIP_CALL( SCIPaddVar(scip, var) );
329  arc_var[i][j] = var;
330  }
331  }
332 
333  /* add arc-routing - tour constraints */
334  char con_name[255];
335  vector< vector<SCIP_CONS*> > arc_con( num_nodes );
336  for (int i = 0; i < num_nodes; ++i)
337  {
338  arc_con[i].resize(i, (SCIP_CONS*)NULL);
339  for (int j = 0; j < i; ++j)
340  {
341  SCIP_CONS* con;
342  SCIPsnprintf(con_name, 255, "A%d_%d", i, j);
343  SCIP_VAR* index = arc_var[i][j];
344  SCIP_Real coeff = -1;
345  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 1, &index, &coeff,
346  -SCIPinfinity(scip), /* lhs */
347  0.0, /* rhs */
348  true, /* initial */
349  false, /* separate */
350  true, /* enforce */
351  true, /* check */
352  true, /* propagate */
353  false, /* local */
354  true, /* modifiable */
355  false, /* dynamic */
356  false, /* removable */
357  false) ); /* stickingatnode */
358  SCIP_CALL( SCIPaddCons(scip, con) );
359  arc_con[i][j] = con;
360  }
361  }
362 
363  /* add arc-routing - degree constraints */
364  for (int i = 1; i < num_nodes; ++i)
365  {
366  SCIP_CONS* con;
367  SCIPsnprintf(con_name, 255, "D%d", i);
368  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 0, 0, 0,
369  2.0, /* lhs */
370  2.0, /* rhs */
371  true, /* initial */
372  false, /* separate */
373  true, /* enforce */
374  true, /* check */
375  true, /* propagate */
376  false, /* local */
377  false, /* modifiable */
378  false, /* dynamic */
379  false, /* removable */
380  false) ); /* stickingatnode */
381  SCIP_CALL( SCIPaddCons(scip, con) );
382  for (int j = 0; j < num_nodes; ++j)
383  {
384  if ( j != i )
385  {
386  SCIP_CALL( SCIPaddCoefLinear(scip, con, i > j ? arc_var[i][j] : arc_var[j][i], 1.0) );
387  }
388  }
389  }
390 
391  /* add set packing constraints (Node 0 is the depot) */
392  vector<SCIP_CONS*> part_con(num_nodes, (SCIP_CONS*)NULL);
393  for (int i = 1; i < num_nodes; ++i)
394  {
395  SCIP_CONS* con = NULL;
396  SCIPsnprintf(con_name, 255, "C%d", i);
397  SCIP_CALL( SCIPcreateConsLinear( scip, &con, con_name, 0, NULL, NULL,
398  1.0, /* lhs */
399  SCIPinfinity(scip), /* rhs */
400  true, /* initial */
401  false, /* separate */
402  true, /* enforce */
403  true, /* check */
404  true, /* propagate */
405  false, /* local */
406  true, /* modifiable */
407  false, /* dynamic */
408  false, /* removable */
409  false /* stickingatnode */ ) );
410  SCIP_CALL( SCIPaddCons(scip, con) );
411  part_con[i] = con;
412  }
413 
414  /* include VRP pricer */
415  ObjPricerVRP* vrp_pricer_ptr = new ObjPricerVRP(scip, VRP_PRICER_NAME, num_nodes, capacity, demand, distance,
416  arc_var, arc_con, part_con);
417 
418  SCIP_CALL( SCIPincludeObjPricer(scip, vrp_pricer_ptr, true) );
419 
420  /* activate pricer */
421  SCIP_CALL( SCIPactivatePricer(scip, SCIPfindPricer(scip, VRP_PRICER_NAME)) );
422 
423  // SCIP_CALL( SCIPwriteOrigProblem(scip, "vrp_init.lp", "lp", FALSE) );
424 
425 
426  /*************
427  * Solve *
428  *************/
429 
430  SCIP_CALL( SCIPsolve(scip) );
431 
432 
433  /**************
434  * Statistics *
435  *************/
436  SCIP_CALL( SCIPprintStatistics(scip, NULL) );
437 
438  SCIP_CALL( SCIPprintBestSol(scip, NULL, FALSE) );
439 
440 
441 
442  /********************
443  * Deinitialization *
444  ********************/
445 
446  SCIP_CALL( SCIPfree(&scip) );
447 
448  BMScheckEmptyMemory();
449 
450  return 0;
451 }
452