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