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-2020 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
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"
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 static
62  const char* filename, /**< filename */
63  int& num_nodes, /**< number of nodes in instance */
64  int& capacity, /**< capacity in instance */
65  vector<int>& demand, /**< array of demands of instance */
66  vector<vector<int> >& dist /**< distances between nodes */
67  )
68 {
69  static const string DIMENSION = "DIMENSION";
70  static const string DEMAND_SECTION = "DEMAND_SECTION";
71  static const string DEPOT_SECTION = "DEPOT_SECTION";
72  static const string EDGE_WEIGHT_TYPE = "EDGE_WEIGHT_TYPE";
73  static const string EUC_2D = "EUC_2D";
74  static const string EXPLICIT = "EXPLICIT";
75  static const string LOWER_DIAG_ROW = "LOWER_DIAG_ROW";
76  static const string EDGE_WEIGHT_FORMAT = "EDGE_WEIGHT_FORMAT";
77  static const string EDGE_WEIGHT_SECTION = "EDGE_WEIGHT_SECTION";
78  static const string NODE_COORD_SECTION = "NODE_COORD_SECTION";
79  static const string CAPACITY = "CAPACITY";
80 
81  ifstream file(filename);
82 
83  if ( ! file )
84  {
85  cerr << "Cannot open file " << filename << endl;
86  return 1;
87  }
88 
89  string edge_weight_type = "";
90  string edge_weight_format = "";
91  vector<int> x;
92  vector<int> y;
93 
94  while ( file )
95  {
96  //--------------------
97  // Read keyword.
98  //--------------------
99  string key;
100  string dummy;
101  file >> key;
102 
103  if ( key == DIMENSION )
104  {
105  file >> dummy;
106  file >> num_nodes;
107 
108  demand.resize(num_nodes, 0); /*lint !e732 !e747*/
109  dist.resize(num_nodes); /*lint !e732 !e747*/
110  for (int i = 0; i < num_nodes; ++i)
111  dist[i].resize(i, 0); /*lint !e732 !e747*/
112  }
113 
114  if ( key == CAPACITY )
115  {
116  file >> dummy;
117  file >> capacity;
118  }
119  else if ( key == EDGE_WEIGHT_TYPE )
120  {
121  file >> dummy;
122  file >> edge_weight_type;
123  if ( edge_weight_type != EUC_2D && edge_weight_type != EXPLICIT )
124  {
125  cerr << "Wrong " << EDGE_WEIGHT_TYPE << " " << edge_weight_type << endl;
126  return 1;
127  }
128  if ( edge_weight_type == EUC_2D )
129  {
130  x.resize(num_nodes, 0); /*lint !e732 !e747*/
131  y.resize(num_nodes, 0); /*lint !e732 !e747*/
132  }
133  }
134  else if ( key == EDGE_WEIGHT_FORMAT )
135  {
136  file >> dummy;
137  file >> edge_weight_format;
138  }
139  else if ( key == EDGE_WEIGHT_FORMAT + ":" )
140  {
141  file >> edge_weight_format;
142  }
143  else if ( key == EDGE_WEIGHT_SECTION )
144  {
145  if ( edge_weight_type != EXPLICIT || edge_weight_format != LOWER_DIAG_ROW )
146  {
147  cerr << "Error. Unsupported edge length type." << endl;
148  return 1;
149  }
150  for (int i = 0; i < num_nodes; ++i)
151  {
152  for (int j = 0; j < i; ++j)
153  {
154  int l;
155  file >> l;
156  dist[i][j] = l; /*lint !e732 !e747*/
157  }
158  }
159  }
160  else if ( key == NODE_COORD_SECTION )
161  {
162  if ( edge_weight_type != EUC_2D )
163  {
164  cerr << "Error. Data file contains " << EDGE_WEIGHT_TYPE << " " << edge_weight_type << " and " << NODE_COORD_SECTION << endl;
165  return 1;
166  }
167  for (int i = 0; i < num_nodes; ++i)
168  {
169  int j, xi, yi;
170  file >> j;
171  file >> xi;
172  file >> yi;
173  if ( j != i+1 )
174  {
175  cerr << "Error reading " << NODE_COORD_SECTION << endl;
176  return 1;
177  }
178  x[i] = xi; /*lint !e732 !e747*/
179  y[i] = yi; /*lint !e732 !e747*/
180  }
181  for (int i = 0; i < num_nodes; ++i)
182  {
183  for (int j = 0; j < i; ++j)
184  {
185  int dx = x[i] - x[j]; /*lint !e732 !e747 !e864*/
186  int dy = y[i] - y[j]; /*lint !e732 !e747 !e864*/
187  dist[i][j] = int( sqrt((double)dx*dx + dy*dy) + 0.5 ); /*lint !e732 !e747 !e790*/
188  }
189  }
190  }
191  else if ( key == DEMAND_SECTION )
192  {
193  for (int i = 0; i < num_nodes; ++i)
194  {
195  int j, d;
196  file >> j;
197  file >> d;
198  if ( j != i+1 )
199  {
200  cerr << "Error reading " << DEMAND_SECTION << endl;
201  return 1;
202  }
203  demand[i] = d; /*lint !e732 !e747*/
204  }
205  }
206  else if ( key == DEPOT_SECTION )
207  {
208  for (int i = 0; i != -1 ;)
209  {
210  file >> i;
211  if ( i != -1 && i != 1 )
212  {
213  cerr << "Error: This file specifies other depots than 1." << endl;
214  return 1;
215  }
216  }
217  }
218  else
219  {
220  (void) getline(file, dummy);
221  }
222  }
223 
224  return 0;
225 }
226 
227 
228 //------------------------------------------------------------
229 static
230 SCIP_RETCODE execmain(int argc, char** argv)
231 {
232  SCIP* scip = NULL;
233 
234  cout << "Solving the vehicle routing problem using SCIP." << endl;
235  cout << "Implemented by Andreas Bley." << endl << endl;
236 
237  if ( argc != 2 && argc != 3 )
238  {
239  cerr << "Usage: vrp [-h] datafile" << endl;
240  cerr << "Options:" << endl;
241  cerr << " -h Uses hop limit instead of capacity limit for tours."<< endl;
242  return SCIP_INVALIDDATA;
243  }
244 
245 
246  /**********************
247  * Setup problem data *
248  **********************/
249 
250  static const char* VRP_PRICER_NAME = "VRP_Pricer";
251 
252  vector<vector<int> > dist;
253  vector<int> demand;
254  int capacity;
255  int num_nodes;
256 
257  if ( read_problem(argv[argc-1], num_nodes, capacity, demand, dist) )
258  {
259  cerr << "Error reading data file " << argv[argc-1] << endl;
260  return SCIP_READERROR;
261  }
262 
263  cout << "Number of nodes: " << num_nodes << endl;
264 
265  if ( argc == 3 )
266  {
267  if ( string("-h") != argv[1] )
268  {
269  cerr << "Unknow option " << argv[2] << endl;
270  return SCIP_PARAMETERUNKNOWN;
271  }
272 
273  int total_demand = 0;
274  for (int i = 1; i< num_nodes; ++i)
275  total_demand += demand[i]; /*lint !e732 !e747*/
276 
277  if( total_demand == 0.0 )
278  {
279  cerr << "Total demand is zero!" << endl;
280  return SCIP_INVALIDDATA;
281  }
282 
283  capacity = (num_nodes - 1) * capacity / total_demand;
284  demand.assign(num_nodes, 1);
285  demand[0] = 0; /*lint !e747*/
286  cout << "Max customers per tour: " << capacity << endl << endl;
287  }
288  else
289  cout << "Max demand per tour: " << capacity << endl << endl;
290 
291  /**************
292  * Setup SCIP *
293  **************/
294 
295  /* initialize SCIP environment */
296  SCIP_CALL( SCIPcreate(&scip) );
297 
298  /***********************
299  * Version information *
300  ***********************/
301 
302  SCIPprintVersion(scip, NULL);
303  SCIPinfoMessage(scip, NULL, "\n");
304 
305  /* include default plugins */
307 
308  /* set verbosity parameter */
309  SCIP_CALL( SCIPsetIntParam(scip, "display/verblevel", 5) );
310  /* SCIP_CALL( SCIPsetBoolParam(scip, "display/lpinfo", TRUE) ); */
311 
312  /* create empty problem */
313  SCIP_CALL( SCIPcreateProb(scip, "VRP", 0, 0, 0, 0, 0, 0, 0) );
314 
315  /* add arc-routing variables */
316  char var_name[255];
317  vector< vector<SCIP_VAR*> > arc_var( num_nodes ); /*lint !e732 !e747*/
318  for (int i = 0; i < num_nodes; ++i)
319  {
320  arc_var[i].resize(i, (SCIP_VAR*) NULL); /*lint !e732 !e747*/
321  for (int j = 0; j < i; ++j)
322  {
323  SCIP_VAR* var;
324  (void) SCIPsnprintf(var_name, 255, "E%d_%d", i, j );
325 
326  SCIP_CALL( SCIPcreateVar(scip,
327  &var, // returns new index
328  var_name, // name
329  0.0, // lower bound
330  2.0, // upper bound
331  dist[i][j], // objective
332  SCIP_VARTYPE_INTEGER, // variable type
333  true, // initial
334  false, // forget the rest ...
335  NULL, NULL, NULL, NULL, NULL) ); /*lint !e732 !e747*/
336  SCIP_CALL( SCIPaddVar(scip, var) );
337  arc_var[i][j] = var; /*lint !e732 !e747*/
338  }
339  }
340 
341  /* add arc-routing - tour constraints */
342  char con_name[255];
343  vector< vector<SCIP_CONS*> > arc_con( num_nodes ); /*lint !e732 !e747*/
344  for (int i = 0; i < num_nodes; ++i)
345  {
346  arc_con[i].resize(i, (SCIP_CONS*)NULL); /*lint !e732 !e747*/
347  for (int j = 0; j < i; ++j)
348  {
349  SCIP_CONS* con;
350  (void) SCIPsnprintf(con_name, 255, "A%d_%d", i, j);
351  SCIP_VAR* idx = arc_var[i][j]; /*lint !e732 !e747*/
352  SCIP_Real coeff = -1;
353  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 1, &idx, &coeff,
354  -SCIPinfinity(scip), /* lhs */
355  0.0, /* rhs */
356  true, /* initial */
357  false, /* separate */
358  true, /* enforce */
359  true, /* check */
360  true, /* propagate */
361  false, /* local */
362  true, /* modifiable */
363  false, /* dynamic */
364  false, /* removable */
365  false) ); /* stickingatnode */
366  SCIP_CALL( SCIPaddCons(scip, con) );
367  arc_con[i][j] = con; /*lint !e732 !e747*/
368  }
369  }
370 
371  /* add arc-routing - degree constraints */
372  for (int i = 1; i < num_nodes; ++i)
373  {
374  SCIP_CONS* con;
375  (void) SCIPsnprintf(con_name, 255, "D%d", i);
376  SCIP_CALL( SCIPcreateConsLinear(scip, &con, con_name, 0, 0, 0,
377  2.0, /* lhs */
378  2.0, /* rhs */
379  true, /* initial */
380  false, /* separate */
381  true, /* enforce */
382  true, /* check */
383  true, /* propagate */
384  false, /* local */
385  false, /* modifiable */
386  false, /* dynamic */
387  false, /* removable */
388  false) ); /* stickingatnode */
389  SCIP_CALL( SCIPaddCons(scip, con) );
390  for (int j = 0; j < num_nodes; ++j)
391  {
392  if ( j != i )
393  {
394  SCIP_CALL( SCIPaddCoefLinear(scip, con, i > j ? arc_var[i][j] : arc_var[j][i], 1.0) ); /*lint !e732 !e747*/
395  }
396  }
397  SCIP_CALL( SCIPreleaseCons(scip, &con) );
398  }
399 
400  /* add set packing constraints (Node 0 is the depot) */
401  vector<SCIP_CONS*> part_con(num_nodes, (SCIP_CONS*)NULL); /*lint !e732 !e747*/
402  for (int i = 1; i < num_nodes; ++i)
403  {
404  SCIP_CONS* con = NULL;
405  (void) SCIPsnprintf(con_name, 255, "C%d", i);
406  SCIP_CALL( SCIPcreateConsLinear( scip, &con, con_name, 0, NULL, NULL,
407  1.0, /* lhs */
408  SCIPinfinity(scip), /* rhs */
409  true, /* initial */
410  false, /* separate */
411  true, /* enforce */
412  true, /* check */
413  true, /* propagate */
414  false, /* local */
415  true, /* modifiable */
416  false, /* dynamic */
417  false, /* removable */
418  false /* stickingatnode */ ) );
419  SCIP_CALL( SCIPaddCons(scip, con) );
420  part_con[i] = con; /*lint !e732 !e747*/
421  }
422 
423  /* include VRP pricer */
424  ObjPricerVRP* vrp_pricer_ptr = new ObjPricerVRP(scip, VRP_PRICER_NAME, num_nodes, capacity, demand, dist,
425  arc_var, arc_con, part_con);
426 
427  SCIP_CALL( SCIPincludeObjPricer(scip, vrp_pricer_ptr, true) );
428 
429  /* activate pricer */
430  SCIP_CALL( SCIPactivatePricer(scip, SCIPfindPricer(scip, VRP_PRICER_NAME)) );
431 
432  // SCIP_CALL( SCIPwriteOrigProblem(scip, "vrp_init.lp", "lp", FALSE) );
433 
434 
435  /*************
436  * Solve *
437  *************/
438 
439  SCIP_CALL( SCIPsolve(scip) );
440 
441 
442  /**************
443  * Statistics *
444  *************/
445  SCIP_CALL( SCIPprintStatistics(scip, NULL) );
446 
447  SCIP_CALL( SCIPprintBestSol(scip, NULL, FALSE) );
448 
449 
450 
451  /********************
452  * Deinitialization *
453  ********************/
454 
455  /* release variables */
456  for (int i = 0; i < num_nodes; ++i)
457  {
458  if ( i > 0 )
459  {
460  SCIP_CALL( SCIPreleaseCons(scip, &part_con[i]) );
461  }
462  for (int j = 0; j < i; ++j)
463  {
464  SCIP_CALL( SCIPreleaseVar(scip, &arc_var[i][j]) );
465  SCIP_CALL( SCIPreleaseCons(scip, &arc_con[i][j]) );
466  }
467  }
468 
469 
470  SCIP_CALL( SCIPfree(&scip) );
471 
473 
474  return SCIP_OKAY;
475 }
476 
477 int main(int argc, char** argv)
478 {
479  return execmain(argc, argv) != SCIP_OKAY ? 1 : 0;
480 }
SCIP_RETCODE SCIPcreateConsLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable, SCIP_Bool stickingatnode)
#define BMScheckEmptyMemory()
Definition: memory.h:147
static SCIP_RETCODE execmain(int argc, char **argv)
Definition: main_vrp.cpp:230
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
#define FALSE
Definition: def.h:73
SCIP_PRICER * SCIPfindPricer(SCIP *scip, const char *name)
Definition: scip_pricer.c:302
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_RETCODE SCIPincludeObjPricer(SCIP *scip, scip::ObjPricer *objpricer, SCIP_Bool deleteobject)
Definition: objpricer.cpp:212
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2527
SCIP_VAR ** x
Definition: circlepacking.c:54
SCIP_RETCODE SCIPprintBestSol(SCIP *scip, FILE *file, SCIP_Bool printzeros)
Definition: scip_sol.c:2371
C++ wrapper for default SCIP plugins.
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:283
SCIPInterval sqrt(const SCIPInterval &x)
Definition: pqueue.h:28
VRP pricer plugin.
#define NULL
Definition: lpi_spx1.cpp:155
C++ wrapper classes for SCIP.
void SCIPprintVersion(SCIP *scip, FILE *file)
Definition: scip_general.c:146
#define SCIP_CALL(x)
Definition: def.h:364
Definition: grphload.c:88
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
int main(int argc, char **argv)
Definition: main_vrp.cpp:477
SCIP_RETCODE SCIPcreateVar(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype, SCIP_Bool initial, SCIP_Bool removable, SCIP_DECL_VARDELORIG((*vardelorig)), SCIP_DECL_VARTRANS((*vartrans)), SCIP_DECL_VARDELTRANS((*vardeltrans)), SCIP_DECL_VARCOPY((*varcopy)), SCIP_VARDATA *vardata)
Definition: scip_var.c:105
SCIP_RETCODE SCIPprintStatistics(SCIP *scip, FILE *file)
SCIP_RETCODE SCIPcreateProb(SCIP *scip, const char *name, SCIP_DECL_PROBDELORIG((*probdelorig)), SCIP_DECL_PROBTRANS((*probtrans)), SCIP_DECL_PROBDELTRANS((*probdeltrans)), SCIP_DECL_PROBINITSOL((*probinitsol)), SCIP_DECL_PROBEXITSOL((*probexitsol)), SCIP_DECL_PROBCOPY((*probcopy)), SCIP_PROBDATA *probdata)
Definition: scip_prob.c:107
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1666
static int read_problem(const char *filename, int &num_nodes, int &capacity, vector< int > &demand, vector< vector< int > > &dist)
Definition: main_vrp.cpp:61
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10590
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1252
#define SCIP_Real
Definition: def.h:163
SCIP_VAR ** y
Definition: circlepacking.c:55
SCIP_RETCODE SCIPactivatePricer(SCIP *scip, SCIP_PRICER *pricer)
Definition: scip_pricer.c:375
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2764
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:315
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1110
SCIP_RETCODE SCIPsetIntParam(SCIP *scip, const char *name, int value)
Definition: scip_param.c:497