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