Scippy

SCIP

Solving Constraint Integer Programs

pricer_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-2019 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file pricer_vrp.cpp
17  * @brief VRP pricer plugin
18  * @author Andreas Bley
19  * @author Marc Pfetsch
20  */
21 
22 #include "pricer_vrp.h"
23 #include "pqueue.h"
24 
25 #include <iostream>
26 #include <map>
27 #include <vector>
28 
29 #include "scip/cons_linear.h"
30 
31 using namespace std;
32 using namespace scip;
33 
34 
35 
36 
37 /** Constructs the pricer object with the data needed
38  *
39  * An alternative is to have a problem data class which allows to access the data.
40  */
42  SCIP* scip, /**< SCIP pointer */
43  const char* p_name, /**< name of pricer */
44  const int p_num_nodes, /**< number of nodes */
45  const int p_capacity, /**< vehicle capacity */
46  const vector< int >& p_demand, /**< demand array */
47  const vector< vector<int> >& p_distance, /**< matrix of distances */
48  const vector< vector<SCIP_VAR*> >& p_arc_var, /**< matrix of arc variables */
49  const vector< vector<SCIP_CONS*> >& p_arc_con, /**< matrix of arc constraints */
50  const vector<SCIP_CONS* >& p_part_con /**< array of partitioning constraints */
51  ):
52  ObjPricer(scip, p_name, "Finds tour with negative reduced cost.", 0, TRUE),
53  _num_nodes(p_num_nodes),
54  _capacity(p_capacity),
55  _demand(p_demand),
56  _distance(p_distance),
57  _arc_var(p_arc_var),
58  _arc_con(p_arc_con),
59  _part_con(p_part_con)
60 {}
61 
62 
63 /** Destructs the pricer object. */
65 {}
66 
67 
68 /** initialization method of variable pricer (called after problem was transformed)
69  *
70  * Because SCIP transformes the original problem in preprocessing, we need to get the references to
71  * the variables and constraints in the transformed problem from the references in the original
72  * problem.
73  */
74 SCIP_DECL_PRICERINIT(ObjPricerVRP::scip_init)
75 {
76  for (int i = 0; i < num_nodes(); ++i)
77  {
78  for (int j = 0; j < i; ++j)
79  {
80  SCIP_CALL( SCIPgetTransformedVar(scip, _arc_var[i][j], &_arc_var[i][j]) ); /*lint !e732 !e747*/
81  SCIP_CALL( SCIPgetTransformedCons(scip, _arc_con[i][j], &_arc_con[i][j]) ); /*lint !e732 !e747*/
82  }
83  }
84  for (int i = 1; i < num_nodes(); ++i)
85  {
86  SCIP_CALL( SCIPgetTransformedCons(scip, _part_con[i], &_part_con[i]) ); /*lint !e732 !e747*/
87  }
88 
89  return SCIP_OKAY;
90 } /*lint !e715*/
91 
92 
93 /** perform pricing
94  *
95  * @todo compute shortest length restricted tour w.r.t. duals
96  */
98  SCIP* scip, /**< SCIP data structure */
99  bool isfarkas /**< whether we perform Farkas pricing */
100  ) const
101 {
102  /* allocate array for reduced costs */
103  vector< vector<SCIP_Real> > red_length(num_nodes()); /*lint !e732 !e747*/
104  for (int i = 0; i < num_nodes(); ++i)
105  red_length[i].resize(i, 0.0); /*lint !e732 !e747*/
106 
107  /* compute reduced-cost arc lengths store only lower triangualar matrix, i.e., red_length[i][j] only for i > j */
108  if ( isfarkas )
109  {
110  for (int i = 0; i < num_nodes(); ++i)
111  {
112  assert( i == 0 || part_con(i) != 0 );
113  for (int j = 0; j < i; ++j)
114  {
115  SCIP_Real r = 0.0;
116  assert( arc_con(i,j) != 0 );
117 
118  r -= SCIPgetDualfarkasLinear(scip, arc_con(i,j));
119  if ( j != 0 )
120  r -= 0.5 * SCIPgetDualfarkasLinear(scip, part_con(j));
121  if ( i != 0 )
122  r -= 0.5 * SCIPgetDualfarkasLinear(scip, part_con(i));
123  red_length[i][j] = r; /*lint !e732 !e747*/
124  }
125  }
126  }
127  else
128  {
129  for (int i = 0; i < num_nodes(); ++i)
130  {
131  assert( i == 0 || part_con(i) != 0 );
132  for (int j = 0; j < i; ++j)
133  {
134  SCIP_Real r = 0.0;
135  assert( arc_con(i,j) != 0 );
136 
137  r -= SCIPgetDualsolLinear(scip, arc_con(i,j));
138  if ( j != 0 )
139  r -= 0.5 * SCIPgetDualsolLinear(scip, part_con(j));
140  if ( i != 0 )
141  r -= 0.5 * SCIPgetDualsolLinear(scip, part_con(i));
142  red_length[i][j] = r; /*lint !e732 !e747*/
143  }
144  }
145  }
146 
147 #ifdef SCIP_OUTPUT
148  if ( isfarkas )
149  {
150  SCIPinfoMessage(scip, NULL, "dual ray solution:\n");
151  for (int i = 0; i < num_nodes(); ++i)
152  {
153  for (int j = 0; j < i; ++j)
154  SCIPinfoMessage(scip, NULL, "arc_%d_%d: %g\n", i, j, SCIPgetDualfarkasLinear(scip, arc_con(i,j)));
155  }
156 
157  for (int i = 1; i < num_nodes(); ++i)
158  SCIPinfoMessage(scip, NULL, "part_%d: %g\n", i, SCIPgetDualfarkasLinear(scip, part_con(i)));
159 
160  for (int i = 0; i < num_nodes(); ++i)
161  {
162  for (int j = 0; j < i; ++j)
163  SCIPinfoMessage(scip, NULL, "length_%d_%d: %g\n", i, j, red_length[i][j]);
164  }
165  }
166  else
167  {
168  SCIPinfoMessage(scip, NULL, "dual solution:\n");
169  for (int i = 0; i < num_nodes(); ++i)
170  {
171  for (int j = 0; j < i; ++j)
172  SCIPinfoMessage(scip, NULL, "arc_%d_%d: %g\n", i, j, SCIPgetDualsolLinear(scip, arc_con(i,j)));
173  }
174 
175  for (int i = 1; i < num_nodes(); ++i)
176  SCIPinfoMessage(scip, NULL, "part_%d: %g\n", i, SCIPgetDualsolLinear(scip, part_con(i)));
177 
178  for (int i = 0; i < num_nodes(); ++i)
179  {
180  for (int j = 0; j < i; ++j)
181  SCIPinfoMessage(scip, NULL, "length_%d_%d: %g\n", i, j, red_length[i][j]);
182  }
183  }
184 #endif
185 
186  /* compute shortest length restricted tour w.r.t. reduced-cost arc length */
187  list<int> tour;
188  SCIP_Real reduced_cost = find_shortest_tour(red_length, tour);
189 
190  /* add tour variable */
191  if ( SCIPisNegative(scip, reduced_cost) )
192  {
193  return add_tour_variable(scip, tour);
194  }
195 
196 #ifdef SCIP_OUTPUT
197  SCIP_CALL( SCIPwriteTransProblem(scip, "vrp.lp", "lp", FALSE) );
198 #endif
199 
200  return SCIP_OKAY;
201 }
202 
203 
204 
205 /** Pricing of additional variables if LP is feasible.
206  *
207  * - get the values of the dual variables you need
208  * - construct the reduced-cost arc lengths from these values
209  * - find the shortest admissible tour with respect to these lengths
210  * - if this tour has negative reduced cost, add it to the LP
211  *
212  * possible return values for *result:
213  * - SCIP_SUCCESS : at least one improving variable was found, or it is ensured that no such variable exists
214  * - SCIP_DIDNOTRUN : the pricing process was aborted by the pricer, there is no guarantee that the current LP solution is optimal
215  */
216 SCIP_DECL_PRICERREDCOST(ObjPricerVRP::scip_redcost)
217 {
218  SCIPdebugMsg(scip, "call scip_redcost ...\n");
219 
220  /* set result pointer, see above */
221  *result = SCIP_SUCCESS;
222 
223  /* call pricing routine */
224  SCIP_CALL( pricing(scip, false) );
225 
226  return SCIP_OKAY;
227 } /*lint !e715*/
228 
229 
230 /** Pricing of additional variables if LP is infeasible.
231  *
232  * - get the values of the dual Farks multipliers you need
233  * - construct the reduced-cost arc lengths from these values
234  * - find the shortest admissible tour with respect to these lengths
235  * - if this tour has negative reduced cost, add it to the LP
236  */
237 SCIP_DECL_PRICERFARKAS(ObjPricerVRP::scip_farkas)
238 {
239  SCIPdebugMsg(scip, "call scip_farkas ...\n");
240 
241  /* call pricing routine */
242  SCIP_CALL( pricing(scip, true) );
243 
244  return SCIP_OKAY;
245 } /*lint !e715*/
246 
247 
248 /** add tour variable to problem */
250  SCIP* scip, /**< SCIP data structure */
251  const list<int>& tour /**< list of nodes in tour */
252  ) const
253 {
254  /* create meaningful variable name */
255  char tmp_name[255];
256  char var_name[255];
257  (void) SCIPsnprintf(var_name, 255, "T");
258  for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it) /*lint !e1702*/
259  {
260  strncpy(tmp_name, var_name, 255);
261  (void) SCIPsnprintf(var_name, 255, "%s_%d", tmp_name, *it);
262  }
263  SCIPdebugMsg(scip, "new variable <%s>\n", var_name);
264 
265  /* create the new variable: Use upper bound of infinity such that we do not have to care about
266  * the reduced costs of the variable in the pricing. The upper bound of 1 is implicitly satisfied
267  * due to the set partitioning constraints.
268  */
269  SCIP_VAR* var;
270  SCIP_CALL( SCIPcreateVar(scip, &var, var_name,
271  0.0, // lower bound
272  SCIPinfinity(scip), // upper bound
273  0.0, // objective
274  SCIP_VARTYPE_CONTINUOUS, // variable type
275  false, false, NULL, NULL, NULL, NULL, NULL) );
276 
277  /* add new variable to the list of variables to price into LP (score: leave 1 here) */
278  SCIP_CALL( SCIPaddPricedVar(scip, var, 1.0) );
279 
280  /* add coefficient into the set partition constraints */
281  for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it) /*lint !e1702*/
282  {
283  assert( 0 <= *it && *it < num_nodes() );
284  SCIP_CALL( SCIPaddCoefLinear(scip, part_con(*it), var, 1.0) );
285  }
286 
287  /* add coefficient into arc routing constraints */
288  int last = 0;
289  for (list<int>::const_iterator it = tour.begin(); it != tour.end(); ++it) /*lint !e1702*/
290  {
291  assert( 0 <= *it && *it < num_nodes() );
292  SCIP_CALL( SCIPaddCoefLinear(scip, arc_con(last, *it), var, 1.0) );
293  last = *it;
294  }
295  SCIP_CALL( SCIPaddCoefLinear(scip, arc_con(last, 0), var, 1.0 ) );
296 
297  /* cleanup */
298  SCIP_CALL( SCIPreleaseVar(scip, &var) );
299 
300  return SCIP_OKAY;
301 }
302 
303 
304 /** Computes a shortest admissible tour with respect to the given lengths. The function must return
305  * the computed tour via the parameter tour and the length (w.r.t. given lengths) of this tour as
306  * return parameter. The returned tour must be the ordered list of customer nodes contained in the
307  * tour (i.e., 2-5-7 for the tour 0-2-5-7-0).
308  */
309 namespace
310 {
311 
312 /* types needed for prioity queue -------------------- */
313 static const SCIP_Real eps = 1e-9;
314 
315 struct PQUEUE_KEY
316 {
317  int demand;
318  SCIP_Real length;
319 
320  PQUEUE_KEY() : demand(0), length(0.0) {}
321 };
322 
323 bool operator< (const PQUEUE_KEY& l1, const PQUEUE_KEY& l2)
324 {
325  if ( l1.demand < l2.demand )
326  return true;
327  if ( l1.demand > l2.demand )
328  return false;
329  if ( l1.length < l2.length-eps )
330  return true;
331  /* not needed, since we return false anyway:
332  if ( l1.length > l2.length+eps )
333  return false;
334  */
335  return false;
336 }
337 
338 typedef int PQUEUE_DATA; // node
339 typedef pqueue<PQUEUE_KEY,PQUEUE_DATA> PQUEUE;
340 typedef PQUEUE::pqueue_item PQUEUE_ITEM;
341 
342 
343 /* types needed for dyn. programming table */
344 struct NODE_TABLE_DATA
345 {
346  SCIP_Real length;
347  int predecessor;
348  PQUEUE::pqueue_item queue_item;
349 
350  NODE_TABLE_DATA( ) : length(0.0), predecessor(-1), queue_item( NULL ) {}
351 };
352 
353 typedef int NODE_TABLE_KEY; // demand
354 typedef std::map< NODE_TABLE_KEY, NODE_TABLE_DATA > NODE_TABLE;
355 }
356 
357 
358 /** return negative reduced cost tour (uses restricted shortest path dynamic programming algorithm)
359  *
360  * The algorithm uses the priority queue implementation in pqueue.h. SCIP's implementation of
361  * priority queues cannot be used, since it currently does not support removal of elements that are
362  * not at the top.
363  */
365  const vector< vector<SCIP_Real> >& length, /**< matrix of lengths */
366  list<int>& tour /**< list of nodes in tour */
367  ) const
368 {
369  tour.clear();
370 
371  SCIPdebugMessage("Enter RSP - capacity: %d\n", capacity());
372 
373  /* begin algorithm */
374  PQUEUE PQ;
375  vector< NODE_TABLE > table(num_nodes()); /*lint !e732 !e747*/
376 
377  /* insert root node (start at node 0) */
378  PQUEUE_KEY queue_key;
379  PQUEUE_DATA queue_data = 0;
380  PQUEUE_ITEM queue_item = PQ.insert(queue_key, queue_data);
381 
382  NODE_TABLE_KEY table_key = 0;
383  NODE_TABLE_DATA table_entry;
384 
385  /* run Dijkstra-like updates */
386  while ( ! PQ.empty() )
387  {
388  /* get front queue entry */
389  queue_item = PQ.top();
390  queue_key = PQ.get_key (queue_item);
391  queue_data = PQ.get_data(queue_item);
392  PQ.pop();
393 
394  /* get corresponding node and node-table key */
395  const int curr_node = queue_data;
396  const SCIP_Real curr_length = queue_key.length;
397  const int curr_demand = queue_key.demand;
398 
399  /* stop as soon as some negative length tour was found */
400  if ( curr_node == 0 && curr_length < -eps )
401  break;
402 
403  /* stop as soon don't create multi-tours */
404  if ( curr_node == 0 && curr_demand != 0 )
405  continue;
406 
407  /* update all active neighbors */
408  for (int next_node = 0; next_node < num_nodes(); ++next_node)
409  {
410  if ( next_node == curr_node )
411  continue;
412  if ( have_edge( next_node, curr_node ) == false )
413  continue;
414 
415  const int next_demand = curr_demand + demand(next_node);
416 
417  if ( next_demand > capacity() )
418  continue;
419 
420  const SCIP_Real next_length = curr_length + ( curr_node > next_node ? /*lint !e732 !e747*/
421  length[curr_node][next_node] : /*lint !e732 !e747*/
422  length[next_node][curr_node] ); /*lint !e732 !e747*/
423 
424  NODE_TABLE& next_table = table[next_node]; /*lint !e732 !e747*/
425 
426  /* check if new table entry would be dominated */
427  bool skip = false;
428  list<NODE_TABLE::iterator> dominated;
429 
430  for (NODE_TABLE::iterator it = next_table.begin(); it != next_table.end() && ! skip; ++it) /*lint !e1702*/
431  {
432  if ( next_demand >= it->first && next_length >= it->second.length - eps )
433  skip = true;
434 
435  if ( next_demand <= it->first && next_length <= it->second.length + eps )
436  dominated.push_front( it );
437  }
438  if ( skip )
439  continue;
440 
441  /* remove dominated table and queue entries */
442  for (list<NODE_TABLE::iterator>::iterator it = dominated.begin(); it != dominated.end(); ++it) /*lint !e1702*/
443  {
444  PQ.remove( (*it)->second.queue_item );
445  next_table.erase( *it );
446  }
447 
448  /* insert new table and queue entry */
449  queue_key.demand = next_demand;
450  queue_key.length = next_length;
451  queue_data = next_node;
452 
453  queue_item = PQ.insert(queue_key, queue_data);
454 
455  table_key = next_demand;
456  table_entry.length = next_length;
457  table_entry.predecessor = curr_node;
458  table_entry.queue_item = queue_item;
459 
460  next_table[table_key] = table_entry;
461 
462 #ifdef SCIP_OUTPUT
463  printf("new entry node = %d demand = %d length = %g pref = %d\n", next_node, next_demand, next_length, curr_node);
464 #endif
465  }
466  }
467 
468  SCIPdebugMessage("Done RSP DP.\n");
469 
470  table_entry.predecessor = -1;
471  table_entry.length = 0;
472  int curr_node = 0;
473 
474  /* find most negative tour */
475  for (NODE_TABLE::iterator it = table[0].begin(); it != table[0].end(); ++it) /*lint !e1702 !e732 !e747*/
476  {
477  if ( it->second.length < table_entry.length )
478  {
479  table_key = it->first;
480  table_entry = it->second;
481  }
482  }
483  SCIP_Real tour_length = table_entry.length;
484 
485  while ( table_entry.predecessor > 0 )
486  {
487  table_key -= demand(curr_node);
488  curr_node = table_entry.predecessor;
489  tour.push_front(curr_node);
490  table_entry = table[curr_node][table_key]; /*lint !e732 !e747*/
491  }
492 
493  SCIPdebugMessage("Leave RSP tour length = %g\n", tour_length);
494 
495  return tour_length;
496 }
#define NULL
Definition: def.h:246
C++ wrapper for variable pricer.
Definition: objpricer.h:42
int num_nodes() const
Definition: pricer_vrp.h:87
SCIP_RETCODE SCIPgetTransformedVar(SCIP *scip, SCIP_VAR *var, SCIP_VAR **transvar)
Definition: scip_var.c:1442
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1251
virtual SCIP_DECL_PRICERREDCOST(scip_redcost)
#define FALSE
Definition: def.h:72
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10253
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
#define TRUE
Definition: def.h:71
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
#define SCIPdebugMessage
Definition: pub_message.h:77
bool have_edge(const int i, const int j) const
Definition: pricer_vrp.h:142
#define SCIPdebugMsg
Definition: scip_message.h:88
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:279
real eps
SCIP_CONS * arc_con(const int i, const int j) const
Definition: pricer_vrp.h:125
int capacity() const
Definition: pricer_vrp.h:93
double find_shortest_tour(const vector< vector< double > > &length, list< int > &tour) const
Definition: pricer_vrp.cpp:364
SCIP_Real SCIPgetDualsolLinear(SCIP *scip, SCIP_CONS *cons)
virtual SCIP_DECL_PRICERINIT(scip_init)
Definition: pqueue.h:28
VRP pricer plugin.
virtual SCIP_DECL_PRICERFARKAS(scip_farkas)
#define SCIP_CALL(x)
Definition: def.h:358
virtual ~ObjPricerVRP()
Definition: pricer_vrp.cpp:64
SCIP_RETCODE SCIPgetTransformedCons(SCIP *scip, SCIP_CONS *cons, SCIP_CONS **transcons)
Definition: scip_cons.c:1688
SCIP_Real SCIPgetDualfarkasLinear(SCIP *scip, SCIP_CONS *cons)
SCIP_RETCODE SCIPwriteTransProblem(SCIP *scip, const char *filename, const char *extension, SCIP_Bool genericnames)
Definition: scip_prob.c:699
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:104
Constraint handler for linear constraints in their most general form, .
class for priority queues
SCIP_RETCODE pricing(SCIP *scip, bool isfarkas) const
Definition: pricer_vrp.cpp:97
SCIP_Real * r
Definition: circlepacking.c:50
SCIP_RETCODE SCIPaddPricedVar(SCIP *scip, SCIP_VAR *var, SCIP_Real score)
Definition: scip_prob.c:1789
ObjPricerVRP(SCIP *scip, const char *p_name, const int p_num_nodes, const int p_capacity, const vector< int > &p_demand, const vector< vector< int > > &p_distance, const vector< vector< SCIP_VAR *> > &p_arc_var, const vector< vector< SCIP_CONS *> > &p_arc_con, const vector< SCIP_CONS * > &p_part_con)
Definition: pricer_vrp.cpp:41
#define SCIP_Real
Definition: def.h:157
SCIP_RETCODE add_tour_variable(SCIP *scip, const list< int > &tour) const
Definition: pricer_vrp.cpp:249
int demand(const int i) const
Definition: pricer_vrp.h:99
SCIP_CONS * part_con(const int i) const
Definition: pricer_vrp.h:134