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