VRPH
1.0
|
00001 00002 // // 00003 // This file is part of the VRPH software package for // 00004 // generating solutions to vehicle routing problems. // 00005 // VRPH was developed by Chris Groer (cgroer@gmail.com). // 00006 // // 00007 // (c) Copyright 2010 Chris Groer. // 00008 // All Rights Reserved. VRPH is licensed under the // 00009 // Common Public License. See LICENSE file for details. // 00010 // // 00012 00013 #include "VRPH.h" 00014 #define VRPH_MAX_CYCLES 500 00015 00016 void VRPH_version() 00017 { 00018 printf("--------------------------------------------\n"); 00019 printf("VRPH, version 1.0.1\nCopyright 2010 Chris Groer\nDistributed under Common Public License 1.0\n"); 00020 printf("--------------------------------------------\n\n"); 00021 } 00022 00023 00024 VRP::VRP(int n) 00025 { 00029 00030 int i,j; 00031 00032 num_nodes=n; 00033 num_original_nodes=n; 00034 total_demand=0; 00035 num_days=0; 00036 00037 next_array = new int[n+2]; 00038 pred_array = new int[n+2]; 00039 route_num = new int[n+2]; 00040 route = new VRPRoute[n+2]; 00041 routed = new bool[n+2]; 00042 best_sol_buff = new int[n+2]; 00043 current_sol_buff = new int[n+2]; 00044 search_space = new int[n+2]; 00045 nodes = new VRPNode[n+2]; 00046 // Add an extra spot for the VRPH_DEPOT and for the dummy node 00047 00048 symmetric=true; 00049 // Set to false only when we encounter FULL_MATRIX file 00050 00051 forbid_tiny_moves=true; 00052 // Default is to allow these moves 00053 00054 d=NULL; 00055 // The distance matrix is allocated when the problem is loaded 00056 fixed=new bool*[n+2]; 00057 fixed[0]=new bool[(n+2)*(n+2)]; 00058 for(i=1;i<n+2;i++) 00059 fixed[i]=fixed[i-1]+(n+2); 00060 for(i=0;i<n+2;i++) 00061 { 00062 routed[i]=false; 00063 for(j=0;j<n+2;j++) 00064 fixed[i][j]=false; 00065 } 00066 00067 // Set these to default values--they may change once 00068 // we read the file. 00069 00070 min_vehicles=-1; 00071 has_service_times=false; 00072 max_route_length=VRP_INFINITY; 00073 orig_max_route_length=VRP_INFINITY; 00074 total_route_length=0.0; 00075 best_known=VRP_INFINITY; 00076 depot_normalized=false; 00077 // Will be set to true if we shift nodes so VRPH_DEPOT is at origin 00078 00079 // These are for record-to-record travel 00080 record = 0.0; 00081 deviation = VRPH_DEFAULT_DEVIATION; 00082 00083 // For keeping track of the statistics 00084 00085 for(i=0;i<NUM_HEURISTICS;i++) 00086 { 00087 num_evaluations[i]=0; 00088 num_moves[i]=0; 00089 00090 } 00091 00092 total_service_time = 0.0; 00093 00094 // Create the solution warehouse 00095 this->solution_wh=new VRPSolutionWarehouse(NUM_ELITE_SOLUTIONS,n); 00096 00097 this->tabu_list=new VRPTabuList(MAX_VRPH_TABU_LIST_SIZE); 00098 this->route_wh=NULL; 00099 00100 // Set this to true only if we have valid coordinates 00101 // This is valid only for plotting the solution 00102 can_display=false; 00103 00104 00105 } 00106 00107 VRP::VRP(int n, int ndays) 00108 { 00112 00113 int i,j; 00114 00115 num_nodes=n; 00116 num_original_nodes=n; 00117 total_demand=0; 00118 num_days=ndays; 00119 00120 next_array = new int[n+2]; 00121 pred_array = new int[n+2]; 00122 route_num = new int[n+2]; 00123 route = new VRPRoute[n+2]; 00124 routed = new bool[n+2]; 00125 best_sol_buff = new int[n+2]; 00126 current_sol_buff = new int[n+2]; 00127 search_space = new int[n+2]; 00128 nodes = new VRPNode[n+2]; 00129 // Add an extra spot for the VRPH_DEPOT and for the dummy node 00130 00131 forbid_tiny_moves=true; 00132 // Default is to forbid these moves 00133 00134 d=NULL; 00135 // The distance matrix is allocated when the problem is loaded 00136 fixed=new bool*[n+2]; 00137 fixed[0]=new bool[(n+2)*(n+2)]; 00138 for(i=1;i<n+2;i++) 00139 fixed[i]=fixed[i-1]+(n+2); 00140 for(i=0;i<n+2;i++) 00141 { 00142 routed[i]=false; 00143 for(j=0;j<n+2;j++) 00144 fixed[i][j]=false; 00145 } 00146 00147 // Set these to default values--they may change once 00148 // we read the file. 00149 00150 min_vehicles=-1; 00151 has_service_times=false; 00152 max_route_length=VRP_INFINITY; 00153 orig_max_route_length=VRP_INFINITY; 00154 total_route_length=0.0; 00155 best_known=VRP_INFINITY; 00156 00157 depot_normalized=false; 00158 // Will be set to true if we shift nodes so VRPH_DEPOT is at origin 00159 00160 // These are for record-to-record travel 00161 record = 0.0; 00162 deviation = VRPH_DEFAULT_DEVIATION; 00163 00164 // For keeping track of the statistics 00165 00166 for(i=0;i<NUM_HEURISTICS;i++) 00167 { 00168 num_evaluations[i]=0; 00169 num_moves[i]=0; 00170 00171 } 00172 00173 total_service_time = 0.0; 00174 00175 // Create the solution warehouse 00176 this->solution_wh=new VRPSolutionWarehouse(NUM_ELITE_SOLUTIONS,n); 00177 00178 this->tabu_list=new VRPTabuList(MAX_VRPH_TABU_LIST_SIZE); 00179 this->route_wh=NULL; 00180 00181 // Now allocate d days worth of storage at each of the nodes 00182 for(i=0;i<=n+1;i++) 00183 { 00184 this->nodes[i].daily_demands=new int[ndays+1]; 00185 this->nodes[i].daily_service_times=new double[ndays+1]; 00186 } 00187 00188 // Set this to true only if we have valid coordinates 00189 // This is valid only for plotting the solution 00190 can_display=false; 00191 00192 } 00193 00194 VRP::~VRP() 00195 { 00199 00200 delete [] this->best_sol_buff; 00201 delete [] this->current_sol_buff; 00202 delete [] this->d[0]; 00203 delete [] this->d; 00204 delete [] this->fixed[0]; 00205 delete [] this->fixed; 00206 delete [] this->next_array; 00207 delete [] this->search_space; 00208 delete [] this->nodes; 00209 delete [] this->pred_array; 00210 delete [] this->route; 00211 delete [] this->route_num; 00212 delete [] this->routed; 00213 delete this->solution_wh; 00214 delete this->tabu_list; 00215 00216 } 00217 00218 00219 // Accessor functions for private data 00220 int VRP::get_num_nodes() 00221 { 00225 00226 return this->num_nodes; 00227 } 00228 00229 double VRP::get_total_route_length() 00230 { 00234 00235 return this->total_route_length; 00236 } 00237 00238 double VRP::get_total_service_time() 00239 { 00243 00244 return this->total_service_time; 00245 } 00246 00247 double VRP::get_best_sol_buff(int *sol_buff) 00248 { 00254 00255 memcpy(sol_buff,this->best_sol_buff,(this->num_nodes+1)*sizeof(int)); 00256 return this->best_total_route_length; 00257 } 00258 00259 double VRP::get_best_total_route_length() 00260 { 00264 00265 return this->best_total_route_length; 00266 } 00267 00268 int VRP::get_total_number_of_routes() 00269 { 00273 00274 return this->total_number_of_routes; 00275 } 00276 00277 int VRP::get_num_original_nodes() 00278 { 00282 00283 return this->num_original_nodes; 00284 } 00285 00286 int VRP::get_num_days() 00287 { 00291 00292 return this->num_days; 00293 } 00294 00295 double VRP::get_best_known() 00296 { 00297 return this->best_known; 00298 00299 } 00300 00301 int VRP::get_max_veh_capacity() 00302 { 00303 return this->max_veh_capacity; 00304 } 00305 00306 double VRP::get_max_route_length() 00307 { 00308 return this->max_route_length; 00309 } 00310 00311 void VRP::set_best_total_route_length(double val) 00312 { 00313 this->best_total_route_length=val; 00314 return; 00315 00316 } 00317 bool VRP::clone(VRP *W) 00318 { 00322 00323 this->balance_parameter=W->balance_parameter; 00324 this->best_known=W->best_known; 00325 this->best_total_route_length=W->best_total_route_length; 00326 // Copy the best solution buffer 00327 memcpy(this->best_sol_buff, W->best_sol_buff, (sizeof(int))*(W->num_nodes+2)); 00328 00329 // Copy the currnet solution buffer 00330 memcpy(this->current_sol_buff, W->current_sol_buff, (sizeof(int))*(W->num_nodes+2)); 00331 00332 this->coord_type=W->coord_type; 00333 this->d=W->d; // OK to just copy the pointers here 00334 this->deviation=W->deviation; 00335 this->display_type=W->display_type; 00336 this->edge_weight_format=W->edge_weight_format; 00337 this->dummy_index=W->dummy_index; 00338 this->has_service_times=W->has_service_times; 00339 this->matrix_size=W->matrix_size; 00340 this->max_route_length=W->max_route_length; 00341 this->neighbor_list_size=W->neighbor_list_size; 00342 this->nodes = W->nodes; // OK to just copy the pointers here? 00343 00344 this->num_nodes=W->num_nodes; 00345 this->total_route_length=W->total_route_length; 00346 this->orig_max_route_length=W->orig_max_route_length; 00347 this->orig_max_veh_capacity=W->orig_max_veh_capacity; 00348 00349 this->problem_type=W->problem_type; 00350 this->record=W->record; 00351 00352 00353 this->search_size=W->search_size; 00354 00355 this->temperature=W->temperature; 00356 this->total_number_of_routes=W->total_number_of_routes; 00357 this->total_service_time=W->total_service_time; 00358 this->max_veh_capacity=W->max_veh_capacity; 00359 this->violation=W->violation; 00360 this->edge_weight_type=W->edge_weight_type; 00361 00362 // Assume that the current_solution has been sent to W-> 00363 W->export_solution_buff(W->current_sol_buff); 00364 this->import_solution_buff(W->current_sol_buff); 00365 return true; 00366 00367 } 00368 00369 00370 void VRP::create_distance_matrix(int type) 00371 { 00378 00379 00380 int i,j,k,n; 00381 n=this->num_nodes; 00382 00383 if(type==VRPH_EXPLICIT) 00384 { 00385 // We have presumably already loaded in the distance matrix 00386 // But still add in the service times!! 00387 if(this->has_service_times) 00388 { 00389 for(i=0;i<=n+1;i++) 00390 { 00391 for(j=0;j<=n+1;j++) 00392 { 00393 //printf("Adding service time of %f\n",this->nodes[j].service_time); 00394 this->d[i][j]+= .5*this->nodes[j].service_time; 00395 this->d[i][j]+= .5*this->nodes[i].service_time; 00396 } 00397 } 00398 } 00399 // We have presumably already loaded in the distance matrix! 00400 return; 00401 } 00402 00403 n= this->num_original_nodes; 00404 k=0; 00405 00406 // Otherwise construct the matrix - we store the whole thing even though 00407 // it is symmetric as we found that this was quite a bit faster... 00408 00409 for(i=0;i<=n+1;i++) 00410 { 00411 for(j=0;j<=n+1;j++) 00412 // CSG: Changed Nov. 20 - add in .5 service time from i and j 00413 this->d[i][j]=VRPDistance(type, this->nodes[i].x,this->nodes[i].y,this->nodes[j].x, 00414 this->nodes[j].y) + 00415 0.5*this->nodes[i].service_time + 00416 0.5*this->nodes[j].service_time; 00417 } 00418 00419 00420 return; 00421 00422 } 00423 00424 void VRP::create_neighbor_lists(int nsize) 00425 { 00430 00431 if(nsize>num_nodes ) 00432 { 00433 fprintf(stderr,"Requested neighbor list size is greater than num_nodes!\n%d>%d\n", 00434 nsize,num_nodes); 00435 report_error("%s: Neighbor list error!!\n",__FUNCTION__); 00436 } 00437 00438 if(nsize>MAX_NEIGHBORLIST_SIZE ) 00439 { 00440 fprintf(stderr,"Requested neighbor list size is greater than MAX_NEIGHBORLIST_SIZE!\n%d>%d\n", 00441 nsize,MAX_NEIGHBORLIST_SIZE); 00442 report_error("%s: Neighbor list error!!\n",__FUNCTION__); 00443 } 00444 00445 00446 int i,j,n,b; 00447 int k; 00448 double dd; 00449 double max; 00450 int maxpos; 00451 VRPNeighborElement *NList; 00452 00453 n= num_nodes; 00454 00455 // Set the neighbor_list_size value 00456 neighbor_list_size=nsize; 00457 00458 NList=new VRPNeighborElement[nsize]; 00459 00460 // First, do the neighbor_list for the VRPH_DEPOT 00461 max=0.0; 00462 // NEW 00463 maxpos=0; 00464 00465 for(i=1;i<=nsize;i++) 00466 { 00467 NList[i-1].val=d[VRPH_DEPOT][i]; 00468 00469 NList[i-1].position = i; 00470 00471 // keep track of the largest value 00472 if(NList[i-1].val > max) 00473 { 00474 max=NList[i-1].val; 00475 maxpos=i-1; 00476 } 00477 00478 } 00479 00480 // Now go through the rest of the nodes. 00481 for(i=nsize+1;i<=n;i++) 00482 { 00483 if(d[VRPH_DEPOT][i]<max) 00484 { 00485 // Replace element in position maxpos 00486 NList[maxpos].val=d[VRPH_DEPOT][i]; 00487 //nodes[i].depot_distance; 00488 NList[maxpos].position=i; 00489 } 00490 // Now find the new max 00491 max=0.0; 00492 for(b=0;b<nsize;b++) 00493 { 00494 if(NList[b].val>max) 00495 { 00496 maxpos=b; 00497 max=NList[b].val; 00498 } 00499 } 00500 } 00501 00502 qsort(NList,nsize,sizeof(VRPNeighborElement), VRPNeighborCompare); 00503 i=0; 00504 for(b=0;b<nsize;b++) 00505 { 00506 nodes[VRPH_DEPOT].neighbor_list[b].val=NList[b].val; 00507 nodes[VRPH_DEPOT].neighbor_list[b].position=NList[b].position; 00508 00509 # if NEIGHBOR_DEBUG 00510 printf("(%d,%d,%f) \n",VRPH_DEPOT,nodes[VRPH_DEPOT].neighbor_list[b].position, 00511 nodes[VRPH_DEPOT].neighbor_list[b].val); 00512 #endif 00513 00514 00515 } 00516 00517 # if NEIGHBOR_DEBUG 00518 printf("\n\n"); 00519 #endif 00520 00521 00522 k=0; 00523 // Done with NeighborList for VRPH_DEPOT. Now do it for the rest of the nodes 00524 for(i=1;i<=n;i++) 00525 { 00526 // Loop through all the nodes and create their Neighbor Lists 00527 // First initialize the NList 00528 for(j=0;j<nsize;j++) 00529 { 00530 NList[j].position=VRP_INFINITY; 00531 NList[j].val=VRP_INFINITY; 00532 } 00533 00534 00535 // VRPH_DEPOT can be in NList, so set NList[0] to depot_distance 00536 NList[0].position=VRPH_DEPOT; 00537 NList[0].val=d[i][VRPH_DEPOT]; 00538 00539 00540 for(j=1;j<i;j++) 00541 { 00542 #if NEIGHBOR_DEBUG 00543 printf("Loop 1; i=%d; in j loop(%d)\n",i,j); 00544 #endif 00545 00546 dd=d[i][j]; 00547 00548 if(j<nsize) 00549 { 00550 // Construct the original nsize long neighbor list 00551 NList[j].val=dd; 00552 NList[j].position=j; 00553 // Update max size 00554 if(NList[j].val>max) 00555 { 00556 max=NList[j].val; 00557 maxpos=j; 00558 00559 } 00560 } 00561 else 00562 { 00563 // j >= nsize 00564 if(dd<max) 00565 { 00566 NList[maxpos].val=dd; 00567 NList[maxpos].position=j; 00568 // Now find the new max element; 00569 max=0.0; 00570 for(b=0;b<nsize;b++) 00571 { 00572 if(NList[b].val>max) 00573 { 00574 max=NList[b].val; 00575 maxpos=b; 00576 } 00577 } 00578 } 00579 } 00580 00581 } 00582 00583 00584 for(j=i+1;j<=n;j++) 00585 { 00586 #if NEIGHBOR_DEBUG 00587 printf("Loop 2; i=%d; in j loop(%d)\n",i,j); 00588 #endif 00589 00590 dd=d[i][j]; 00591 00592 if(j<=nsize)// was < 00593 { 00594 // Construct the original nsize long neighbor list 00595 NList[j-1].val=dd; 00596 NList[j-1].position=j; 00597 // Update max size 00598 if(NList[j-1].val>max) 00599 { 00600 max=NList[j-1].val; 00601 maxpos=j-1; 00602 00603 } 00604 } 00605 else 00606 { 00607 if(dd<max) 00608 { 00609 // Replace maxpos with this one 00610 NList[maxpos].val=dd; 00611 NList[maxpos].position=j; 00612 // Now find the new max element; 00613 max=0.0; 00614 for(b=0;b<nsize;b++) 00615 { 00616 if(NList[b].val>max) 00617 { 00618 max=NList[b].val; 00619 maxpos=b; 00620 } 00621 } 00622 } 00623 } 00624 00625 } 00626 00627 00628 #if NEIGHBOR_DEBUG 00629 printf("Node %d neighbor list before sorting\n",i); 00630 for(b=0;b<nsize;b++) 00631 { 00632 printf("NList[%d]=%d,%f[%f]\n",b,NList[b].position,NList[b].val,d[i][NList[b].position]); 00633 } 00634 #endif 00635 // Now sort the NList and stick the resulting list into 00636 // the neighbor list for node i 00637 qsort (NList, nsize, sizeof(VRPNeighborElement), VRPNeighborCompare); 00638 00639 for(b=0;b<nsize;b++) 00640 { 00641 nodes[i].neighbor_list[b].position=NList[b].position; 00642 nodes[i].neighbor_list[b].val=NList[b].val; 00643 if(i== nodes[i].neighbor_list[b].position) 00644 { 00645 fprintf(stderr,"ERROR:: Node %d is in it's own neighbor list!!\n", i); 00646 fprintf(stderr,"i=%d; %d[%d],%f[%f])\n",i,nodes[i].neighbor_list[b].position, 00647 NList[b].position,nodes[i].neighbor_list[b].val,NList[b].val); 00648 00649 report_error("%s: Error creating neighbor lists\n",__FUNCTION__); 00650 } 00651 #if NEIGHBOR_DEBUG 00652 00653 printf("NList[%d]=%d,%f[%f]",b,nodes[i].neighbor_list[b].position,nodes[i].neighbor_list[b].val,d[i][nodes[i].neighbor_list[b].position]); 00654 00655 #endif 00656 } 00657 00658 // Hackk to make the VRPH_DEPOT in everyone's list 00659 //nodes[i].neighbor_list[nsize-1].position=VRPH_DEPOT; 00660 //nodes[i].neighbor_list[nsize-1].val=d[i][VRPH_DEPOT]; 00661 00662 } 00663 00664 delete [] NList; 00665 00666 return; 00667 00668 } 00669 00670 00671 bool VRP::check_feasibility(VRPViolation *VV) 00672 { 00678 00679 int i; 00680 bool is_feasible=true; 00681 // We will set this to false if we find a violation 00682 00683 VV->capacity_violation=-VRP_INFINITY; 00684 VV->length_violation=-VRP_INFINITY; 00685 00686 this->normalize_route_numbers(); 00687 00688 for(i=1;i<=this->total_number_of_routes;i++) 00689 { 00690 if(this->route[i].length>this->orig_max_route_length) 00691 { 00692 // Length violation 00693 if(this->route[i].length - this->orig_max_route_length > VV->length_violation) 00694 VV->length_violation = this->route[i].length - this->orig_max_route_length ; 00695 00696 is_feasible=false; 00697 } 00698 00699 if(this->route[i].load>this->orig_max_veh_capacity) 00700 { 00701 // Load violation 00702 if(this->route[i].load - this->orig_max_veh_capacity > VV->capacity_violation) 00703 VV->capacity_violation = this->route[i].load - this->orig_max_veh_capacity ; 00704 00705 is_feasible=false; 00706 00707 } 00708 00709 } 00710 00711 return is_feasible; 00712 00713 } 00714 00715 00716 00717 00718 void VRP::refresh_routes() 00719 { 00724 00725 00726 int i, n, cnt; 00727 double len=0; 00728 double rlen=0; 00729 double tot_len=0; 00730 int next_node; 00731 int current_node, current_route, route_start, current_start, current_end; 00732 int total_load = 0; 00733 int current_load = 0; 00734 00735 00736 n=num_nodes; 00737 00738 i = 1; 00739 cnt = 0; 00740 route_start = -next_array[VRPH_DEPOT]; 00741 00742 current_node = route_start; 00743 current_route = route_num[current_node]; 00744 00745 current_start = route[current_route].start; 00746 current_end = route[current_route].end; 00747 00748 total_load+=nodes[current_node].demand; 00749 current_load+=nodes[current_node].demand; 00750 len+=d[VRPH_DEPOT][current_node]; 00751 rlen+=d[VRPH_DEPOT][current_node]; 00752 00753 while(route_start != 0 && i<num_nodes+1) 00754 { 00755 // When we finally get a route beginning at 0, this is the last route 00756 // and there is no next route, so break out 00757 if(next_array[current_node]==current_node) 00758 { 00759 // We've entered a loop 00760 fprintf(stderr,"(2)Self loop found in next array(%d)\n",current_node); 00761 report_error("%s: Error in refresh_routes()\n",__FUNCTION__); 00762 } 00763 00764 if(next_array[current_node]==0) 00765 { 00766 // We're at the end of the Solution! 00767 00768 len+=d[current_node][VRPH_DEPOT]; 00769 rlen+=d[current_node][VRPH_DEPOT]; 00770 current_route=route_num[current_node]; 00771 route[current_route].length = rlen; 00772 route[current_route].load = current_load; 00773 total_route_length=len; 00774 00775 // We're done now 00776 return; 00777 } 00778 00779 if(next_array[current_node]>0) 00780 { 00781 // Next node is somewhere in the middle of a route 00782 00783 next_node = next_array[current_node]; 00784 len+=d[current_node][next_node]; 00785 rlen+=d[current_node][next_node]; 00786 current_node=next_node; 00787 total_load+=nodes[current_node].demand; 00788 current_load+=nodes[current_node].demand; 00789 cnt++; 00790 } 00791 else 00792 { 00793 // We must have a non-positive "next" node indicating the beginning of a new route 00794 //len+=nodes[current_node].depot_distance; 00795 len+=d[current_node][VRPH_DEPOT]; 00796 //rlen+=nodes[current_node].depot_distance; 00797 rlen+=d[current_node][VRPH_DEPOT]; 00798 tot_len+=rlen; 00799 current_route=route_num[current_node]; 00800 current_end = route[current_route].end; 00801 00802 route[current_route].length = rlen; 00803 route[current_route].load = current_load; 00804 00805 i++; 00806 route_start = - (next_array[current_node]); 00807 current_route = route_num[route_start]; 00808 current_start = route[current_route].start; 00809 current_end = route[current_route].end; 00810 if(route_start != current_start) 00811 { 00812 fprintf(stderr,"Route %d: %d != %d\n",current_route, route_start, current_start); 00813 report_error("%s: Error in refresh_routes()\n",__FUNCTION__); 00814 } 00815 00816 current_node = route_start; 00817 total_load+=nodes[current_node].demand; 00818 // reset current_load to 0 00819 current_load=nodes[current_node].demand; 00820 len+=d[VRPH_DEPOT][current_node]; 00821 rlen=d[VRPH_DEPOT][current_node]; 00822 cnt++; 00823 } 00824 } 00825 00826 00827 return; 00828 00829 00830 } 00831 void VRP::create_pred_array() 00832 { 00836 int i,j; 00837 00838 i=VRPH_DEPOT; 00839 j=next_array[i]; 00840 while(j!=VRPH_DEPOT) 00841 { 00842 if(j>0) 00843 pred_array[j]=i; 00844 else 00845 pred_array[-j]=-i; 00846 00847 i=VRPH_ABS(j); 00848 j=next_array[i]; 00849 00850 } 00851 // The VRPH_DEPOT 00852 pred_array[j]=-i; 00853 return; 00854 00855 00856 00857 00858 00859 /*for(i=0;i<=num_nodes;i++) 00860 { 00861 j=next_array[i]; 00862 if(j>0) 00863 pred_array[j] = i; 00864 else 00865 pred_array[-j] = -i; 00866 } 00867 return;*/ 00868 } 00869 00870 00871 00872 00873 00874 bool VRP::get_segment_info(int a, int b, struct VRPSegment *S) 00875 { 00885 00886 00887 if(a==b) 00888 { 00889 // Degenerate case... 00890 S->segment_start=b; 00891 S->segment_end=a; 00892 S->len=0; //nodes[a].service_time??; 00893 S->load=nodes[a].demand; 00894 S->num_custs=1; 00895 return true; 00896 00897 } 00898 00899 S->len=0; 00900 S->segment_start=a; 00901 S->segment_end=b; 00902 S->num_custs=0; 00903 00904 00905 // Check special cases 00906 if(a==VRPH_DEPOT) 00907 { 00908 // Safe to assume we want the segment 0-i-j-...-k-b-... 00909 00910 S->segment_start=route[route_num[b]].start; 00911 S->segment_end=b; 00912 S->len+=d[VRPH_DEPOT][S->segment_start]; 00913 00914 } 00915 00916 if(b==VRPH_DEPOT) 00917 { 00918 // Safe to assume we want the segment ...a-i-j-...-k-0 00919 00920 S->segment_start=a; 00921 S->segment_end=route[route_num[a]].end; 00922 00923 } 00924 00925 // Now calculate the length, load, and # customer on this segment 00926 00927 int current_node = S->segment_start; 00928 int next_node; 00929 00930 S->load=nodes[current_node].demand; 00931 if(current_node!=dummy_index) 00932 S->num_custs++;; 00933 00934 while(current_node != S->segment_end) 00935 { 00936 next_node = VRPH_MAX(next_array[current_node],0); 00937 S->len+=d[current_node][next_node]; 00938 current_node=next_node; 00939 S->load+=nodes[next_node].demand; 00940 if(current_node!= dummy_index) 00941 S->num_custs += 1; 00942 } 00943 00944 if(b==VRPH_DEPOT) 00945 S->len+=d[S->segment_end][VRPH_DEPOT]; 00946 00947 00948 00949 return true; 00950 00951 } 00952 00953 int VRP::get_string_end(int a, int len) 00954 { 00960 00961 int ctr=1; 00962 int current_node; 00963 00964 current_node=a; 00965 while(ctr<len) 00966 { 00967 current_node= next_array[current_node]; 00968 if(current_node<0) 00969 { 00970 // The string of length len beginning at a is not contained in a single route 00971 return -1; 00972 } 00973 ctr++; 00974 00975 } 00976 00977 return current_node; 00978 } 00979 00980 void VRP::reverse_route(int i) 00981 { 00987 00988 00989 int current_node, start_node, last_node, prev_route, next_route, temp; 00990 int orig_start, orig_end ; 00991 00992 if(i<=0) 00993 { 00994 // Illegal! 00995 fprintf(stderr,"Reversing route of negative index?? i=%d\n",i); 00996 report_error("%s: Can't reverse route!\n",__FUNCTION__); 00997 } 00998 00999 orig_end = route[i].end; 01000 orig_start = route[i].start; 01001 01002 start_node=orig_start; 01003 current_node=start_node; 01004 01005 // We begin with the first node in the route - we will change its next_array[] 01006 // value at the very end 01007 01008 temp = next_array[current_node]; 01009 prev_route = -pred_array[current_node]; 01010 // So next[prev_route] should equal start_node 01011 // We will change this so that next[prev_route]=-last_node 01012 01013 pred_array[current_node] = temp; 01014 current_node = temp; 01015 while( (temp = next_array[current_node]) >0) 01016 { 01017 next_array[current_node] = pred_array[current_node]; 01018 pred_array[current_node] = temp; 01019 current_node = temp; 01020 } 01021 01022 // current_node is now the end of the original route, so we can finish the reversal 01023 01024 temp= pred_array[current_node]; 01025 last_node = current_node; 01026 next_route = -next_array[last_node]; 01027 01028 route[i].end=orig_start; 01029 route[i].start=orig_end; 01030 01031 #if REVERSE_DEBUG 01032 printf("start_node: %d; last_node: %d; prev_route: %d; next_route: %d\n",start_node,last_node, prev_route,next_route); 01033 #endif 01034 01035 // Final set of updates 01036 next_array[prev_route]=-last_node; 01037 pred_array[next_route]=-start_node; 01038 next_array[start_node]=-next_route; 01039 pred_array[last_node]= -prev_route; 01040 next_array[last_node]=temp; 01041 next_array[prev_route]=-last_node; 01042 01043 // Need to update length if asymmetric 01044 if(!this->symmetric) 01045 { 01046 // Update length as it possibly changed! 01047 this->refresh_routes(); 01048 } 01049 01050 return; 01051 } 01052 01053 01054 01055 bool VRP::postsert_dummy(int i) 01056 { 01057 // This function inserts dummy node after i 01058 01059 int pre_i, post_i, dummy; 01060 int start, end, start_i, end_i; 01061 int n, i_route; 01062 //double tu, uv, tv, iu, ju, ij, u_loss, i_gain, i_change, i_length, u_length; 01063 01064 // Special cases first 01065 01066 if(i<=VRPH_DEPOT || i>matrix_size) 01067 report_error("%s: input doesn't make sense\n",__FUNCTION__); 01068 01069 n= num_nodes; 01070 dummy= dummy_index; 01071 01072 i_route= route_num[i]; 01073 01074 start_i= route[i_route].start; 01075 end_i= route[i_route].end; 01076 01077 start=start_i; 01078 if(end_i==i) 01079 // i is last in the route 01080 end=dummy; 01081 else 01082 end=end_i; 01083 01084 // pre_i is what used to be before i 01085 pre_i= pred_array[i]; 01086 // post_i is what used to be after i 01087 post_i= next_array[i]; 01088 01089 next_array[i] = dummy; 01090 next_array[dummy] = post_i; 01091 01092 // Now need to update pred_array as well! 01093 01094 // u's predecessor is now i 01095 pred_array[dummy]=i; 01096 01097 // The element who used to be after i is now preceded by u 01098 if(post_i>=0) 01099 pred_array[post_i]=dummy; 01100 else 01101 // post_i 01102 pred_array[-post_i]=-dummy; 01103 01104 //start_array[dummy]=start; 01105 //end_array[dummy]=end; 01106 01107 // Update i_route information 01108 route_num[dummy]=i_route; 01109 route[i_route].end=end; 01110 route[i_route].start=start; 01111 01112 return true; 01113 } 01114 01115 bool VRP::presert_dummy(int i) 01116 { 01117 // This function inserts a dummy node BEFORE node i in 01118 // whatever route node i is currently in 01119 01120 int pre_i, post_i; 01121 int start, end, start_i, end_i, dummy; 01122 int n, i_route; 01123 01124 01125 n= num_nodes ; 01126 dummy = dummy_index; // dummy node has index (n+1) 01127 01128 if(i<=VRPH_DEPOT) 01129 report_error("%s: bad index\n",__FUNCTION__); 01130 01131 i_route = route_num[i]; 01132 01133 start_i= route[i_route].start; 01134 end_i= route[i_route].end; 01135 01136 01137 // Get the start and end of i's route 01138 start=start_i; 01139 end=end_i; 01140 01141 if(start==i) 01142 // i is first in the route, so dummy is now the first node 01143 start=dummy; 01144 01145 // pre_i is what used to be before i 01146 pre_i= pred_array[i]; 01147 // post_i is what used to be after i 01148 post_i= next_array[i]; 01149 01150 01151 // dummy is now followed by i 01152 next_array[dummy]=i; 01153 pred_array[i]=dummy; 01154 pred_array[dummy]=pre_i; 01155 01156 // The element who used to be after i is now preceded by dummy 01157 if(pre_i>0)// was >=!! 01158 next_array[pre_i]=dummy; 01159 else 01160 // post_i 01161 next_array[VRPH_ABS(pre_i)]=-dummy; 01162 01163 // Update i_route information 01164 route_num[dummy]=i_route; 01165 route[i_route].end=end; 01166 route[i_route].start=start; 01167 01168 // Add in the relevant Data fields for the dummy node!! 01169 01170 return true; 01171 01172 } 01173 bool VRP::remove_dummy() 01174 { 01175 int pre_d, post_d, d_route, dummy, d_start, d_end, n; 01176 01177 n= num_nodes; 01178 dummy= dummy_index; 01179 01180 post_d= next_array[dummy]; 01181 pre_d= pred_array[dummy]; 01182 01183 if(post_d > dummy || post_d < -dummy || pre_d > dummy || pre_d < -dummy) 01184 { 01185 fprintf(stderr,"post_d= %d; pre_d=%d\n",post_d,pre_d); 01186 report_error("%s: invalid indices\n",__FUNCTION__); 01187 } 01188 01189 d_route= route_num[dummy]; 01190 d_start= route[d_route].start; 01191 d_end= route[d_route].end; 01192 01193 if(d_start==dummy) 01194 { 01195 if(post_d<0) 01196 { 01197 report_error("%s: post_d error in removal\n",__FUNCTION__); 01198 } 01199 01200 route[d_route].start=post_d; 01201 } 01202 01203 if(d_end==dummy) 01204 { 01205 if(pre_d<0) 01206 { 01207 report_error("%s: pre_d error in removal\n",__FUNCTION__); 01208 } 01209 route[d_route].end=pre_d; 01210 } 01211 01212 next_array[VRPH_ABS(pre_d)]=post_d; 01213 01214 if(d_start==dummy) 01215 next_array[VRPH_ABS(pre_d)]=-post_d; 01216 01217 01218 pred_array[VRPH_ABS(post_d)]=pre_d; 01219 01220 if(d_end==dummy) 01221 pred_array[VRPH_ABS(post_d)]=-pre_d; 01222 01223 01224 return true; 01225 } 01226 bool VRP::create_default_routes() 01227 { 01234 01235 01236 int i,n; 01237 bool is_feasible=true; 01238 // No violations yet... 01239 violation.capacity_violation = 0; 01240 violation.length_violation = 0; 01241 01242 01243 total_route_length=0; 01244 01245 n = this->num_original_nodes; 01246 // First, create the initial routes: VRPH_DEPOT - node - VRPH_DEPOT for all nodes 01247 01248 routed[VRPH_DEPOT]=true; 01249 01250 next_array[VRPH_DEPOT] = -1; 01251 for(i=1;i<=n;i++) 01252 { 01253 next_array[i] = -(i+1); 01254 total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]); 01255 01256 route_num[i]=i; 01257 route[i].start=i; 01258 route[i].end=i; 01259 route[i].load= nodes[i].demand; 01260 route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]; 01261 01262 01263 // Check capacities 01264 if(route[i].load > max_veh_capacity) 01265 is_feasible=false; 01266 if(route[i].length > max_route_length) 01267 is_feasible=false; 01268 01269 route[i].num_customers=1; 01270 01271 routed[i]=true; 01272 01273 } 01274 01275 01276 next_array[n]=VRPH_DEPOT; 01277 01278 route_num[VRPH_DEPOT]=0; 01279 01280 // Now create the associated pred_array implied by the newly created next_array 01281 create_pred_array(); 01282 01283 total_number_of_routes=n; 01284 01285 01286 01287 01288 01289 if(is_feasible == false) 01290 { 01291 // We had infeasible routes - record the violation in the Solution and 01292 // return false; 01293 01294 01295 for(i=1;i<=n;i++) 01296 { 01297 routed[i]=false; 01298 // Set routed status to false since default routes are infeasible 01299 01300 // Check capacities 01301 if(route[i].load > max_veh_capacity) 01302 { 01303 // Update the worst violations 01304 01305 printf("Default routes load violation: %d > %d\n",route[i].load, 01306 max_veh_capacity); 01307 01308 if( (route[i].load - max_veh_capacity) > 01309 violation.capacity_violation) 01310 { 01311 violation.capacity_violation = 01312 (route[i].load - max_veh_capacity); 01313 } 01314 01315 } 01316 if(route[i].length > max_route_length) 01317 { 01318 // Update the worst violations 01319 01320 printf("Default routes length violation: %f > %f\n",route[i].length, 01321 max_route_length); 01322 01323 if( (route[i].length -max_route_length ) > 01324 violation.length_violation) 01325 { 01326 violation.length_violation = 01327 (route[i].length -max_route_length ); 01328 } 01329 01330 01331 } 01332 01333 01334 } 01335 01336 return false; 01337 } 01338 else 01339 { 01340 // All routes were feasible 01341 for(i=1;i<=n;i++) 01342 { 01343 routed[i]=true; 01344 } 01345 return true; 01346 } 01347 01348 01349 } 01350 01351 01352 bool VRP::create_default_routes(int day) 01353 { 01360 01361 int i,n; 01362 bool is_feasible=true; 01363 // No violations yet... 01364 violation.capacity_violation = 0; 01365 violation.length_violation = 0; 01366 total_route_length=0; 01367 01368 n = this->num_original_nodes; 01369 this->num_nodes=n; 01370 01371 // First, create the initial routes: VRPH_DEPOT - node - VRPH_DEPOT for all nodes 01372 01373 routed[VRPH_DEPOT]=true; 01374 next_array[VRPH_DEPOT] = -1; 01375 for(i=1;i<=n;i++) 01376 { 01377 01378 next_array[i] = -(i+1); 01379 total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]); 01380 01381 route_num[i]=i; 01382 route[i].start=i; 01383 route[i].end=i; 01384 route[i].load= nodes[i].demand; 01385 route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]; 01386 route[i].num_customers=1; 01387 routed[i]=true; 01388 01389 } 01390 01391 01392 next_array[n]=VRPH_DEPOT; 01393 01394 route_num[VRPH_DEPOT]=0; 01395 01396 // Now create the associated pred_array implied by the newly created next_array 01397 create_pred_array(); 01398 01399 total_number_of_routes=n; 01400 01401 // Now eject the nodes that don't require service on this day 01402 for(i=1;i<=n;i++) 01403 { 01404 if(this->nodes[i].daily_demands[day]==-1)// indicates no service required 01405 this->eject_node(i); 01406 } 01407 01408 // Now check for feasibility; 01409 this->normalize_route_numbers(); 01410 01411 // Check capacities 01412 for(i=1;i<=this->total_number_of_routes;i++) 01413 { 01414 if(route[i].load > max_veh_capacity) 01415 is_feasible=false; 01416 if(route[i].length > max_route_length) 01417 is_feasible=false; 01418 } 01419 01420 if(is_feasible == false) 01421 { 01422 // We had infeasible routes - record the violation in the Solution and 01423 // return false; 01424 01425 01426 for(i=1;i<=n;i++) 01427 { 01428 routed[i]=false; 01429 // Set routed status to false since default routes are infeasible 01430 01431 // Check capacities 01432 if(route[i].load > max_veh_capacity) 01433 { 01434 // Update the worst violations 01435 01436 printf("Default routes load violation: %d > %d\n",route[i].load, 01437 max_veh_capacity); 01438 01439 if( (route[i].load - max_veh_capacity) > 01440 violation.capacity_violation) 01441 { 01442 violation.capacity_violation = 01443 (route[i].load - max_veh_capacity); 01444 } 01445 01446 } 01447 if(route[i].length > max_route_length) 01448 { 01449 // Update the worst violations 01450 01451 printf("Default routes length violation: %f > %f\n",route[i].length, 01452 max_route_length); 01453 01454 if( (route[i].length -max_route_length ) > 01455 violation.length_violation) 01456 { 01457 violation.length_violation = 01458 (route[i].length -max_route_length ); 01459 } 01460 } 01461 } 01462 return false; 01463 } 01464 else 01465 return true; 01466 01467 } 01468 01469 01470 int VRP::count_num_routes() 01471 { 01475 01476 int current,next; 01477 int num=0; 01478 01479 current=VRPH_DEPOT; 01480 next=-1; 01481 01482 while(next!=VRPH_DEPOT) 01483 { 01484 next= next_array[current]; 01485 if(next<0) 01486 { 01487 // We're in a new route 01488 num++; 01489 current=-next; 01490 } 01491 else 01492 current=next; 01493 } 01494 01495 return num; 01496 01497 01498 } 01499 01500 bool VRP::perturb() 01501 { 01505 01506 int i, n, pre, post; 01507 int current,next; 01508 n= num_nodes; 01509 VRPNeighborElement *v; 01510 VRPMove M; 01511 01512 Postsert Postsert; 01513 Presert Presert; 01514 01515 v=new VRPNeighborElement[n]; 01516 01517 current=VRPH_ABS(next_array[VRPH_DEPOT]); 01518 i=0; 01519 while(current!=VRPH_DEPOT) 01520 { 01521 pre= VRPH_MAX(VRPH_DEPOT,pred_array[current]); 01522 post= VRPH_MAX(VRPH_DEPOT, next_array[current]); 01523 01524 // Define s[i]= d(pre(i),i) + d(i,next(i))-d(pre(i),next(i)) 01525 v[i].val=((double) nodes[current].demand) / 01526 (VRPH_EPSILON+d[pre][current]+d[current][post]-d[pre][post]); 01527 v[i].position=current; 01528 i++; 01529 next=VRPH_ABS(next_array[current]); 01530 current=next; 01531 } 01532 01533 01534 // Now sort the v[i]'s in ascending order of the val's 01535 01536 qsort(v,n,sizeof(VRPNeighborElement), VRPNeighborCompare); 01537 01538 int m=(int) (VRPH_MIN(30,n/5)); 01539 01540 01541 // Now insert the first m nodes from the sorted list into new locations. 01542 01543 int a,b,c,j,k,node1=0, node2=0,b_route,b_load,k_demand, jj, ll; 01544 double best_savings,b_len,jk,kl,jl; 01545 01546 01547 for(j=0;j<m;j++) 01548 { 01549 best_savings=VRP_INFINITY; 01550 01551 k=v[j].position; 01552 k_demand= nodes[k].demand; 01553 01554 // Node k will be moved 01555 01556 // Former neighboring nodes of k 01557 jj=VRPH_MAX(pred_array[k],0); 01558 ll=VRPH_MAX(next_array[k],0); 01559 01560 jk= d[jj][k]; 01561 kl= d[k][ll]; 01562 jl= d[jj][ll]; 01563 01564 01565 // Node k will be moved - look for a new location 01566 // OLD SITUATION: jj-k-ll 01567 b=VRPH_ABS(next_array[VRPH_DEPOT]); 01568 while(b!=VRPH_DEPOT) 01569 { 01570 b_route= route_num[b]; 01571 b_load= route[b].load; 01572 b_len= route[b].length; 01573 01574 a=VRPH_MAX(pred_array[b],0); 01575 c=VRPH_MAX(next_array[b],0); 01576 01577 if(a!=k && b!=k && c!=k)// This guarantees a new location 01578 { 01579 01580 if(a!=VRPH_DEPOT) 01581 { 01582 if(Postsert.evaluate(this,k,a,&M)==true && M.savings<best_savings) 01583 { 01584 best_savings=M.savings; 01585 node1=a;node2=b; 01586 } 01587 } 01588 else 01589 { 01590 01591 if(Presert.evaluate(this,k,b,&M)==true && M.savings<best_savings) 01592 { 01593 best_savings=M.savings; 01594 node1=a;node2=b; 01595 } 01596 01597 01598 01599 } 01600 01601 // Now try the edge b-c 01602 if(b!=VRPH_DEPOT) 01603 { 01604 if(Postsert.evaluate(this,k,b,&M)==true && M.savings<best_savings) 01605 { 01606 best_savings=M.savings; 01607 node1=b;node2=c; 01608 } 01609 } 01610 else 01611 { 01612 // b is the VRPH_DEPOT 01613 if(Presert.evaluate(this,k,c,&M)==true && M.savings<best_savings) 01614 { 01615 best_savings=M.savings; 01616 node1=b;node2=c; 01617 } 01618 01619 01620 } 01621 01622 } 01623 // Advance b 01624 b=VRPH_ABS(next_array[b]); 01625 01626 } 01627 if(best_savings!=VRP_INFINITY) 01628 { 01629 // Now insert k in between node1 and node2 01630 if(node1!=VRPH_DEPOT) 01631 Postsert.move(this,k,node1); 01632 else 01633 Presert.move(this,k,node2); 01634 } 01635 01636 } 01637 01638 01639 01640 delete [] v; 01641 return true; 01642 } 01643 01644 01645 bool VRP::eject_node(int j) 01646 { 01651 01652 if(j<=VRPH_DEPOT || routed[j]==false ) 01653 { 01654 fprintf(stderr,"Tried to eject index %d\n",j); 01655 report_error("%s: VRPH_DEPOT or unrouted index used in eject node\n",__FUNCTION__); 01656 } 01657 01658 01659 int k=j; 01660 01661 int c, e, k_route, c_route, e_route, k_start, k_end,flag; 01662 double change,ck,ke,ce; 01663 01664 flag=0; 01665 01666 // k is no longer routed 01667 routed[k]=false; 01668 01669 // Get k's current location - between c and e 01670 01671 c= pred_array[k]; 01672 e= next_array[k]; 01673 01674 k_route= route_num[k]; 01675 c_route= route_num[VRPH_ABS(c)]; 01676 e_route= route_num[VRPH_ABS(e)]; 01677 01678 k_start= route[k_route].start; 01679 k_end= route[k_route].end; 01680 01681 if(k_start==k && k_end!=k) 01682 { 01683 // c is in a different route 01684 route[k_route].start=VRPH_ABS(e); 01685 next_array[VRPH_ABS(c)]=-VRPH_ABS(e); 01686 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c); 01687 01688 } 01689 01690 if(k_end==k && k_start!=k) 01691 { 01692 // e is in a different route 01693 route[k_route].end=c; 01694 next_array[VRPH_ABS(c)]=-VRPH_ABS(e); 01695 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c); 01696 } 01697 01698 if(k_end==k && k_start==k) 01699 { 01700 // k is its own route 01701 next_array[VRPH_ABS(c)]=-VRPH_ABS(e); 01702 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c); 01703 flag=1; 01704 } 01705 01706 if(k_start!=k && k_end!=k) 01707 { 01708 // k was interior--easy case 01709 next_array[c]=e; 01710 pred_array[e]=c; 01711 01712 } 01713 01714 01715 01716 // Now adjust the objective function value 01717 // and the route length and load, etc. 01718 01719 // old route c-k-e 01720 // new route c-e 01721 01722 if(e<0) 01723 e=VRPH_DEPOT; 01724 if(c<0) 01725 c=VRPH_DEPOT; 01726 01727 ce= this->d[c][e]; 01728 ck= this->d[c][k]; 01729 ke= this->d[k][e]; 01730 01731 change=ce-(ck+ke); 01732 total_route_length+=change; 01733 route[k_route].length+=change; 01734 route[k_route].load-= nodes[k].demand; 01735 01736 // Removing a node from route k_route 01737 route[k_route].num_customers--; 01738 num_nodes--; 01739 01740 if(flag) 01741 // We removed a singleton route 01742 this->total_number_of_routes--; 01743 01744 route_num[k]=-1; 01745 normalize_route_numbers(); 01746 01747 return true; 01748 } 01749 01750 bool VRP::eject_route(int r, int *route_buff) 01751 { 01756 01757 int current,cnt; 01758 01759 // First copy the route into the buffer 01760 current=this->route[r].start; 01761 cnt=0; 01762 while(current>0) 01763 { 01764 route_buff[cnt]=current; 01765 cnt++; 01766 current=this->next_array[current]; 01767 } 01768 route_buff[cnt]=-1; 01769 // Terminate with a -1 01770 01771 01772 01773 // Now eject the nodes 01774 for(int i=0;i<cnt;i++) 01775 this->eject_node(route_buff[i]); 01776 01777 01778 return true; 01779 01780 } 01781 bool VRP::check_move(VRPMove *M, int rules) 01782 { 01783 01788 01789 01790 double savings; 01791 01792 savings=M->savings; 01793 01794 01795 if(this->forbid_tiny_moves) 01796 { 01797 // See if it is a "meaningless" move 01798 if(savings>-VRPH_EPSILON && savings < VRPH_EPSILON) 01799 return false; 01800 } 01801 01802 01803 01804 if( (rules & VRPH_FREE) == VRPH_FREE ) 01805 { 01806 // use with care! 01807 return true; 01808 } 01809 01810 01811 01812 01813 if( (rules & VRPH_DOWNHILL) == VRPH_DOWNHILL) 01814 { 01815 if(savings<-VRPH_EPSILON ) 01816 return true; 01817 else 01818 return false; 01819 01820 } 01821 01822 01823 01824 01825 if( rules & VRPH_RECORD_TO_RECORD ) 01826 { 01827 if(savings<=-VRPH_EPSILON) 01828 return true; 01829 01830 // o/w savings is positive but must be less than deviation*record 01831 01832 01833 01834 if(!has_service_times) 01835 { 01836 if(total_route_length+savings<= (1+deviation)*record) 01837 return true; 01838 else 01839 return false; 01840 01841 } 01842 else 01843 { 01844 // We have service times so remove the service time from the 01845 // deviation calculation 01846 01847 if( (total_route_length - total_service_time) + savings <= 01848 ((1+deviation)*(record-total_service_time))) 01849 return true; 01850 else 01851 return false; 01852 } 01853 } 01854 01855 if( rules & VRPH_SIMULATED_ANNEALING ) 01856 { 01857 if(M->evaluated_savings==true) 01858 return true; // We already checked the random acceptance 01859 01860 // Otherwise, determine whether to accept or not. 01861 if( exp(- (M->savings / this->temperature)) > lcgrand(10) ) 01862 return true; 01863 else 01864 return false; 01865 01866 } 01867 01868 01869 01870 report_error("%s: didn't return yet!\n",__FUNCTION__); 01871 01872 return false; 01873 01874 01875 } 01876 01877 01878 bool VRP::is_feasible(VRPMove *M, int rules) 01879 { 01885 01886 for(int i=0;i<M->num_affected_routes;i++) 01887 { 01888 if( (M->route_lens[i]>this->max_route_length) || (M->route_loads[i]>this->max_veh_capacity) ) 01889 return false; 01890 } 01891 01892 return true; 01893 01894 } 01895 bool VRP::inject_node(int j) 01896 { 01901 01902 if(j==VRPH_DEPOT) 01903 report_error("%s: Can't inject VRPH_DEPOT!!\n",__FUNCTION__); 01904 01905 int edge[4]; 01906 double costs[4]; 01907 this->find_cheapest_insertion(j,edge,costs,VRPH_USE_NEIGHBOR_LIST); 01908 01909 this->insert_node(j,edge[0],edge[1]); 01910 01911 return true; 01912 01913 } 01914 01915 bool VRP::insert_node(int j, int i, int k) 01916 { 01921 01922 if(routed[j] || !routed[i] || !routed[k]) 01923 { 01924 fprintf(stderr,"insert_node(%d,%d,%d)\n",j,i,k); 01925 report_error("%s: Improper nodes used in insert_node\n",__FUNCTION__); 01926 } 01927 01928 double increase; 01929 int r; 01930 routed[j]=true; 01931 01932 if(i==k && k==VRPH_DEPOT) 01933 { 01934 01935 int last_node=VRPH_ABS(pred_array[VRPH_DEPOT]); 01936 this->next_array[last_node]=-j; 01937 this->pred_array[j]=-last_node; 01938 this->next_array[j]=VRPH_DEPOT; 01939 this->pred_array[VRPH_DEPOT]=-j; 01940 01941 increase=this->d[0][j]+this->d[j][0]; 01942 this->total_number_of_routes++; 01943 this->route_num[j]=this->total_number_of_routes; 01944 this->route[total_number_of_routes].length=increase; 01945 this->route[total_number_of_routes].load=nodes[j].demand; 01946 this->route[total_number_of_routes].num_customers=1; 01947 this->route[total_number_of_routes].start=j; 01948 this->route[total_number_of_routes].end=j; 01949 this->num_nodes++; 01950 this->total_route_length+=increase; 01951 01952 return true; 01953 01954 01955 01956 } 01957 01958 if(i!=VRPH_DEPOT && k!=VRPH_DEPOT) 01959 { 01960 01961 // i-k is interior edge 01962 if(VRPH_MAX(VRPH_DEPOT,next_array[i]) != k || VRPH_MAX(VRPH_DEPOT,pred_array[k])!=i) 01963 { 01964 fprintf(stderr,"edge doesn't exist: next(i(%d)) != k(%d)\n",i,k); 01965 report_error("%s\n",__FUNCTION__); 01966 } 01967 01968 increase=d[i][j]+d[j][k]-d[i][k]; 01969 r=route_num[i]; 01970 01971 01972 num_nodes++; 01973 next_array[i]=j; 01974 pred_array[j]=i; 01975 next_array[j]=k; 01976 pred_array[k]=j; 01977 route_num[j]=r; 01978 route[r].length+=increase; 01979 route[r].load+= nodes[j].demand; 01980 route[r].num_customers++; 01981 total_route_length+=increase; 01982 01983 return true; 01984 } 01985 01986 if(i==VRPH_DEPOT) 01987 { 01988 01989 increase=d[i][j]+d[j][k]-d[i][k]; 01990 01991 r=route_num[k]; 01992 01993 01994 int pre=VRPH_ABS(pred_array[k]); 01995 01996 // pre-VRPH_DEPOT-j-k situation... 01997 num_nodes++; 01998 01999 next_array[pre]=-j; 02000 pred_array[j]=-pre; 02001 next_array[j]=k; 02002 pred_array[k]=j; 02003 route_num[j]=r; 02004 route[r].start=j; 02005 route[r].length+=increase; 02006 route[r].load+= nodes[j].demand; 02007 route[r].num_customers++; 02008 02009 total_route_length+=increase; 02010 return true; 02011 02012 } 02013 02014 if(k==VRPH_DEPOT) 02015 { 02016 02017 increase=d[i][j]+d[j][k]-d[i][k]; 02018 r=route_num[i]; 02019 02020 int post=VRPH_ABS(next_array[i]); 02021 // i-j-VRPH_DEPOT-post situation... 02022 02023 num_nodes++; 02024 02025 next_array[i]=j; 02026 pred_array[j]=i; 02027 next_array[j]=-post; 02028 pred_array[post]=-j; 02029 route_num[j]=r; 02030 route[r].length+=increase; 02031 route[r].load+= nodes[j].demand; 02032 route[r].end=j; 02033 route[r].num_customers++; 02034 total_route_length+=increase; 02035 return true; 02036 02037 } 02038 02039 return false; 02040 02041 02042 } 02043 02044 void VRP::find_cheapest_insertion(int j, int *edge, double *costs, int rules) 02045 { 02059 02060 int h,i,k,m, best_route, new_route, next_node; 02061 double ij,ik,jk,min_feasible_increase, increase, min_increase; 02062 02063 this->normalize_route_numbers(); 02064 best_route = -1; 02065 k=-1; 02066 02067 // Set the increase to be a singleton route since this must be feasible 02068 02069 min_feasible_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT]; 02070 min_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT]; 02071 edge[0]=edge[1]=edge[2]=edge[3]=VRPH_DEPOT; 02072 02073 02074 if(!(rules & VRPH_USE_NEIGHBOR_LIST)) 02075 { 02076 i=VRPH_DEPOT; 02077 next_node=VRPH_ABS(next_array[i]); 02078 int cnt=0; 02079 while(next_node != VRPH_DEPOT) 02080 { 02081 cnt++; 02082 if(next_node > 0) 02083 { 02084 k=next_node; 02085 02086 // Try to put j between i and k 02087 02088 ij= d[i][j]; 02089 ik= d[i][k]; 02090 jk= d[j][k]; 02091 increase=ij+jk-ik; 02092 02093 if(increase<min_increase) 02094 { 02095 min_increase=increase; 02096 edge[2]=i; 02097 edge[3]=k; 02098 } 02099 if(increase<min_feasible_increase) 02100 { 02101 // Check feasibility 02102 new_route= route_num[k]; 02103 02104 if( (route[new_route].length+increase <= max_route_length) && 02105 (route[new_route].load + nodes[j].demand <= max_veh_capacity) ) 02106 { 02107 edge[0]=i; 02108 edge[1]=k; 02109 best_route=new_route; 02110 min_feasible_increase=increase; 02111 02112 } 02113 02114 } 02115 // Move on to the next node. 02116 i=k; 02117 02118 } 02119 else 02120 { 02121 // next_node <= 0, implying the edges i-VRPH_DEPOT, VRPH_DEPOT-VRPH_ABS(next_node) 02122 // First, consider i-VRPH_DEPOT 02123 k=VRPH_DEPOT; 02124 02125 ij= d[i][j]; 02126 ik= d[i][k]; 02127 jk= d[j][k]; 02128 increase=ij+jk-ik; 02129 02130 if(increase<min_increase) 02131 { 02132 min_increase=increase; 02133 edge[2]=i; 02134 edge[3]=k; 02135 } 02136 02137 if(increase<min_feasible_increase) 02138 { 02139 // Check feasibility 02140 new_route= route_num[i]; 02141 02142 if( (route[new_route].length+increase <= max_route_length) && 02143 (route[new_route].load + nodes[j].demand <= max_veh_capacity) ) 02144 { 02145 edge[0]=i; 02146 edge[1]=k; 02147 best_route=new_route; 02148 min_feasible_increase=increase; 02149 } 02150 } 02151 02152 // Now consider the edge VRPH_DEPOT-VRPH_ABS(next_node) 02153 i=VRPH_DEPOT; 02154 k=VRPH_ABS(next_node); 02155 02156 ij= d[i][j]; 02157 ik= d[i][k]; 02158 jk= d[j][k]; 02159 increase=ij+jk-ik; 02160 02161 if(increase<min_increase) 02162 { 02163 min_increase=increase; 02164 edge[2]=i; 02165 edge[3]=k; 02166 } 02167 02168 if(increase<min_feasible_increase) 02169 { 02170 new_route= route_num[k]; 02171 02172 if( (route[new_route].length+increase <= max_route_length )&& 02173 (route[new_route].load + nodes[j].demand <= max_veh_capacity ) ) 02174 { 02175 edge[0]=i; 02176 edge[1]=k; 02177 best_route=new_route; 02178 min_feasible_increase=increase; 02179 } 02180 02181 } 02182 i=k; 02183 } 02184 02185 02186 next_node= next_array[i]; 02187 02188 } 02189 // Consider the final edge 02190 k=VRPH_DEPOT; 02191 ij= d[i][j]; 02192 ik= d[i][k]; 02193 jk= d[j][k]; 02194 increase=ij+jk-ik; 02195 02196 if(increase<min_increase) 02197 { 02198 min_increase=increase; 02199 edge[2]=i; 02200 edge[3]=k; 02201 } 02202 02203 if(increase<min_feasible_increase) 02204 { 02205 // Check feasibility 02206 new_route= route_num[i]; 02207 02208 if( (route[new_route].length+increase <= max_route_length) && 02209 (route[new_route].load + nodes[j].demand <= max_veh_capacity) ) 02210 { 02211 edge[0]=i; 02212 edge[1]=k; 02213 best_route=new_route; 02214 min_feasible_increase=increase; 02215 } 02216 } 02217 02218 costs[0]=min_feasible_increase; 02219 costs[1]=min_increase; 02220 return; 02221 02222 } 02223 else 02224 { 02225 // Only search the neighbor list for positions 02226 // We already computed the VRPH_DEPOT-j-VRPH_DEPOT position. 02227 for(m=0;m<neighbor_list_size;m++) 02228 { 02229 i=nodes[j].neighbor_list[m].position; 02230 if(routed[i]) 02231 { 02232 h=VRPH_MAX(VRPH_DEPOT,pred_array[i]); 02233 increase=d[h][j]+d[j][i]-d[h][i]; 02234 02235 if(increase<min_increase) 02236 { 02237 min_increase=increase; 02238 edge[2]=h; 02239 edge[3]=i; 02240 } 02241 02242 if(increase<min_feasible_increase) 02243 { 02244 // Check feasibility 02245 if(i!=VRPH_DEPOT) 02246 new_route= route_num[i]; 02247 else 02248 new_route= route_num[k]; 02249 02250 02251 if( (route[new_route].length+increase <= max_route_length) && 02252 (route[new_route].load + nodes[j].demand <= max_veh_capacity) ) 02253 { 02254 edge[0]=h; 02255 edge[1]=i; 02256 best_route=new_route; 02257 min_feasible_increase=increase; 02258 } 02259 } 02260 02261 k=VRPH_MAX(VRPH_DEPOT,next_array[i]); 02262 increase=d[i][j]+d[j][k]-d[i][k]; 02263 02264 if(increase<min_increase) 02265 { 02266 min_increase=increase; 02267 edge[2]=i; 02268 edge[3]=k; 02269 } 02270 02271 if(increase<min_feasible_increase) 02272 { 02273 // Check feasibility 02274 if(i!=VRPH_DEPOT) 02275 new_route= route_num[i]; 02276 else 02277 new_route= route_num[k]; 02278 02279 02280 if( (route[new_route].length+increase <= max_route_length) && 02281 (route[new_route].load + nodes[j].demand <= max_veh_capacity) ) 02282 { 02283 edge[0]=i; 02284 edge[1]=k; 02285 best_route=new_route; 02286 min_feasible_increase=increase; 02287 } 02288 } 02289 } 02290 } 02291 02292 costs[0]=min_feasible_increase; 02293 costs[1]=min_increase; 02294 return; 02295 } 02296 } 02297 int VRP::inject_set(int num, int *nodelist, int rules, int attempts) 02298 { 02305 02306 if(rules != VRPH_RANDOM_SEARCH && rules != VRPH_REGRET_SEARCH) 02307 report_error("%s: invalid rules\n",__FUNCTION__); 02308 02309 int i,j,k; 02310 double best_obj=VRP_INFINITY; 02311 int *best_sol, *orderings,*best_ordering, *start_sol; 02312 int best_index=0; 02313 02314 best_sol=orderings=best_ordering=start_sol=NULL; 02315 02316 for(i=0;i<num;i++) 02317 { 02318 if(nodelist[i]==VRPH_DEPOT) 02319 { 02320 fprintf(stderr,"nodelist[%d] of %d=VRPH_DEPOT\n",i,num); 02321 report_error("%s: Cannot inject VRPH_DEPOT!\n",__FUNCTION__); 02322 } 02323 } 02324 02325 best_sol=new int[3+(this->num_nodes)+num]; 02326 start_sol=new int[3+(this->num_nodes)+num]; 02327 this->export_solution_buff(start_sol); 02328 this->import_solution_buff(start_sol); 02329 02330 if(rules==VRPH_RANDOM_SEARCH) 02331 { 02332 // Try "attempts" different random orderings to inject the nodes 02333 orderings=new int[num]; 02334 best_ordering=new int[num]; 02335 int edge[4]; 02336 double costs[2]; 02337 for(i=0;i<num;i++) 02338 orderings[i]=i; 02339 02340 for(i=0;i<attempts;i++) 02341 { 02342 // Create a random permutation 02343 random_permutation(orderings,num); 02344 for(j=0;j<num;j++) 02345 { 02346 if(nodelist[orderings[j]]==VRPH_DEPOT) 02347 { 02348 fprintf(stderr,"inject_set: Bizarre nodelist[%d]:\n",j); 02349 for(k=0;k<num;k++) 02350 fprintf(stderr,"(%d,%d) ",orderings[k],nodelist[orderings[k]]); 02351 02352 report_error("%s: Found VRPH_DEPOT in nodelist!\n",__FUNCTION__); 02353 } 02354 this->inject_node(nodelist[orderings[j]]); 02355 } 02356 02357 // Look for a new best 02358 // Look for a new best 02359 if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON) 02360 { 02361 best_index=i; 02362 best_obj=this->total_route_length; 02363 this->export_solution_buff(best_sol); 02364 memcpy(best_ordering,orderings,num*sizeof(int)); 02365 } 02366 02367 // Now return to the original solution... 02368 for(j=0;j<num;j++) 02369 this->eject_node(nodelist[orderings[j]]); 02370 } 02371 // Now revisit the best ordering found and try to improve it! 02372 for(j=0;j<num;j++) 02373 { 02374 this->find_cheapest_insertion(nodelist[best_ordering[j]],edge,costs,VRPH_USE_NEIGHBOR_LIST); 02375 this->insert_node(nodelist[best_ordering[j]],edge[0],edge[1]); 02376 } 02377 } 02378 02379 if(rules==VRPH_REGRET_SEARCH) 02380 { 02381 // Try "attempts" different random orderings to inject the nodes 02382 orderings=new int[num]; 02383 best_ordering=new int[num]; 02384 02385 int *current_list; 02386 current_list=new int[num]; 02387 int edge[4]; 02388 double costs[2]; 02389 for(i=0;i<num;i++) 02390 orderings[i]=i; 02391 02392 for(i=0;i<attempts;i++) 02393 { 02394 // Create a random permutation 02395 random_permutation(orderings,num); 02396 for(j=0;j<num;j++) 02397 current_list[j]=nodelist[orderings[j]]; 02398 02399 // It is unfortunately possible to enter some cycles here - this is a cheap 02400 // way of getting out 02401 int cycle_ctr=0; 02402 for(j=0;j<num;j++) 02403 { 02404 cycle_ctr++; 02405 if(cycle_ctr==VRPH_MAX_CYCLES) 02406 { 02407 fprintf(stderr,"Cycle encountered in REGRET SEARCH!\nReverting to original solution\n"); 02408 if(best_obj<VRP_INFINITY) 02409 { 02410 this->import_solution_buff(best_sol); 02411 delete [] best_sol; 02412 delete [] orderings; 02413 delete [] best_ordering; 02414 delete [] start_sol; 02415 delete [] current_list; 02416 return 0; 02417 } 02418 else 02419 report_error("%s: Couldn't escape cycle in REGRET search!\n",__FUNCTION__); 02420 } 02421 02422 // Find the cheapest way to insert the present node 02423 this->find_cheapest_insertion(current_list[j],edge,costs,VRPH_USE_NEIGHBOR_LIST); 02424 // Now loop over the nodes we already injected and compute their ejection costs 02425 double max_ejection_cost=-VRP_INFINITY; 02426 double current_cost; 02427 int node_to_eject=-1; 02428 for(k=0;k<j;k++) 02429 { 02430 // Make sure the node we are testing isn't involved in the cheapest insertion 02431 if(current_list[k]!=edge[0] && current_list[k]!=edge[1]) 02432 { 02433 if(( current_cost = this->ejection_cost(current_list[k])) > max_ejection_cost) 02434 { 02435 max_ejection_cost=current_cost; 02436 node_to_eject=current_list[k]; 02437 if(!routed[node_to_eject]) 02438 { 02439 fprintf(stderr,"%d NOT ROUTED!!\n",node_to_eject); 02440 report_error("%s: Error in inject_set\n",__FUNCTION__); 02441 } 02442 } 02443 } 02444 } 02445 02446 if(node_to_eject==-1 || costs[0]>= max_ejection_cost) 02447 { 02448 // Insert 02449 this->insert_node(current_list[j],edge[0],edge[1]); 02450 } 02451 else 02452 { 02453 // Eject and Insert 02454 this->eject_node(node_to_eject); 02455 this->insert_node(current_list[j],edge[0],edge[1]); 02456 current_list[j]=node_to_eject; 02457 j--; 02458 } 02459 } 02460 02461 // Look for a new best 02462 if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON) 02463 { 02464 best_index=i; 02465 best_obj=this->total_route_length; 02466 this->export_solution_buff(best_sol); 02467 memcpy(best_ordering,orderings,num*sizeof(int)); 02468 } 02469 02470 // Now return to the original solution... 02471 this->import_solution_buff(start_sol); 02472 } 02473 delete [] current_list; 02474 } 02475 02476 this->import_solution_buff(best_sol); 02477 02478 delete [] best_sol; 02479 delete [] orderings; 02480 delete [] best_ordering; 02481 delete [] start_sol; 02482 02483 02484 return best_index; 02485 } 02486 02487 void VRP::eject_neighborhood(int j, int num, int *nodelist) 02488 { 02495 02496 int i,k,cnt; 02497 int *ejected; 02498 02499 ejected=new int[this->num_nodes+1]; 02500 memset(ejected,0,(this->num_nodes+1)*sizeof(int)); 02501 02502 cnt=0; 02503 02504 // Always eject j 02505 nodelist[cnt]=j; 02506 ejected[j]=1; 02507 cnt++; 02508 02509 i=0; 02510 while(cnt<num) 02511 { 02512 // Find a node near j to eject 02513 i=(int)(lcgrand(17)*2*num); 02514 02515 // Grab it from the neighborlist 02516 i=VRPH_MIN(MAX_NEIGHBORLIST_SIZE-1,i); 02517 02518 // Get the index 02519 k=this->nodes[j].neighbor_list[i].position; 02520 if(ejected[k]==0) 02521 { 02522 if(k!=VRPH_DEPOT) 02523 { 02524 nodelist[cnt]=k; 02525 ejected[k]=1; 02526 cnt++; 02527 } 02528 } 02529 02530 } 02531 02532 // We now have cnt nodes to eject in nodelist 02533 for(i=0;i<cnt;i++) 02534 { 02535 if(nodelist[i]==VRPH_DEPOT) 02536 { 02537 fprintf(stderr,"Trying to eject VRPH_DEPOT??? j=%d\n",j); 02538 report_error("%s: Error in eject_neighborhood\n",__FUNCTION__); 02539 } 02540 02541 this->eject_node(nodelist[i]); 02542 } 02543 02544 #if 0 02545 this->verify_routes("After ejecting neighborhood\n"); 02546 #endif 02547 delete [] ejected; 02548 02549 return; 02550 } 02551 02552 02553 02554 02555 void VRP::normalize_route_numbers() 02556 { 02562 02563 int current, end, n, R, i; 02564 int current_route; 02565 int *indices; 02566 R= count_num_routes(); 02567 n= this->num_original_nodes; 02568 02569 02570 // First check to see if we are already normalized 02571 02572 // Make the arrays 1-based to avoid insanity... 02573 indices=new int[n+1]; 02574 02575 for(i=0;i<=n;i++) 02576 indices[i]=1; 02577 02578 int ctr=0; 02579 i=VRPH_ABS(next_array[VRPH_DEPOT]); 02580 while(i!=VRPH_DEPOT) 02581 { 02582 if(route_num[i] <= R && routed[i]==true) 02583 { 02584 // route_num[i] is less than R, so we have to 02585 // avoid this index. 02586 indices[route_num[i]]=0; 02587 } 02588 else 02589 ctr++; 02590 02591 i=VRPH_ABS(next_array[i]); 02592 } 02593 02594 if(ctr==0) 02595 { 02596 delete [] indices; 02597 return;// route numbers are already normalized since no route_nums were >R. 02598 } 02599 02600 // Now all of the available indices that are less than or equal to R 02601 // are marked with a 1. The unavailable indices are set to zero. 02602 02603 int next_index=1; 02604 while(indices[next_index]==0) 02605 next_index++; 02606 02607 // So next_index is the first index available for use. 02608 02609 // Get the first route and change the information 02610 current=VRPH_ABS(next_array[VRPH_DEPOT]); 02611 02612 while(current!=VRPH_DEPOT) 02613 { 02614 current_route= route_num[current]; 02615 end= route[current_route].end; 02616 02617 // Get the next available index 02618 while(indices[next_index]==0) 02619 next_index++; 02620 02621 // Modify the current index if necessary 02622 if(current_route>R) 02623 { 02624 // Update the route_num 02625 i=VRPH_ABS(next_array[VRPH_DEPOT]); 02626 while(i!=VRPH_DEPOT) 02627 { 02628 if(route_num[i]==current_route) 02629 route_num[i]=next_index; 02630 i=VRPH_ABS(next_array[i]); 02631 02632 } 02633 // Update the route information 02634 route[next_index].start = route[current_route].start; 02635 route[next_index].end = route[current_route].end; 02636 route[next_index].length = route[current_route].length; 02637 route[next_index].load = route[current_route].load; 02638 route[next_index].num_customers = route[current_route].num_customers; 02639 02640 next_index++; 02641 02642 02643 } 02644 02645 // Set current to be the first node in the next route 02646 current=VRPH_ABS(next_array[end]); 02647 02648 02649 } 02650 02651 delete [] indices; 02652 02653 return; 02654 02655 02656 02657 } 02658 02659 02660 bool VRP::create_search_neighborhood(int j, int rules) 02661 { 02666 02667 int i,k, cnt; 02668 // Define the search space 02669 02670 02671 if( rules & VRPH_USE_NEIGHBOR_LIST ) 02672 { 02673 02674 // Search only those nodes that are in the neighbor list 02675 search_size=0; 02676 cnt=0; 02677 for(i=0;i<neighbor_list_size;i++) 02678 { 02679 // Consider node k 02680 k=nodes[j].neighbor_list[i].position; 02681 if( routed[k] == true) 02682 { 02683 // The node is routed 02684 if(rules & VRPH_INTER_ROUTE_ONLY) 02685 { 02686 // Make sure k is in a different route 02687 if(route_num[k]!=route_num[j]) 02688 { 02689 search_space[cnt] = k; 02690 cnt++; 02691 02692 } 02693 } 02694 else 02695 { 02696 if(rules & VRPH_INTRA_ROUTE_ONLY) 02697 { 02698 // Make sure k is in the same route 02699 02700 if(route_num[k]==route_num[j]) 02701 { 02702 // The node is in a different route - add to the search space 02703 search_space[cnt] = k; 02704 cnt++; 02705 02706 } 02707 } 02708 else 02709 { 02710 // No intra/inter restrictions 02711 search_space[cnt] = k; 02712 cnt++; 02713 02714 02715 } 02716 } 02717 } 02718 } 02719 search_size=cnt; 02720 goto randomize; 02721 } 02722 02723 02724 // Otherwise no neighborlist 02725 if( (rules & VRPH_INTRA_ROUTE_ONLY) ) 02726 { 02727 // Search_space is just the route itself 02728 search_space[0]=VRPH_DEPOT; 02729 search_space[1]=this->route[route_num[j]].start; 02730 for(i=2;i<=this->route[route_num[j]].num_customers;i++) 02731 search_space[i]=this->next_array[search_space[i-1]]; 02732 02733 search_size=this->route[route_num[j]].num_customers+1;//add 1 for depot 02734 02735 goto randomize; 02736 02737 /* 02738 02739 search_size=0; 02740 i=0; 02741 search_space[i]=VRPH_DEPOT; 02742 i++; 02743 search_space[i]=route[route_num[j]].start; 02744 for(;;) 02745 { 02746 search_space[i]=next_array[search_space[i]]; 02747 if(search_space[i]<=VRPH_DEPOT) 02748 { 02749 search_size=i; 02750 goto randomize; 02751 } 02752 i++; 02753 }*/ 02754 02755 } 02756 02757 if( (rules & VRPH_INTER_ROUTE_ONLY) ) 02758 { 02759 // Only nodes in a different route 02760 i=0; 02761 search_space[i]=VRPH_DEPOT; 02762 search_size=1; 02763 k=VRPH_ABS(next_array[search_space[i]]); 02764 i++; 02765 02766 for(;;) 02767 { 02768 if(k==VRPH_DEPOT) 02769 // end of solution 02770 goto randomize; 02771 02772 if(route_num[k]!=route_num[j]) 02773 { 02774 search_space[i]=k; 02775 i++; 02776 search_size++; 02777 } 02778 k=VRPH_ABS(next_array[k]); 02779 } 02780 02781 } 02782 02783 02784 // Search space will consist of all positions in routes moving in the 02785 // given direction 02786 02787 if((rules & VRPH_FORWARD)) 02788 { 02789 search_size=0; 02790 02791 i=0; 02792 search_space[i]=VRPH_ABS(next_array[route[route_num[j]].end]); 02793 search_size++; 02794 02795 // First node in next route in the solution 02796 for(;;) 02797 { 02798 search_space[i+1]=VRPH_ABS(next_array[search_space[i]]); 02799 search_size++; 02800 if(search_space[i+1]==VRPH_DEPOT) 02801 goto randomize; 02802 i++; 02803 02804 } 02805 02806 } 02807 02808 if((rules & VRPH_BACKWARD)) 02809 { 02810 search_size=0; 02811 02812 i=0; 02813 search_space[i]=VRPH_ABS(pred_array[route[route_num[j]].start]); 02814 search_size++; 02815 02816 // Last node in previous route in the solution 02817 for(;;) 02818 { 02819 search_space[i+1]=VRPH_ABS(pred_array[search_space[i]]); 02820 search_size++; 02821 if(search_space[i+1]==VRPH_DEPOT) 02822 goto randomize; 02823 i++; 02824 02825 } 02826 02827 } 02828 02829 // No rules given - search space is set of all nodes 02830 search_size=0; 02831 i=0; 02832 02833 search_space[i]=VRPH_ABS(next_array[VRPH_DEPOT]); 02834 search_size++; 02835 02836 for(;;) 02837 { 02838 search_space[i+1]=VRPH_ABS(next_array[search_space[i]]); 02839 search_size++; 02840 if(search_space[i+1]==VRPH_DEPOT) 02841 goto randomize; 02842 i++; 02843 } 02844 02845 randomize: 02846 02847 // Now permute it if the rules is VRPH_RANDOMIZED 02848 if(rules & VRPH_RANDOMIZED) 02849 random_permutation(search_space, search_size); 02850 02851 return true; 02852 02853 02854 02855 } 02856 02857 02858 double VRP::insertion_cost(int u, int a, int b) 02859 { 02864 02865 // a-b is an existing edge 02866 // u is the node we are trying to insert between a and b 02867 // Returns the increased cost to the route containing a-b of 02868 // inserting u in between a and b 02869 02870 if(a==b && b==VRPH_DEPOT) 02871 return d[VRPH_DEPOT][u]+d[u][VRPH_DEPOT];// assume feasible!! 02872 02873 if(a==u || u==b || u==VRPH_DEPOT) 02874 report_error("%s: overlap or VRPH_DEPOT found\n",__FUNCTION__); 02875 02876 if(routed[a]==false || routed[b]==false) 02877 { 02878 fprintf(stderr,"%d,%d: not routed!\n",a,b); 02879 report_error("%s: Unrouted nodes in insertion_cost\n",__FUNCTION__); 02880 } 02881 02882 if( (a!=VRPH_DEPOT) && VRPH_MAX(next_array[a],VRPH_DEPOT)!=b) 02883 { 02884 fprintf(stderr,"%d,%d: not an edge!\n",a,b); 02885 report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__); 02886 } 02887 02888 if( (a==VRPH_DEPOT) && VRPH_MAX(pred_array[b],VRPH_DEPOT)!=a) 02889 { 02890 fprintf(stderr,"%d,%d: not an edge!\n",a,b); 02891 report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__); 02892 } 02893 02894 02895 int new_route; 02896 if(a==VRPH_DEPOT) 02897 new_route=route_num[b]; 02898 else 02899 new_route=route_num[a]; 02900 02901 if(nodes[u].demand+route[new_route].load>max_veh_capacity) 02902 return VRP_INFEASIBLE; 02903 02904 double increase=d[a][u] + d[u][b] - d[a][b]; 02905 if(route[new_route].length+increase>max_route_length) 02906 return VRP_INFEASIBLE; 02907 02908 return increase; 02909 02910 } 02911 02912 02913 02914 double VRP::ejection_cost(int u) 02915 { 02920 02921 if( u==VRPH_DEPOT) 02922 report_error("%s: Cannot eject the VRPH_DEPOT!\n",__FUNCTION__); 02923 02924 if(!routed[u]) 02925 return -VRP_INFINITY; 02926 02927 // t-u-v --> t-v 02928 02929 return this->d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][u]+d[u][VRPH_MAX(next_array[u],VRPH_DEPOT)] - 02930 d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][VRPH_MAX(next_array[u],VRPH_DEPOT)]; 02931 02932 02933 } 02934 void VRP::clean_route(int r, int heuristics) 02935 { 02940 02941 int i,j, r_start, r_end; 02942 double start_val, end_val, start_rlen, end_rlen; 02943 02944 OnePointMove OPM; 02945 TwoPointMove TPM; 02946 TwoOpt TO; 02947 ThreeOpt ThreeO; 02948 OrOpt OrO; 02949 02950 int rules= VRPH_INTRA_ROUTE_ONLY+VRPH_DOWNHILL+VRPH_FIRST_ACCEPT+VRPH_SAVINGS_ONLY; 02951 02952 02953 start_improving: 02954 02955 start_val = route[r].length; 02956 02957 #if CLEAN_DEBUG 02958 printf("CLEAN::start_val=%f\n",start_val); 02959 #endif 02960 02961 end_val = - VRP_INFINITY; 02962 r_start=route[r].start; 02963 r_end=route[r].end; 02964 02965 if((heuristics & ONE_POINT_MOVE)==ONE_POINT_MOVE) 02966 { 02967 // Try the ONE_POINT_MOVE for all nodes in route r; 02968 r_start=route[r].start; 02969 02970 opm_search: 02971 start_rlen=route[r].length; 02972 i=r_start; 02973 j=VRPH_MAX(next_array[i],0); 02974 while(i!=VRPH_DEPOT) 02975 { 02976 while(OPM.search(this,i,rules)); 02977 i=j; 02978 j=VRPH_MAX(next_array[j],0); 02979 } 02980 02981 end_rlen = route[r].length; 02982 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON)) 02983 goto opm_search; 02984 02985 #if CLEAN_DEBUG 02986 printf("CLEAN::OPM end_rlen=%f\n",end_rlen); 02987 #endif 02988 02989 } 02990 02991 if((heuristics & TWO_POINT_MOVE)==TWO_POINT_MOVE) 02992 { 02993 // Try the TWO_POINT_MOVE for all nodes in route r; 02994 tpm_search: 02995 r_start=route[r].start; 02996 start_rlen=route[r].length; 02997 i=r_start; 02998 j=VRPH_MAX(next_array[i],0); 02999 while(i!=VRPH_DEPOT) 03000 { 03001 while(TPM.search(this,i,rules)); 03002 i=VRPH_MAX(next_array[i],VRPH_DEPOT); 03003 } 03004 03005 end_rlen = route[r].length; 03006 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON)) 03007 goto tpm_search; 03008 03009 #if CLEAN_DEBUG 03010 printf("CLEAN::TPM end_rlen=%f\n",end_rlen); 03011 #endif 03012 03013 03014 } 03015 03016 if((heuristics & TWO_OPT)==TWO_OPT) 03017 { 03018 03019 to_search: 03020 r_start=route[r].start; 03021 start_rlen=route[r].length; 03022 i=r_start; 03023 j=VRPH_MAX(next_array[i],0); 03024 while(i!=VRPH_DEPOT) 03025 { 03026 while(TO.search(this,i,rules)); 03027 i=VRPH_MAX(next_array[i],0); 03028 } 03029 03030 end_rlen = route[r].length; 03031 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON)) 03032 goto to_search; 03033 03034 #if CLEAN_DEBUG 03035 printf("CLEAN::TO end_rlen=%f\n",end_rlen); 03036 #endif 03037 03038 03039 } 03040 03041 if((heuristics & OR_OPT)==OR_OPT) 03042 { 03043 // Use strings of length 2 and 3 03044 or_search: 03045 r_start=route[r].start; 03046 start_rlen=route[r].length; 03047 i=r_start; 03048 j=VRPH_MAX(next_array[i],0); 03049 while(i!=VRPH_DEPOT) 03050 { 03051 OrO.search(this,i,3,rules); 03052 i=VRPH_MAX(next_array[i],0); 03053 } 03054 03055 i=r_start; 03056 j=VRPH_MAX(next_array[i],0); 03057 while(i!=VRPH_DEPOT) 03058 { 03059 OrO.search(this,i,2,rules); 03060 i=VRPH_MAX(next_array[i],0); 03061 } 03062 end_rlen = route[r].length; 03063 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON)) 03064 goto or_search; 03065 03066 #if CLEAN_DEBUG 03067 printf("CLEAN::TO end_rlen=%f\n",end_rlen); 03068 #endif 03069 03070 } 03071 03072 if((heuristics & THREE_OPT)==THREE_OPT) 03073 { 03074 //show_route(r); 03075 03076 while(ThreeO.route_search(this,r,rules)) 03077 03078 end_val = route[r].length; 03079 #if CLEAN_DEBUG 03080 printf("CLEAN::3O end_val=%f\n",end_val); 03081 #endif 03082 03083 03084 } 03085 03086 end_val= route[r].length; 03087 03088 #if CLEAN_DEBUG 03089 printf("CLEAN::final end_val=%f\n",end_val); 03090 #endif 03091 03092 03093 if(VRPH_ABS(start_val-end_val)>VRPH_EPSILON) 03094 goto start_improving; 03095 else 03096 return; 03097 03098 } 03099 03100 bool VRP::before(int a, int b) 03101 { 03107 03108 int i; 03109 03110 if(a==VRPH_DEPOT || b==VRPH_DEPOT) 03111 report_error("%s: before called with VRPH_DEPOT\n",__FUNCTION__); 03112 03113 if(route_num[a]!=route_num[b]) 03114 { 03115 fprintf(stderr,"Ordering error: before called with %d and %d not in the same route!\n",a,b); 03116 report_error("%s: differnet routes\n",__FUNCTION__); 03117 } 03118 03119 if(next_array[a]==b) 03120 return true; 03121 if(next_array[b]==a) 03122 return false; 03123 03124 i=a; 03125 while(i>0 && i!=b) 03126 i=next_array[i]; 03127 03128 // At the end of this loop, if i<=0, then we're at the end of a route 03129 // and haven't encountered b ==> must have a after b in the route 03130 // If i==b, then we know that b follows a 03131 03132 if(i==b) 03133 return true; 03134 else 03135 return false; 03136 03137 } 03138 03139 03140 void VRP::update(VRPMove *M) 03141 { 03145 03146 if(M->num_affected_routes==0) // Nothing to do! 03147 return; 03148 03149 int i; 03150 03151 for(i=0; i < M->num_affected_routes; i++) 03152 { 03153 // Update length 03154 route[M->route_nums[i]].length = M->route_lens[i]; 03155 03156 // Update load 03157 route[M->route_nums[i]].load = M->route_loads[i]; 03158 03159 // Update # of customers 03160 route[M->route_nums[i]].num_customers = M->route_custs[i]; 03161 } 03162 03163 // Now update total_route_length 03164 total_route_length = M->new_total_route_length; 03165 03166 // # of routes 03167 total_number_of_routes = M->total_number_of_routes; 03168 03169 return; 03170 03171 } 03172 03173 03174 03175 03176 void VRP::compute_route_center(int r) 03177 { 03182 03183 int current_node = route[r].start; 03184 double total_x=0; 03185 double total_y=0; 03186 03187 while(current_node != VRPH_DEPOT) 03188 { 03189 total_x += nodes[current_node].x; 03190 total_y += nodes[current_node].y; 03191 current_node=VRPH_MAX(VRPH_DEPOT, next_array[current_node]); 03192 } 03193 03194 route[r].x_center = total_x / ((double) route[r].num_customers); 03195 route[r].y_center = total_y / ((double) route[r].num_customers); 03196 03197 } 03198 void VRP::find_neighboring_routes() 03199 { 03205 03206 int i,j; 03207 VRPNeighborElement **rd; 03208 03209 // Compute the route centers 03210 normalize_route_numbers(); 03211 for(i=1;i<=total_number_of_routes;i++) 03212 { 03213 compute_route_center(i); 03214 } 03215 03216 // Create the matrix 03217 rd = new VRPNeighborElement *[total_number_of_routes + 1]; 03218 rd[0]= new VRPNeighborElement[(total_number_of_routes +1) * (total_number_of_routes + 1)]; 03219 for(i=1;i<total_number_of_routes+1;i++) 03220 rd[i]=rd[i-1] + (total_number_of_routes+1); 03221 03222 // Fill the matrix with inter-route distances 03223 for(i=1;i<=this->total_number_of_routes;i++) 03224 { 03225 rd[i][0].position=VRP_INFINITY; 03226 rd[i][0].val=VRP_INFINITY; 03227 for(j=1;j<=this->total_number_of_routes;j++) 03228 { 03229 rd[i][j].position = j; 03230 rd[i][j].val = VRPDistance(VRPH_EUC_2D,route[i].x_center, route[i].y_center, 03231 route[j].x_center,route[j].y_center); 03232 } 03233 } 03234 03235 // Now sort each row in put the top MAX_NEIGHBORING_ROUTES in the route structs. 03236 03237 03238 for(i=1;i<=total_number_of_routes;i++) 03239 qsort(rd[i],total_number_of_routes+1,sizeof(VRPNeighborElement),VRPNeighborCompare); 03240 03241 03242 03243 for(i=1;i<=total_number_of_routes;i++) 03244 { 03245 for(j=0;j<MAX_NEIGHBORING_ROUTES;j++) 03246 { 03247 route[i].neighboring_routes[j]=rd[i][j+1].position; 03248 // don't include the route itself! 03249 } 03250 } 03251 03252 delete [] rd[0]; 03253 delete [] rd; 03254 03255 03256 03257 return; 03258 } 03259 03260 void VRP::capture_best_solution() 03261 { 03265 03266 if( (this->total_route_length < this->best_total_route_length) && 03267 (VRPH_ABS(this->total_route_length - this->best_total_route_length) > VRPH_EPSILON) ) 03268 { 03269 this->best_total_route_length=this->total_route_length; 03270 this->export_solution_buff(this->best_sol_buff); 03271 03272 } 03273 03274 03275 if(this->total_route_length < this->solution_wh->worst_obj || 03276 this->solution_wh->num_sols < this->solution_wh->max_size) 03277 { 03278 VRPSolution this_sol(this->num_nodes); 03279 03280 this_sol.obj=this->total_route_length; 03281 this_sol.in_IP=false; 03282 03283 // Export buffer 03284 this->export_canonical_solution_buff(this_sol.sol); 03285 03286 this->solution_wh->add_sol(&this_sol, 0); 03287 // We don't know any information to help us know where to insert 03288 } 03289 03290 return; 03291 03292 03293 } 03294 03295 void VRP::update_solution_wh() 03296 { 03300 03301 VRPSolution this_sol(this->num_nodes); 03302 03303 this_sol.obj=this->total_route_length; 03304 this_sol.in_IP=false; 03305 03306 // Export buffer 03307 this->export_canonical_solution_buff(this_sol.sol); 03308 03309 this->solution_wh->add_sol(&this_sol, 0); 03310 // We don't know any information to help us know where to insert 03311 03312 03313 return ; 03314 03315 03316 } 03317 03318 03319 void VRP::update_route(int j, VRPRoute *R) 03320 { 03326 03327 int i, current; 03328 double st=0;// for the service time 03329 03330 R->x[0]=this->nodes[VRPH_DEPOT].x; 03331 R->y[0]=this->nodes[VRPH_DEPOT].y; 03332 R->start=route[j].start; 03333 R->end=route[j].end; 03334 R->num_customers=route[j].num_customers; 03335 R->load=route[j].load; 03336 R->length=route[j].length; 03337 R->obj_val=this->total_route_length-this->total_service_time; 03338 03339 03340 if(R->start < R->end) 03341 { 03342 current=R->start; 03343 // Just output the normal ordering 03344 R->ordering[0]=current; 03345 R->x[1]=this->nodes[current].x; 03346 R->y[1]=this->nodes[current].y; 03347 st+=this->nodes[current].service_time; 03348 03349 for(i=1; i<R->num_customers; i++) 03350 { 03351 current=this->next_array[current]; 03352 st+=this->nodes[current].service_time; 03353 R->ordering[i]=current; 03354 R->x[i+1]=this->nodes[current].x; 03355 R->y[i+1]=this->nodes[current].y; 03356 } 03357 03358 R->total_service_time=st; 03359 03360 return; 03361 } 03362 03363 // Otherwise start > end - store the reverse 03364 current=R->end; 03365 R->ordering[0]=current; 03366 R->x[1]=this->nodes[current].x; 03367 R->y[1]=this->nodes[current].y; 03368 st+=this->nodes[current].service_time; 03369 03370 03371 for(i=1; i<R->num_customers; i++) 03372 { 03373 current=this->pred_array[current]; 03374 st+=this->nodes[current].service_time; 03375 R->ordering[i]=current; 03376 R->x[i+1]=this->nodes[current].x; 03377 R->y[i+1]=this->nodes[current].y; 03378 03379 } 03380 03381 R->total_service_time=st; 03382 03383 return; 03384 03385 } 03386 03387 03388 03389 double VRP::split(double p) 03390 { 03398 03399 if(p>.5) 03400 report_error("%s: p must be less than .5\n"); 03401 03402 int i,j,k; 03403 double beta; 03404 struct double_int *thetas; 03405 03406 thetas=new double_int[this->num_nodes]; 03407 03408 for(i=1;i<=this->num_nodes;i++) 03409 { 03410 thetas[i-1].k=i; 03411 thetas[i-1].d=this->nodes[i].theta; 03412 } 03413 03414 // Sort the list of thetas 03415 qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare); 03416 03417 // Select a random beta in (min_theta, max_theta ) and see how the split is 03418 // Stop when we get a satisfactory split 03419 03420 03421 for(;;) 03422 { 03423 beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta); 03424 03425 // Now find out how many nodes we get in each half 03426 k=0; 03427 for(j=0;j<this->num_nodes;j++) 03428 { 03429 // Count the nodes above the line 03430 if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x)) 03431 k++; 03432 } 03433 03434 if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes)))) 03435 break; // The split is good 03436 } 03437 03438 // Eject them - note that we assume to have begun with a full solution... 03439 for(i=1;i<=this->num_original_nodes;i++) 03440 { 03441 if(routed[i]) 03442 { 03443 if(this->nodes[i].y >= tan(beta)*(this->nodes[i].x)) 03444 this->eject_node(i); 03445 } 03446 } 03447 03448 03449 delete [] thetas; 03450 return beta; 03451 } 03452 03453 03454 int VRP::split_routes(double p, int **ejected_routes, double *t) 03455 { 03466 03467 if(p>.5) 03468 report_error("%s: p must be less than .5\n"); 03469 03470 int i,j,k; 03471 double beta; 03472 struct double_int *thetas; 03473 03474 thetas=new double_int[this->num_nodes]; 03475 03476 for(i=1;i<=this->num_nodes;i++) 03477 { 03478 thetas[i-1].k=i; 03479 thetas[i-1].d=this->nodes[i].theta; 03480 } 03481 03482 // Sort the list of thetas 03483 qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare); 03484 03485 // Select a random beta in (min_theta, max_theta ) and see how the split is 03486 // Stop when we get a satisfactory split 03487 03488 03489 for(;;) 03490 { 03491 beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta); 03492 03493 // Now find out how many nodes we get in each half 03494 k=0; 03495 for(j=0;j<this->num_nodes;j++) 03496 { 03497 // Count the nodes above the line 03498 if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x)) 03499 k++; 03500 } 03501 03502 if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes)))) 03503 { 03504 break; // The split is good 03505 } 03506 } 03507 03508 // Eject the routes that do not have any nodes in the larger portion 03509 03510 int num_ejected=0; 03511 for(i=1;i<=this->total_number_of_routes;i++) 03512 { 03513 int start=this->route[i].start; 03514 int current=start; 03515 bool will_eject=true; 03516 03517 while(current!=VRPH_DEPOT) 03518 { 03519 03520 if(this->nodes[current].y <= tan(beta)*(this->nodes[current].x)) 03521 { 03522 // We have at least one node inside the large part 03523 will_eject=false; 03524 break; 03525 } 03526 03527 current=VRPH_MAX(VRPH_DEPOT,this->next_array[current]); 03528 } 03529 if(will_eject) 03530 { 03531 // Eject route i from the solution - place the ejected nodes in 03532 // the ejected_routes[][] array 03533 this->eject_route(i, ejected_routes[num_ejected]); 03534 num_ejected++; 03535 03536 } 03537 03538 } 03539 03540 delete [] thetas; 03541 03542 *t=beta; 03543 return num_ejected; 03544 } 03545 03546 03547 void VRP::fix_edge(int start, int end) 03548 { 03556 03557 this->fixed[start][end]=true; 03558 this->fixed[end][start]=true; 03559 03560 // Handle the dummy_node 03561 if(start==VRPH_DEPOT) 03562 { 03563 this->fixed[dummy_index][start]=true; 03564 this->fixed[start][dummy_index]=true; 03565 } 03566 03567 if(end==VRPH_DEPOT) 03568 { 03569 this->fixed[dummy_index][end]=true; 03570 this->fixed[end][dummy_index]=true; 03571 } 03572 03573 03574 } 03575 03576 void VRP::unfix_edge(int start, int end) 03577 { 03581 03582 if(this->fixed[start][end]) 03583 report_error("%s: Edge %d-%d is not already fixed!\n",__FUNCTION__,start,end); 03584 03585 this->fixed[start][end]=false; 03586 this->fixed[end][start]=false; 03587 03588 // Handle the dummy_node 03589 if(start==VRPH_DEPOT) 03590 { 03591 this->fixed[dummy_index][start]=false; 03592 this->fixed[start][dummy_index]=false; 03593 } 03594 03595 if(end==VRPH_DEPOT) 03596 { 03597 this->fixed[dummy_index][end]=false; 03598 this->fixed[end][dummy_index]=false; 03599 } 03600 03601 03602 } 03603 03604 void VRP::unfix_all() 03605 { 03609 03610 for(int i=0;i<=matrix_size;i++) 03611 for(int j=0;j<=matrix_size;j++) 03612 fixed[i][j]=false; 03613 } 03614 03615 void VRP::fix_string(int *node_string, int k) 03616 { 03623 03624 03625 for(int i=0;i<k-1;i++) 03626 this->fix_edge(node_string[i], node_string[i+1]); 03627 03628 03629 } 03630 03631 03632 03633 03634 03635 int VRP::read_fixed_edges(const char *filename) 03636 { 03649 03650 03651 int i, a, b, k; 03652 FILE *in; 03653 03654 if( (in=fopen(filename,"r"))==NULL) 03655 { 03656 fprintf(stderr,"Error opening %s for reading\n",filename); 03657 report_error("%s: Bad file in read_fixed_edges\n",__FUNCTION__); 03658 } 03659 03660 fscanf(in,"%d\n",&k); 03661 03662 for(i=0;i<k;i++) 03663 { 03664 fscanf(in,"%d %d\n",&a, &b); 03665 if(a<0 || b<0 || a>this->num_original_nodes || b>this->num_original_nodes) 03666 { 03667 fprintf(stderr,"Tried to fix edge %d-%d\n",a,b); 03668 report_error("%s: Error in read_fixed_edges\n",__FUNCTION__); 03669 } 03670 this->fix_edge(a,b); 03671 } 03672 03673 // Return the number of fixed edges read in 03674 return k; 03675 } 03676 void VRP::list_fixed_edges(int *fixed_list) 03677 { 03684 03685 int pos=0; 03686 int current=VRPH_DEPOT; 03687 int next=VRPH_ABS(next_array[current]); 03688 03689 while(next!=VRPH_DEPOT) 03690 { 03691 if(fixed[current][next]) 03692 { 03693 fixed_list[pos]=current; 03694 fixed_list[pos+1]=next; 03695 pos+=2; 03696 } 03697 03698 current=next; 03699 next=next_array[current]; 03700 03701 if(next<0) 03702 { 03703 if(fixed[current][VRPH_DEPOT]) 03704 { 03705 fixed_list[pos]=current; 03706 fixed_list[pos+1]=VRPH_DEPOT; 03707 pos+=2; 03708 } 03709 03710 if(fixed[VRPH_DEPOT][-next]) 03711 { 03712 fixed_list[pos]=VRPH_DEPOT; 03713 fixed_list[pos+1]=-next; 03714 pos+=2; 03715 } 03716 03717 next=-next; 03718 } 03719 } 03720 03721 03722 03723 } 03724 void VRP::perturb_locations(double c) 03725 { 03732 03733 int i, pre, post; 03734 double v,theta; 03735 03736 // First export the solution buffer 03737 this->export_solution_buff(this->current_sol_buff); 03738 03739 for(i=1;i<=this->num_nodes;i++) 03740 { 03741 pre= VRPH_MAX(VRPH_DEPOT,this->pred_array[i]); 03742 post=VRPH_MAX(VRPH_DEPOT,this->next_array[i]); 03743 03744 v= (this->d[pre][i]+this->d[i][post]) - (this->nodes[i].service_time + this->nodes[post].service_time); 03745 v= v*c; 03746 03747 // Now alter the location of node i by shifting v units in a random direction theta 03748 theta= VRPH_PI*lcgrand(8); 03749 03750 this->nodes[i].x += v*cos(theta); 03751 this->nodes[i].y += v*sin(theta); 03752 03753 } 03754 03755 // Now re-create the distance matrix 03756 03757 this->create_distance_matrix(this->edge_weight_type); 03758 03759 // Now re-import the solution with the new distances 03760 this->import_solution_buff(this->current_sol_buff); 03761 03762 return; 03763 03764 } 03765 void VRP::add_route(int *route_buff) 03766 { 03772 03773 // We will do this by exporting the current solution to a buffer, appending the 03774 // new route to this buffer, and then importing the new resulting solution. 03775 // Probably not the fastest way, but we have many fields to update when adding a 03776 // new route! 03777 03778 int *temp_buff; 03779 03780 this->verify_routes("Before adding route\n"); 03781 03782 temp_buff=new int[this->num_original_nodes+2]; 03783 this->export_solution_buff(temp_buff); 03784 03785 03786 int old_num=this->num_nodes; 03787 int i=0; 03788 while(route_buff[i]!=-1) 03789 { 03790 temp_buff[old_num+1+i]=route_buff[i]; 03791 if(i==0) 03792 temp_buff[old_num+1+i]=-temp_buff[old_num+1+i]; 03793 03794 temp_buff[0]++; 03795 i++; 03796 } 03797 03798 temp_buff[old_num+1+i]=VRPH_DEPOT; 03799 03800 // Now import the new solution 03801 this->import_solution_buff(temp_buff); 03802 03803 this->verify_routes("After adding route\n"); 03804 03805 delete [] temp_buff; 03806 03807 } 03808 void VRP::append_route(int *sol_buff, int *route_buff) 03809 { 03816 03817 int i,j,current_num; 03818 03819 current_num=sol_buff[0]; // The # of nodes in the sol_buff solution 03820 03821 // Now find out how many nodes are in the route_buff - ends in a -1; 03822 03823 j=0; 03824 while(route_buff[j]!=-1) 03825 j++; 03826 03827 // The route has j nodes - we assume they are all different from those in the solution! 03828 // Increment the # of nodes in sol_buff; 03829 sol_buff[0]+=j; 03830 03831 sol_buff[current_num+1]=-route_buff[0]; 03832 for(i=1;i<j;i++) 03833 sol_buff[current_num+1+i]=route_buff[i]; 03834 03835 03836 // End at the VRPH_DEPOT; 03837 sol_buff[current_num+1+j]=VRPH_DEPOT; 03838 03839 return; 03840 03841 03842 } 03843 03844 int VRP::intersect_solutions(int *new_sol, int **route_list, int *sol1, int *sol2, int min_routes) 03845 { 03856 03857 int i,j,k,m,num_routes; 03858 int *rnums; 03859 03860 rnums=new int[this->num_original_nodes+1]; // Plenty big... 03861 03862 // Find the routes in common 03863 j = this->find_common_routes(sol1, sol2, rnums); 03864 03865 if(j==0) 03866 { 03867 // No routes in common - just put sol1 into new_sol 03868 memcpy(new_sol,sol1, (this->num_original_nodes+2)*sizeof(int)); 03869 delete[] rnums; 03870 return 0; 03871 } 03872 03873 // Otherwise we have some routes in common - copy them to the routes[] array 03874 // Import sol1 to do this 03875 this->import_solution_buff(sol1); 03876 num_routes=this->total_number_of_routes; 03877 k=0; 03878 for(i=0;i<j;i++) 03879 { 03880 // Copy route i from sol1 into the route_list[] array 03881 route_list[i][0]=this->route[rnums[i]].start; 03882 03883 m=0; 03884 while(route_list[i][m]!=this->route[rnums[i]].end) 03885 { 03886 m++; 03887 route_list[i][m]=this->next_array[route_list[i][m-1]]; 03888 } 03889 // Add a -1 to indicate the end of the route 03890 m++; 03891 route_list[i][m]=-1; 03892 03893 k++; 03894 num_routes--; 03895 if(num_routes==min_routes) 03896 break; 03897 } 03898 03899 // We actually copied k of the j routes in the intersection 03900 // k may be less than j 03901 03902 // Now eject the routes from sol1 (already imported) 03903 int *junk; 03904 junk=new int[this->num_original_nodes]; 03905 // won't use this 03906 for(i=0;i<k;i++) 03907 this->eject_route(rnums[i], junk); 03908 03909 // Now export the solution after ejections to new_sol 03910 this->export_canonical_solution_buff(new_sol); 03911 03912 // Now import the new (possibly partial) solution 03913 this->import_solution_buff(new_sol); 03914 03915 delete [] rnums; 03916 delete [] junk; 03917 03918 return k; 03919 03920 } 03921 bool VRP::osman_insert(int k, double alpha) 03922 { 03930 03931 int i,j,l,ll,m,mm, best_l, best_m; 03932 double savings, best_savings, ik,kj,ij; 03933 Postsert postsert; 03934 Presert presert; 03935 VRPMove M; 03936 03937 if(k==VRPH_DEPOT) 03938 return false; 03939 03940 i=VRPH_MAX(this->pred_array[k],VRPH_DEPOT); 03941 j=VRPH_MAX(this->next_array[k],VRPH_DEPOT); 03942 03943 03944 ik=this->d[i][k]; 03945 kj=this->d[k][j]; 03946 ij=this->d[i][j]; 03947 03948 best_savings=VRP_INFINITY; 03949 best_l=-1; 03950 best_m=-1; 03951 03952 l=VRPH_DEPOT; 03953 m=abs(this->next_array[l]); 03954 while(m!=VRPH_DEPOT) 03955 { 03956 if(m>0) 03957 { 03958 savings = (ik + kj - this->d[l][m]) - alpha*(this->d[l][k] + this->d[k][m] - ij); 03959 if(savings<best_savings && l!=i && m!=j) 03960 { 03961 // Now check feasibility-put k before m 03962 if(presert.evaluate(this,k,m,&M)==true) 03963 { 03964 best_savings=savings; 03965 best_l=l; 03966 best_m=m; 03967 } 03968 } 03969 } 03970 else 03971 { 03972 // m<0 indicating the end of a route 03973 03974 // consider the edge l-VRPH_DEPOT 03975 mm=VRPH_DEPOT; 03976 savings = (ik + kj - this->d[l][mm]) - alpha*(this->d[l][k] + this->d[k][mm] - ij); 03977 if(savings<best_savings && l!=i && mm!=j) 03978 { 03979 // put k after l 03980 if(postsert.evaluate(this,k,l,&M)==true) 03981 { 03982 best_savings=savings; 03983 best_l=l; 03984 best_m=mm; 03985 } 03986 } 03987 03988 // Now consider the edge VRPH_DEPOT-abs(m) 03989 ll=VRPH_DEPOT; 03990 m=abs(m); 03991 savings = (ik + kj - this->d[ll][m]) - alpha*(this->d[ll][k] + this->d[k][m] - ij); 03992 if(savings<best_savings && ll!=i && m!=j) 03993 { 03994 // put k before m 03995 if(presert.evaluate(this,k,m,&M)==true) 03996 { 03997 best_savings=savings; 03998 best_l=ll; 03999 best_m=m; 04000 } 04001 } 04002 } 04003 // Now advance to the next edge 04004 l=m; 04005 m=this->next_array[l]; 04006 } 04007 04008 // We now have the best location to move node k by 04009 // placing in between best_l and best_m; 04010 04011 if(best_savings==VRP_INFINITY) 04012 { 04013 // No locations found 04014 return false; 04015 } 04016 04017 if(best_l!=VRPH_DEPOT) 04018 { 04019 if(postsert.move(this,k,best_l)==false) 04020 report_error("%s: postsert.move is false\n",__FUNCTION__); 04021 } 04022 else 04023 { 04024 if(presert.move(this,k,best_m)==false) 04025 report_error("%s: presert.move is false\n",__FUNCTION__); 04026 } 04027 04028 return true; 04029 04030 } 04031 04032 int VRP::osman_perturb(int num, double alpha) 04033 { 04039 04040 int tot=0; 04041 int k; 04042 int num_attempts=0; 04043 04044 while(tot<num) 04045 { 04046 // Select a random node to move 04047 k=VRPH_MAX((int)(this->num_original_nodes * lcgrand(10)),1);//to avoid inserting VRPH_DEPOT 04048 04049 if(routed[k]) 04050 { 04051 04052 // Try to insert 04053 if(this->osman_insert(k,alpha)==true) 04054 tot++; 04055 04056 num_attempts++; 04057 if(num_attempts>2*this->num_original_nodes) 04058 return tot;// Only moved tot nodes 04059 } 04060 } 04061 04062 // We were able to move num nodes 04063 04064 return num; 04065 04066 04067 } 04068 04069 bool VRP::check_fixed_edges(const char *message) 04070 { 04076 04077 04078 04079 int i,j; 04080 04081 for(i=0;i<=this->num_original_nodes;i++) 04082 { 04083 for(j=0;j<=this->num_original_nodes;j++) 04084 { 04085 if(this->fixed[i][j]) 04086 { 04087 // Make sure i-j or j-i exists 04088 if(i!=VRPH_DEPOT) 04089 { 04090 if(VRPH_MAX(this->next_array[i],VRPH_DEPOT)!=j && VRPH_MAX(this->pred_array[i],VRPH_DEPOT)!=j) 04091 { 04092 fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j); 04093 fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[i],VRPH_DEPOT),i, 04094 VRPH_MAX(this->next_array[i],VRPH_DEPOT)); 04095 fprintf(stderr,"%s", message); 04096 04097 if(this->fixed[j][i]) 04098 fprintf(stderr,"%d-%d also fixed\n",j,i); 04099 else 04100 fprintf(stderr,"%d-%d NOT fixed!\n",j,i); 04101 04102 return false; 04103 } 04104 } 04105 04106 // Make sure i-j or j-i exists 04107 if(j!=VRPH_DEPOT) 04108 { 04109 if(VRPH_MAX(this->next_array[j],VRPH_DEPOT)!=i && VRPH_MAX(this->pred_array[j],VRPH_DEPOT)!=i) 04110 { 04111 fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j); 04112 fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[j],VRPH_DEPOT),j, 04113 VRPH_MAX(this->next_array[j],VRPH_DEPOT)); 04114 fprintf(stderr,"%s", message); 04115 if(this->fixed[j][i]) 04116 fprintf(stderr,"%d-%d also fixed\n",j,i); 04117 else 04118 fprintf(stderr,"%d-%d NOT fixed!\n",j,i); 04119 return false; 04120 04121 } 04122 } 04123 } 04124 04125 } 04126 04127 } 04128 04129 return true; 04130 } 04131 int VRP::find_common_routes(int *sol1, int *sol2, int *route_nums) 04132 { 04138 04139 int *h11, *h12, *h21, *h22; // to store the hash values 04140 double *L1, *L2; 04141 int i, j, cnt, r1, r2; // to store the # of routes. 04142 VRPRoute rte(num_nodes); 04143 04144 // First, import sol1 04145 this->import_solution_buff(sol1); 04146 04147 r1=this->total_number_of_routes; 04148 h11=new int[r1+1]; h12=new int[r1+1];// 1-based arrays 04149 L1=new double[r1+1]; 04150 for(i=1;i<=r1;i++) 04151 { 04152 // Import route i into rte 04153 this->update_route(i,&rte); 04154 // Hash the route 04155 h11[i]=rte.hash(SALT_1); 04156 h12[i]=rte.hash(SALT_2); 04157 L1[i]=rte.length; 04158 } 04159 04160 // Now import sol2 04161 this->import_solution_buff(sol2); 04162 04163 r2=this->total_number_of_routes; 04164 h21=new int[r2+1]; h22=new int[r2+1];// 1-based arrays 04165 L2=new double[r2+1]; 04166 for(i=1;i<=r2;i++) 04167 { 04168 this->update_route(i,&rte); 04169 h21[i]=rte.hash(SALT_1); 04170 h22[i]=rte.hash(SALT_2); 04171 L2[i]=rte.length; 04172 } 04173 04174 04175 04176 // Now compare the hash values 04177 cnt=0; 04178 for(i=1;i<=r1;i++) 04179 { 04180 for(j=1;j<=r2;j++) 04181 { 04182 if(h11[i]==h21[j] && h12[i]==h22[j] && VRPH_ABS(L1[i]-L2[j])<VRPH_EPSILON) 04183 { 04184 // route i in sol1 is the same as route j in sol2 04185 route_nums[cnt]=i; 04186 cnt++; 04187 } 04188 } 04189 } 04190 04191 04192 delete [] h11; delete [] h12; delete [] h21; delete [] h22; 04193 delete [] L1; delete [] L2; 04194 return cnt; 04195 04196 } 04197 04198 void VRP::set_daily_demands(int day) 04199 { 04205 04206 int i; 04207 04208 if(day>0) 04209 { 04210 for(i=0;i<=this->num_original_nodes;i++) 04211 { 04212 if(this->nodes[i].daily_demands[day]>=0) 04213 this->nodes[i].demand=this->nodes[i].daily_demands[day]; 04214 else 04215 this->nodes[i].demand=-1; 04216 04217 } 04218 } 04219 else 04220 { 04221 // Day is 0 - use the mean demand 04222 for(i=0;i<=this->num_original_nodes;i++) 04223 { 04224 int mean_demand=0; 04225 int k=0; 04226 for(int j=1;j<=this->num_days;j++) 04227 { 04228 if(this->nodes[i].daily_demands[j]>=0) 04229 { 04230 mean_demand+=this->nodes[i].daily_demands[j]; 04231 k++; 04232 } 04233 } 04234 if(k>0) 04235 this->nodes[i].demand=(int)((double)mean_demand/(double)k); 04236 else 04237 this->nodes[i].demand=-1; 04238 04239 } 04240 04241 } 04242 } 04243 04244 void VRP::set_daily_service_times(int day) 04245 { 04250 04251 int i; 04252 04253 for(i=0;i<=this->num_original_nodes;i++) 04254 { 04255 this->nodes[i].service_time=this->nodes[i].daily_service_times[day]; 04256 } 04257 04258 // We also have to recompute the distance matrix if the service times are not identical 04259 // across the days 04260 this->create_distance_matrix(this->edge_weight_type); 04261 } 04262 04263 04264 void VRP::update_arrival_times() 04265 { 04269 04270 int routenum, next, current; 04271 double t; 04272 04273 // Set the arrival times to -1 so that the only ones with positive 04274 // arrival_time values are those that are actually visited 04275 for(int i=1;i<=this->num_original_nodes;i++) 04276 this->nodes[i].arrival_time=-1; 04277 04278 04279 for(int i=1;i<=this->num_original_nodes;i++) 04280 { 04281 04282 if(routed[i]) 04283 { 04284 // Get the route number 04285 routenum=this->route_num[i]; 04286 04287 // Start at the beginning of this route 04288 current=this->route[routenum].start; 04289 t=this->d[VRPH_DEPOT][current]; 04290 04291 while(current!=i) 04292 { 04293 next=VRPH_MAX(VRPH_DEPOT, this->next_array[current]); 04294 t+=this->d[current][next]; 04295 current=next; 04296 04297 } 04298 04299 // Now subtract the service time at t as this is included in the distance 04300 // matrix 04301 t-=this->nodes[i].service_time; 04302 this->nodes[i].arrival_time=t; 04303 } 04304 } 04305 04306 } 04307 04308 bool VRP::check_tabu_status(VRPMove *M, int *old_sol) 04309 { 04317 04318 04319 #if VRPH_TABU_DEBUG 04320 printf("Checking tabu status\n"); 04321 #endif 04322 04323 04324 int i,j; 04325 04326 // We will always accept a move that reduces the # of routes 04327 if(M->num_affected_routes>1) 04328 { 04329 for(i=0;i<M->num_affected_routes;i++) 04330 { 04331 if(M->route_custs[i]==0) 04332 return true; 04333 } 04334 } 04335 04336 VRPRoute r(this->num_nodes); 04337 int num_tabu_routes=0; 04338 04339 #if VRPH_TABU_DEBUG 04340 printf("Checking tabu status of %d routes\n",M->num_affected_routes); 04341 #endif 04342 04343 for(i=0;i<M->num_affected_routes;i++) 04344 { 04345 this->update_route(M->route_nums[i],&r); 04346 r.hash_val = r.hash(SALT_1); 04347 r.hash_val2 = r.hash(SALT_2); 04348 04349 for(j=0;j<this->tabu_list->num_entries;j++) 04350 { 04351 if(r.hash_val==this->tabu_list->hash_vals1[j] && r.hash_val2==this->tabu_list->hash_vals2[j]) 04352 // The move is tabu! 04353 num_tabu_routes++; 04354 } 04355 04356 } 04357 04358 04359 if(num_tabu_routes>0) 04360 { 04361 // At least one tabu route 04362 this->import_solution_buff(old_sol); 04363 #if VRPH_TABU_DEBUG 04364 printf("Move is Tabu! Reverting to old solution\n"); 04365 #endif 04366 return false; 04367 } 04368 04369 #if VRPH_TABU_DEBUG 04370 printf("Move is not Tabu. Updating list of %d routes\n",M->num_affected_routes); 04371 #endif 04372 04373 // The move is not tabu - update the tabu list 04374 for(i=0;i<M->num_affected_routes;i++) 04375 { 04376 this->update_route(M->route_nums[i],&r); 04377 04378 this->tabu_list->update_list(&r); 04379 } 04380 04381 return true; 04382 04383 } 04384 04385 04386 void VRP::print_stats() 04387 { 04392 04393 printf(" (Moves, Evaluations)\n"); 04394 printf(" One Point Move: (%010d, %010d)\n", 04395 this->num_moves[ONE_POINT_MOVE_INDEX], this->num_evaluations[ONE_POINT_MOVE_INDEX]); 04396 printf(" Two Point Move: (%010d, %010d)\n", 04397 this->num_moves[TWO_POINT_MOVE_INDEX], this->num_evaluations[TWO_POINT_MOVE_INDEX]); 04398 printf(" Three Point Move: (%010d, %010d)\n", 04399 this->num_moves[THREE_POINT_MOVE_INDEX], this->num_evaluations[THREE_POINT_MOVE_INDEX]); 04400 printf(" Two-opt Move: (%010d, %010d)\n", 04401 this->num_moves[TWO_OPT_INDEX], this->num_evaluations[TWO_OPT_INDEX]); 04402 printf(" Three-opt Move: (%010d, %010d)\n", 04403 this->num_moves[THREE_OPT_INDEX], this->num_evaluations[THREE_OPT_INDEX]); 04404 printf(" Or-opt Move: (%010d, %010d)\n", 04405 this->num_moves[OR_OPT_INDEX], this->num_evaluations[OR_OPT_INDEX]); 04406 printf("Cross-Exchange Move: (%010d, %010d)\n\n", 04407 this->num_moves[CROSS_EXCHANGE_INDEX], this->num_evaluations[CROSS_EXCHANGE_INDEX]); 04408 04409 return; 04410 04411 } 04412 04413 04414 void VRP::reset() 04415 { 04419 04420 this->solution_wh->liquidate(); 04421 for(int j=1;j<=this->num_original_nodes;j++) 04422 this->routed[j]=false; 04423 } 04424