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-2016 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> >& dist /**< 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); /*lint !e732 !e747*/
108  dist.resize(num_nodes); /*lint !e732 !e747*/
109  for (int i = 0; i < num_nodes; ++i)
110  dist[i].resize(i, 0); /*lint !e732 !e747*/
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); /*lint !e732 !e747*/
130  y.resize(num_nodes, 0); /*lint !e732 !e747*/
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  dist[i][j] = l; /*lint !e732 !e747*/
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; /*lint !e732 !e747*/
178  y[i] = yi; /*lint !e732 !e747*/
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]; /*lint !e732 !e747 !e864*/
185  int dy = y[i] - y[j]; /*lint !e732 !e747 !e864*/
186  dist[i][j] = int( sqrt((double)dx*dx + dy*dy) + 0.5 ); /*lint !e732 !e747 !e790*/
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; /*lint !e732 !e747*/
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  (void) getline(file, dummy);
220  }
221  }
222 
223  return 0;
224 }
225 
226 
227 //------------------------------------------------------------
228 SCIP_RETCODE execmain(int argc, char** argv)
229 {
230  SCIP* scip = NULL;
231 
232  cout << "Solving the vehicle routing problem using SCIP." << endl;
233  cout << "Implemented by Andreas Bley." << endl << endl;
234 
235  if ( argc != 2 && argc != 3 )
236  {
237  cerr << "Usage: vrp [-h] datafile" << endl;
238  cerr << "Options:" << endl;
239  cerr << " -h Uses hop limit instead of capacity limit for tours."<< endl;
240  return SCIP_INVALIDDATA;
241  }
242 
243 
244  /**********************
245  * Setup problem data *
246  **********************/
247 
248  static const char* VRP_PRICER_NAME = "VRP_Pricer";
249 
250  vector<vector<int> > dist;
251  vector<int> demand;
252  int capacity;
253  int num_nodes;
254 
255  if ( read_problem(argv[argc-1], num_nodes, capacity, demand, dist) )
256  {
257  cerr << "Error reading data file " << argv[argc-1] << endl;
258  return SCIP_READERROR;
259  }
260 
261  cout << "Number of nodes: " << num_nodes << endl;
262 
263  if ( argc == 3 )
264  {
265  if ( string("-h") != argv[1] )
266  {
267  cerr << "Unknow option " << argv[2] << endl;
268  return SCIP_PARAMETERUNKNOWN;
269  }
270 
271  int total_demand = 0;
272  for (int i = 1; i< num_nodes; ++i)
273  total_demand += demand[i]; /*lint !e732 !e747*/
274 
275  if( total_demand == 0.0 )
276  {
277  cerr << "Total demand is zero!" << endl;
278  return SCIP_INVALIDDATA;
279  }
280 
281  capacity = (num_nodes - 1) * capacity / total_demand;
282  demand.assign(num_nodes, 1);
283  demand[0] = 0; /*lint !e747*/
284  cout << "Max customers per tour: " << capacity << endl << endl;
285  }
286  else
287  cout << "Max demand per tour: " << capacity << endl << endl;
288 
289  /**************
290  * Setup SCIP *
291  **************/
292 
293  /* initialize SCIP environment */
294  SCIP_CALL( SCIPcreate(&scip) );
295 
296  /***********************
297  * Version information *
298  ***********************/
299 
300  SCIPprintVersion(scip, NULL);
301  SCIPinfoMessage(scip, NULL, "\n");
302 
303  /* include default plugins */
304  SCIP_CALL( SCIPincludeDefaultPlugins(scip) );
305 
306  /* set verbosity parameter */
307  SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 5) );
308  /* SCIP_CALL( SCIPsetBoolParam(scip, "display/lpinfo", TRUE) ); */
309 
310  /* create empty problem */
311  SCIP_CALL( SCIPcreateProb(scip, "VRP", 0, 0, 0, 0, 0, 0, 0) );
312 
313  /* add arc-routing variables */
314  char var_name[255];
315  vector< vector<SCIP_VAR*> > arc_var( num_nodes ); /*lint !e732 !e747*/
316  for (int i = 0; i < num_nodes; ++i)
317  {
318  arc_var[i].resize(i, (SCIP_VAR*) NULL); /*lint !e732 !e747*/
319  for (int j = 0; j < i; ++j)
320  {
321  SCIP_VAR* var;
322  (void) SCIPsnprintf(var_name, 255, "E%d_%d", i, j );
323 
324  SCIP_CALL( SCIPcreateVar(scip,
325  &var, // returns new index
326  var_name, // name
327  0.0, // lower bound
328  2.0, // upper bound
329  dist[i][j], // objective
330  SCIP_VARTYPE_INTEGER, // variable type
331  true, // initial
332  false, // forget the rest ...
333  NULL, NULL, NULL, NULL, NULL) ); /*lint !e732 !e747*/
334  SCIP_CALL( SCIPaddVar(scip, var) );
335  arc_var[i][j] = var; /*lint !e732 !e747*/
336  }
337  }
338 
339  /* add arc-routing - tour constraints */
340  char con_name[255];
341  vector< vector<SCIP_CONS*> > arc_con( num_nodes ); /*lint !e732 !e747*/
342  for (int i = 0; i < num_nodes; ++i)
343  {
344  arc_con[i].resize(i, (SCIP_CONS*)NULL); /*lint !e732 !e747*/
345  for (int j = 0; j < i; ++j)
346  {
347  SCIP_CONS* con;
348  (void) SCIPsnprintf(con_name, 255, "A%d_%d", i, j);
349  SCIP_VAR* idx = arc_var[i][j]; /*lint !e732 !e747*/
350  SCIP_Real coeff = -1;
351  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 1, &idx, &coeff,
352  -SCIPinfinity(scip), /* lhs */
353  0.0, /* rhs */
354  true, /* initial */
355  false, /* separate */
356  true, /* enforce */
357  true, /* check */
358  true, /* propagate */
359  false, /* local */
360  true, /* modifiable */
361  false, /* dynamic */
362  false, /* removable */
363  false) ); /* stickingatnode */
364  SCIP_CALL( SCIPaddCons(scip, con) );
365  arc_con[i][j] = con; /*lint !e732 !e747*/
366  }
367  }
368 
369  /* add arc-routing - degree constraints */
370  for (int i = 1; i < num_nodes; ++i)
371  {
372  SCIP_CONS* con;
373  (void) SCIPsnprintf(con_name, 255, "D%d", i);
374  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 0, 0, 0,
375  2.0, /* lhs */
376  2.0, /* rhs */
377  true, /* initial */
378  false, /* separate */
379  true, /* enforce */
380  true, /* check */
381  true, /* propagate */
382  false, /* local */
383  false, /* modifiable */
384  false, /* dynamic */
385  false, /* removable */
386  false) ); /* stickingatnode */
387  SCIP_CALL( SCIPaddCons(scip, con) );
388  for (int j = 0; j < num_nodes; ++j)
389  {
390  if ( j != i )
391  {
392  SCIP_CALL( SCIPaddCoefLinear(scip, con, i > j ? arc_var[i][j] : arc_var[j][i], 1.0) ); /*lint !e732 !e747*/
393  }
394  }
395  }
396 
397  /* add set packing constraints (Node 0 is the depot) */
398  vector<SCIP_CONS*> part_con(num_nodes, (SCIP_CONS*)NULL); /*lint !e732 !e747*/
399  for (int i = 1; i < num_nodes; ++i)
400  {
401  SCIP_CONS* con = NULL;
402  (void) SCIPsnprintf(con_name, 255, "C%d", i);
403  SCIP_CALL( SCIPcreateConsLinear( scip, &con, con_name, 0, NULL, NULL,
404  1.0, /* lhs */
405  SCIPinfinity(scip), /* rhs */
406  true, /* initial */
407  false, /* separate */
408  true, /* enforce */
409  true, /* check */
410  true, /* propagate */
411  false, /* local */
412  true, /* modifiable */
413  false, /* dynamic */
414  false, /* removable */
415  false /* stickingatnode */ ) );
416  SCIP_CALL( SCIPaddCons(scip, con) );
417  part_con[i] = con; /*lint !e732 !e747*/
418  }
419 
420  /* include VRP pricer */
421  ObjPricerVRP* vrp_pricer_ptr = new ObjPricerVRP(scip, VRP_PRICER_NAME, num_nodes, capacity, demand, dist,
422  arc_var, arc_con, part_con);
423 
424  SCIP_CALL( SCIPincludeObjPricer(scip, vrp_pricer_ptr, true) );
425 
426  /* activate pricer */
427  SCIP_CALL( SCIPactivatePricer(scip, SCIPfindPricer(scip, VRP_PRICER_NAME)) );
428 
429  // SCIP_CALL( SCIPwriteOrigProblem(scip, "vrp_init.lp", "lp", FALSE) );
430 
431 
432  /*************
433  * Solve *
434  *************/
435 
436  SCIP_CALL( SCIPsolve(scip) );
437 
438 
439  /**************
440  * Statistics *
441  *************/
442  SCIP_CALL( SCIPprintStatistics(scip, NULL) );
443 
444  SCIP_CALL( SCIPprintBestSol(scip, NULL, FALSE) );
445 
446 
447 
448  /********************
449  * Deinitialization *
450  ********************/
451  SCIP_CALL( SCIPfree(&scip) );
452 
453  BMScheckEmptyMemory();
454 
455  return SCIP_OKAY;
456 }
457 
458 int main(int argc, char** argv)
459 {
460  return execmain(argc, argv) != SCIP_OKAY ? 1 : 0;
461 }
std
Definition: pqueue.h:28
read_problem
int read_problem(const char *filename, int &num_nodes, int &capacity, vector< int > &demand, vector< vector< int > > &dist)
Definition: main_vrp.cpp:60
pricer_vrp.h
VRP pricer plugin.
main
int main(int argc, char **argv)
Definition: main_vrp.cpp:458
execmain
SCIP_RETCODE execmain(int argc, char **argv)
Definition: main_vrp.cpp:228