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