1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 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 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 */
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)41 ObjPricerVRP::ObjPricerVRP(
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. */
~ObjPricerVRP()64 ObjPricerVRP::~ObjPricerVRP()
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 */
SCIP_DECL_PRICERINIT(ObjPricerVRP::scip_init)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 */
pricing(SCIP * scip,bool isfarkas) const97 SCIP_RETCODE ObjPricerVRP::pricing(
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 */
SCIP_DECL_PRICERREDCOST(ObjPricerVRP::scip_redcost)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 */
SCIP_DECL_PRICERFARKAS(ObjPricerVRP::scip_farkas)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 */
add_tour_variable(SCIP * scip,const list<int> & tour) const249 SCIP_RETCODE ObjPricerVRP::add_tour_variable(
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
PQUEUE_KEY__anon143de6db0111::PQUEUE_KEY320 PQUEUE_KEY() : demand(0), length(0.0) {}
321 };
322
operator <(const PQUEUE_KEY & l1,const PQUEUE_KEY & l2)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
NODE_TABLE_DATA__anon143de6db0111::NODE_TABLE_DATA350 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 */
find_shortest_tour(const vector<vector<SCIP_Real>> & length,list<int> & tour) const364 SCIP_Real ObjPricerVRP::find_shortest_tour(
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 }
497