00001
00002
00003
00004
00005
00006
00007
00008
00009
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\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
00047
00048 symmetric=true;
00049
00050
00051 forbid_tiny_moves=true;
00052
00053
00054 d=NULL;
00055
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
00068
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
00078
00079
00080 record = 0.0;
00081 deviation = VRPH_DEFAULT_DEVIATION;
00082
00083
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
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
00101
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
00130
00131 forbid_tiny_moves=true;
00132
00133
00134 d=NULL;
00135
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
00148
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
00159
00160
00161 record = 0.0;
00162 deviation = VRPH_DEFAULT_DEVIATION;
00163
00164
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
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
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
00189
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
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
00327 memcpy(this->best_sol_buff, W->best_sol_buff, (sizeof(int))*(W->num_nodes+2));
00328
00329
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;
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;
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
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
00382
00383 if(type==VRPH_EXPLICIT)
00384
00385 return;
00386
00387 n= this->num_original_nodes;
00388 k=0;
00389
00390
00391
00392
00393 for(i=0;i<=n+1;i++)
00394 {
00395 for(j=0;j<=n+1;j++)
00396 this->d[i][j]=VRPDistance(type, this->nodes[i].x,this->nodes[i].y,this->nodes[j].x,
00397 this->nodes[j].y) + this->nodes[j].service_time ;
00398
00399 }
00400
00401
00402 return;
00403
00404 }
00405
00406 void VRP::create_neighbor_lists(int nsize)
00407 {
00412
00413 if(nsize>num_nodes )
00414 {
00415 fprintf(stderr,"Requested neighbor list size is greater than num_nodes!\n%d>%d\n",
00416 nsize,num_nodes);
00417 report_error("%s: Neighbor list error!!\n",__FUNCTION__);
00418 }
00419
00420 if(nsize>MAX_NEIGHBORLIST_SIZE )
00421 {
00422 fprintf(stderr,"Requested neighbor list size is greater than MAX_NEIGHBORLIST_SIZE!\n%d>%d\n",
00423 nsize,MAX_NEIGHBORLIST_SIZE);
00424 report_error("%s: Neighbor list error!!\n",__FUNCTION__);
00425 }
00426
00427
00428 int i,j,n,b;
00429 int k;
00430 double dd;
00431 double max;
00432 int maxpos;
00433 VRPNeighborElement *NList;
00434
00435 n= num_nodes;
00436
00437
00438 neighbor_list_size=nsize;
00439
00440 NList=new VRPNeighborElement[nsize];
00441
00442
00443 max=0.0;
00444
00445 maxpos=0;
00446
00447 for(i=1;i<=nsize;i++)
00448 {
00449 NList[i-1].val=d[VRPH_DEPOT][i];
00450
00451 NList[i-1].position = i;
00452
00453
00454 if(NList[i-1].val > max)
00455 {
00456 max=NList[i-1].val;
00457 maxpos=i-1;
00458 }
00459
00460 }
00461
00462
00463 for(i=nsize+1;i<=n;i++)
00464 {
00465 if(d[VRPH_DEPOT][i]<max)
00466 {
00467
00468 NList[maxpos].val=d[VRPH_DEPOT][i];
00469
00470 NList[maxpos].position=i;
00471 }
00472
00473 max=0.0;
00474 for(b=0;b<nsize;b++)
00475 {
00476 if(NList[b].val>max)
00477 {
00478 maxpos=b;
00479 max=NList[b].val;
00480 }
00481 }
00482 }
00483
00484 qsort(NList,nsize,sizeof(VRPNeighborElement), VRPNeighborCompare);
00485 i=0;
00486 for(b=0;b<nsize;b++)
00487 {
00488 nodes[VRPH_DEPOT].neighbor_list[b].val=NList[b].val;
00489 nodes[VRPH_DEPOT].neighbor_list[b].position=NList[b].position;
00490
00491 # if NEIGHBOR_DEBUG
00492 printf("(%d,%d,%f) \n",VRPH_DEPOT,nodes[VRPH_DEPOT].neighbor_list[b].position,
00493 nodes[VRPH_DEPOT].neighbor_list[b].val);
00494 #endif
00495
00496
00497 }
00498
00499 # if NEIGHBOR_DEBUG
00500 printf("\n\n");
00501 #endif
00502
00503
00504 k=0;
00505
00506 for(i=1;i<=n;i++)
00507 {
00508
00509
00510 for(j=0;j<nsize;j++)
00511 {
00512 NList[j].position=VRP_INFINITY;
00513 NList[j].val=VRP_INFINITY;
00514 }
00515
00516
00517
00518 NList[0].position=VRPH_DEPOT;
00519 NList[0].val=d[i][VRPH_DEPOT];
00520
00521
00522 for(j=1;j<i;j++)
00523 {
00524 #if NEIGHBOR_DEBUG
00525 printf("Loop 1; i=%d; in j loop(%d)\n",i,j);
00526 #endif
00527
00528 dd=d[i][j];
00529
00530 if(j<nsize)
00531 {
00532
00533 NList[j].val=dd;
00534 NList[j].position=j;
00535
00536 if(NList[j].val>max)
00537 {
00538 max=NList[j].val;
00539 maxpos=j;
00540
00541 }
00542 }
00543 else
00544 {
00545
00546 if(dd<max)
00547 {
00548 NList[maxpos].val=dd;
00549 NList[maxpos].position=j;
00550
00551 max=0.0;
00552 for(b=0;b<nsize;b++)
00553 {
00554 if(NList[b].val>max)
00555 {
00556 max=NList[b].val;
00557 maxpos=b;
00558 }
00559 }
00560 }
00561 }
00562
00563 }
00564
00565
00566 for(j=i+1;j<=n;j++)
00567 {
00568 #if NEIGHBOR_DEBUG
00569 printf("Loop 2; i=%d; in j loop(%d)\n",i,j);
00570 #endif
00571
00572 dd=d[i][j];
00573
00574 if(j<=nsize)
00575 {
00576
00577 NList[j-1].val=dd;
00578 NList[j-1].position=j;
00579
00580 if(NList[j-1].val>max)
00581 {
00582 max=NList[j-1].val;
00583 maxpos=j-1;
00584
00585 }
00586 }
00587 else
00588 {
00589 if(dd<max)
00590 {
00591
00592 NList[maxpos].val=dd;
00593 NList[maxpos].position=j;
00594
00595 max=0.0;
00596 for(b=0;b<nsize;b++)
00597 {
00598 if(NList[b].val>max)
00599 {
00600 max=NList[b].val;
00601 maxpos=b;
00602 }
00603 }
00604 }
00605 }
00606
00607 }
00608
00609
00610 #if NEIGHBOR_DEBUG
00611 printf("Node %d neighbor list before sorting\n",i);
00612 for(b=0;b<nsize;b++)
00613 {
00614 printf("NList[%d]=%d,%f[%f]\n",b,NList[b].position,NList[b].val,d[i][NList[b].position]);
00615 }
00616 #endif
00617
00618
00619 qsort (NList, nsize, sizeof(VRPNeighborElement), VRPNeighborCompare);
00620
00621 for(b=0;b<nsize;b++)
00622 {
00623 nodes[i].neighbor_list[b].position=NList[b].position;
00624 nodes[i].neighbor_list[b].val=NList[b].val;
00625 if(i== nodes[i].neighbor_list[b].position)
00626 {
00627 fprintf(stderr,"ERROR:: Node %d is in it's own neighbor list!!\n", i);
00628 fprintf(stderr,"i=%d; %d[%d],%f[%f])\n",i,nodes[i].neighbor_list[b].position,
00629 NList[b].position,nodes[i].neighbor_list[b].val,NList[b].val);
00630
00631 report_error("%s: Error creating neighbor lists\n",__FUNCTION__);
00632 }
00633 #if NEIGHBOR_DEBUG
00634
00635 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]);
00636
00637 #endif
00638 }
00639
00640
00641
00642
00643
00644 }
00645
00646 delete [] NList;
00647
00648 return;
00649
00650 }
00651
00652
00653 bool VRP::check_feasibility(VRPViolation *VV)
00654 {
00660
00661 int i;
00662 bool is_feasible=true;
00663
00664
00665 VV->capacity_violation=-VRP_INFINITY;
00666 VV->length_violation=-VRP_INFINITY;
00667
00668 this->normalize_route_numbers();
00669
00670 for(i=1;i<=this->total_number_of_routes;i++)
00671 {
00672 if(this->route[i].length>this->orig_max_route_length)
00673 {
00674
00675 if(this->route[i].length - this->orig_max_route_length > VV->length_violation)
00676 VV->length_violation = this->route[i].length - this->orig_max_route_length ;
00677
00678 is_feasible=false;
00679 }
00680
00681 if(this->route[i].load>this->orig_max_veh_capacity)
00682 {
00683
00684 if(this->route[i].load - this->orig_max_veh_capacity > VV->capacity_violation)
00685 VV->capacity_violation = this->route[i].load - this->orig_max_veh_capacity ;
00686
00687 is_feasible=false;
00688
00689 }
00690
00691 }
00692
00693 return is_feasible;
00694
00695 }
00696
00697
00698
00699
00700 void VRP::refresh_routes()
00701 {
00706
00707
00708 int i, n, cnt;
00709 double len=0;
00710 double rlen=0;
00711 double tot_len=0;
00712 int next_node;
00713 int current_node, current_route, route_start, current_start, current_end;
00714 int total_load = 0;
00715 int current_load = 0;
00716
00717
00718 n=num_nodes;
00719
00720 i = 1;
00721 cnt = 0;
00722 route_start = -next_array[VRPH_DEPOT];
00723
00724 current_node = route_start;
00725 current_route = route_num[current_node];
00726
00727 current_start = route[current_route].start;
00728 current_end = route[current_route].end;
00729
00730 total_load+=nodes[current_node].demand;
00731 current_load+=nodes[current_node].demand;
00732 len+=d[VRPH_DEPOT][current_node];
00733 rlen+=d[VRPH_DEPOT][current_node];
00734
00735 while(route_start != 0 && i<num_nodes+1)
00736 {
00737
00738
00739 if(next_array[current_node]==current_node)
00740 {
00741
00742 fprintf(stderr,"(2)Self loop found in next array(%d)\n",current_node);
00743 report_error("%s: Error in refresh_routes()\n",__FUNCTION__);
00744 }
00745
00746 if(next_array[current_node]==0)
00747 {
00748
00749
00750 len+=d[current_node][VRPH_DEPOT];
00751 rlen+=d[current_node][VRPH_DEPOT];
00752 current_route=route_num[current_node];
00753 route[current_route].length = rlen;
00754 route[current_route].load = current_load;
00755 total_route_length=len;
00756
00757
00758 return;
00759 }
00760
00761 if(next_array[current_node]>0)
00762 {
00763
00764
00765 next_node = next_array[current_node];
00766 len+=d[current_node][next_node];
00767 rlen+=d[current_node][next_node];
00768 current_node=next_node;
00769 total_load+=nodes[current_node].demand;
00770 current_load+=nodes[current_node].demand;
00771 cnt++;
00772 }
00773 else
00774 {
00775
00776
00777 len+=d[current_node][VRPH_DEPOT];
00778
00779 rlen+=d[current_node][VRPH_DEPOT];
00780 tot_len+=rlen;
00781 current_route=route_num[current_node];
00782 current_end = route[current_route].end;
00783
00784 route[current_route].length = rlen;
00785 route[current_route].load = current_load;
00786
00787 i++;
00788 route_start = - (next_array[current_node]);
00789 current_route = route_num[route_start];
00790 current_start = route[current_route].start;
00791 current_end = route[current_route].end;
00792 if(route_start != current_start)
00793 {
00794 fprintf(stderr,"Route %d: %d != %d\n",current_route, route_start, current_start);
00795 report_error("%s: Error in refresh_routes()\n",__FUNCTION__);
00796 }
00797
00798 current_node = route_start;
00799 total_load+=nodes[current_node].demand;
00800
00801 current_load=nodes[current_node].demand;
00802 len+=d[VRPH_DEPOT][current_node];
00803 rlen=d[VRPH_DEPOT][current_node];
00804 cnt++;
00805 }
00806 }
00807
00808
00809 return;
00810
00811
00812 }
00813 void VRP::create_pred_array()
00814 {
00818 int i,j;
00819
00820 i=VRPH_DEPOT;
00821 j=next_array[i];
00822 while(j!=VRPH_DEPOT)
00823 {
00824 if(j>0)
00825 pred_array[j]=i;
00826 else
00827 pred_array[-j]=-i;
00828
00829 i=VRPH_ABS(j);
00830 j=next_array[i];
00831
00832 }
00833
00834 pred_array[j]=-i;
00835 return;
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850 }
00851
00852
00853
00854
00855
00856 bool VRP::get_segment_info(int a, int b, struct VRPSegment *S)
00857 {
00867
00868
00869 if(a==b)
00870 {
00871
00872 S->segment_start=b;
00873 S->segment_end=a;
00874 S->len=0;
00875 S->load=nodes[a].demand;
00876 S->num_custs=1;
00877 return true;
00878
00879 }
00880
00881 S->len=0;
00882 S->segment_start=a;
00883 S->segment_end=b;
00884 S->num_custs=0;
00885
00886
00887
00888 if(a==VRPH_DEPOT)
00889 {
00890
00891
00892 S->segment_start=route[route_num[b]].start;
00893 S->segment_end=b;
00894 S->len+=d[VRPH_DEPOT][S->segment_start];
00895
00896 }
00897
00898 if(b==VRPH_DEPOT)
00899 {
00900
00901
00902 S->segment_start=a;
00903 S->segment_end=route[route_num[a]].end;
00904
00905 }
00906
00907
00908
00909 int current_node = S->segment_start;
00910 int next_node;
00911
00912 S->load=nodes[current_node].demand;
00913 if(current_node!=dummy_index)
00914 S->num_custs++;;
00915
00916 while(current_node != S->segment_end)
00917 {
00918 next_node = VRPH_MAX(next_array[current_node],0);
00919 S->len+=d[current_node][next_node];
00920 current_node=next_node;
00921 S->load+=nodes[next_node].demand;
00922 if(current_node!= dummy_index)
00923 S->num_custs += 1;
00924 }
00925
00926 if(b==VRPH_DEPOT)
00927 S->len+=d[S->segment_end][VRPH_DEPOT];
00928
00929
00930
00931 return true;
00932
00933 }
00934
00935 int VRP::get_string_end(int a, int len)
00936 {
00942
00943 int ctr=1;
00944 int current_node;
00945
00946 current_node=a;
00947 while(ctr<len)
00948 {
00949 current_node= next_array[current_node];
00950 if(current_node<0)
00951 {
00952
00953 return -1;
00954 }
00955 ctr++;
00956
00957 }
00958
00959 return current_node;
00960 }
00961
00962 void VRP::reverse_route(int i)
00963 {
00969
00970
00971 int current_node, start_node, last_node, prev_route, next_route, temp;
00972 int orig_start, orig_end ;
00973
00974 if(i<=0)
00975 {
00976
00977 fprintf(stderr,"Reversing route of negative index?? i=%d\n",i);
00978 report_error("%s: Can't reverse route!\n",__FUNCTION__);
00979 }
00980
00981 orig_end = route[i].end;
00982 orig_start = route[i].start;
00983
00984 start_node=orig_start;
00985 current_node=start_node;
00986
00987
00988
00989
00990 temp = next_array[current_node];
00991 prev_route = -pred_array[current_node];
00992
00993
00994
00995 pred_array[current_node] = temp;
00996 current_node = temp;
00997 while( (temp = next_array[current_node]) >0)
00998 {
00999 next_array[current_node] = pred_array[current_node];
01000 pred_array[current_node] = temp;
01001 current_node = temp;
01002 }
01003
01004
01005
01006 temp= pred_array[current_node];
01007 last_node = current_node;
01008 next_route = -next_array[last_node];
01009
01010 route[i].end=orig_start;
01011 route[i].start=orig_end;
01012
01013 #if REVERSE_DEBUG
01014 printf("start_node: %d; last_node: %d; prev_route: %d; next_route: %d\n",start_node,last_node, prev_route,next_route);
01015 #endif
01016
01017
01018 next_array[prev_route]=-last_node;
01019 pred_array[next_route]=-start_node;
01020 next_array[start_node]=-next_route;
01021 pred_array[last_node]= -prev_route;
01022 next_array[last_node]=temp;
01023 next_array[prev_route]=-last_node;
01024
01025
01026 if(!this->symmetric)
01027 {
01028
01029 this->refresh_routes();
01030 }
01031
01032 return;
01033 }
01034
01035
01036
01037 bool VRP::postsert_dummy(int i)
01038 {
01039
01040
01041 int pre_i, post_i, dummy;
01042 int start, end, start_i, end_i;
01043 int n, i_route;
01044
01045
01046
01047
01048 if(i<=VRPH_DEPOT || i>matrix_size)
01049 report_error("%s: input doesn't make sense\n",__FUNCTION__);
01050
01051 n= num_nodes;
01052 dummy= dummy_index;
01053
01054 i_route= route_num[i];
01055
01056 start_i= route[i_route].start;
01057 end_i= route[i_route].end;
01058
01059 start=start_i;
01060 if(end_i==i)
01061
01062 end=dummy;
01063 else
01064 end=end_i;
01065
01066
01067 pre_i= pred_array[i];
01068
01069 post_i= next_array[i];
01070
01071 next_array[i] = dummy;
01072 next_array[dummy] = post_i;
01073
01074
01075
01076
01077 pred_array[dummy]=i;
01078
01079
01080 if(post_i>=0)
01081 pred_array[post_i]=dummy;
01082 else
01083
01084 pred_array[-post_i]=-dummy;
01085
01086
01087
01088
01089
01090 route_num[dummy]=i_route;
01091 route[i_route].end=end;
01092 route[i_route].start=start;
01093
01094 return true;
01095 }
01096
01097 bool VRP::presert_dummy(int i)
01098 {
01099
01100
01101
01102 int pre_i, post_i;
01103 int start, end, start_i, end_i, dummy;
01104 int n, i_route;
01105
01106
01107 n= num_nodes ;
01108 dummy = dummy_index;
01109
01110 if(i<=VRPH_DEPOT)
01111 report_error("%s: bad index\n",__FUNCTION__);
01112
01113 i_route = route_num[i];
01114
01115 start_i= route[i_route].start;
01116 end_i= route[i_route].end;
01117
01118
01119
01120 start=start_i;
01121 end=end_i;
01122
01123 if(start==i)
01124
01125 start=dummy;
01126
01127
01128 pre_i= pred_array[i];
01129
01130 post_i= next_array[i];
01131
01132
01133
01134 next_array[dummy]=i;
01135 pred_array[i]=dummy;
01136 pred_array[dummy]=pre_i;
01137
01138
01139 if(pre_i>0)
01140 next_array[pre_i]=dummy;
01141 else
01142
01143 next_array[VRPH_ABS(pre_i)]=-dummy;
01144
01145
01146 route_num[dummy]=i_route;
01147 route[i_route].end=end;
01148 route[i_route].start=start;
01149
01150
01151
01152 return true;
01153
01154 }
01155 bool VRP::remove_dummy()
01156 {
01157 int pre_d, post_d, d_route, dummy, d_start, d_end, n;
01158
01159 n= num_nodes;
01160 dummy= dummy_index;
01161
01162 post_d= next_array[dummy];
01163 pre_d= pred_array[dummy];
01164
01165 if(post_d > dummy || post_d < -dummy || pre_d > dummy || pre_d < -dummy)
01166 {
01167 fprintf(stderr,"post_d= %d; pre_d=%d\n",post_d,pre_d);
01168 report_error("%s: invalid indices\n",__FUNCTION__);
01169 }
01170
01171 d_route= route_num[dummy];
01172 d_start= route[d_route].start;
01173 d_end= route[d_route].end;
01174
01175 if(d_start==dummy)
01176 {
01177 if(post_d<0)
01178 {
01179 report_error("%s: post_d error in removal\n",__FUNCTION__);
01180 }
01181
01182 route[d_route].start=post_d;
01183 }
01184
01185 if(d_end==dummy)
01186 {
01187 if(pre_d<0)
01188 {
01189 report_error("%s: pre_d error in removal\n",__FUNCTION__);
01190 }
01191 route[d_route].end=pre_d;
01192 }
01193
01194 next_array[VRPH_ABS(pre_d)]=post_d;
01195
01196 if(d_start==dummy)
01197 next_array[VRPH_ABS(pre_d)]=-post_d;
01198
01199
01200 pred_array[VRPH_ABS(post_d)]=pre_d;
01201
01202 if(d_end==dummy)
01203 pred_array[VRPH_ABS(post_d)]=-pre_d;
01204
01205
01206 return true;
01207 }
01208 bool VRP::create_default_routes()
01209 {
01216
01217
01218 int i,n;
01219 bool is_feasible=true;
01220
01221 violation.capacity_violation = 0;
01222 violation.length_violation = 0;
01223
01224
01225 total_route_length=0;
01226
01227 n = this->num_original_nodes;
01228
01229
01230 routed[VRPH_DEPOT]=true;
01231
01232 next_array[VRPH_DEPOT] = -1;
01233 for(i=1;i<=n;i++)
01234 {
01235 next_array[i] = -(i+1);
01236 total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]);
01237
01238 route_num[i]=i;
01239 route[i].start=i;
01240 route[i].end=i;
01241 route[i].load= nodes[i].demand;
01242 route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT];
01243
01244
01245
01246 if(route[i].load > max_veh_capacity)
01247 is_feasible=false;
01248 if(route[i].length > max_route_length)
01249 is_feasible=false;
01250
01251 route[i].num_customers=1;
01252
01253 routed[i]=true;
01254
01255 }
01256
01257
01258 next_array[n]=VRPH_DEPOT;
01259
01260 route_num[VRPH_DEPOT]=0;
01261
01262
01263 create_pred_array();
01264
01265 total_number_of_routes=n;
01266
01267
01268
01269
01270
01271 if(is_feasible == false)
01272 {
01273
01274
01275
01276
01277 for(i=1;i<=n;i++)
01278 {
01279 routed[i]=false;
01280
01281
01282
01283 if(route[i].load > max_veh_capacity)
01284 {
01285
01286
01287 printf("Default routes load violation: %d > %d\n",route[i].load,
01288 max_veh_capacity);
01289
01290 if( (route[i].load - max_veh_capacity) >
01291 violation.capacity_violation)
01292 {
01293 violation.capacity_violation =
01294 (route[i].load - max_veh_capacity);
01295 }
01296
01297 }
01298 if(route[i].length > max_route_length)
01299 {
01300
01301
01302 printf("Default routes length violation: %f > %f\n",route[i].length,
01303 max_route_length);
01304
01305 if( (route[i].length -max_route_length ) >
01306 violation.length_violation)
01307 {
01308 violation.length_violation =
01309 (route[i].length -max_route_length );
01310 }
01311
01312
01313 }
01314
01315
01316 }
01317
01318 return false;
01319 }
01320 else
01321 {
01322
01323 for(i=1;i<=n;i++)
01324 {
01325 routed[i]=true;
01326 }
01327 return true;
01328 }
01329
01330
01331 }
01332
01333
01334 bool VRP::create_default_routes(int day)
01335 {
01342
01343 int i,n;
01344 bool is_feasible=true;
01345
01346 violation.capacity_violation = 0;
01347 violation.length_violation = 0;
01348 total_route_length=0;
01349
01350 n = this->num_original_nodes;
01351 this->num_nodes=n;
01352
01353
01354
01355 routed[VRPH_DEPOT]=true;
01356 next_array[VRPH_DEPOT] = -1;
01357 for(i=1;i<=n;i++)
01358 {
01359
01360 next_array[i] = -(i+1);
01361 total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]);
01362
01363 route_num[i]=i;
01364 route[i].start=i;
01365 route[i].end=i;
01366 route[i].load= nodes[i].demand;
01367 route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT];
01368 route[i].num_customers=1;
01369 routed[i]=true;
01370
01371 }
01372
01373
01374 next_array[n]=VRPH_DEPOT;
01375
01376 route_num[VRPH_DEPOT]=0;
01377
01378
01379 create_pred_array();
01380
01381 total_number_of_routes=n;
01382
01383
01384 for(i=1;i<=n;i++)
01385 {
01386 if(this->nodes[i].daily_demands[day]==-1)
01387 this->eject_node(i);
01388 }
01389
01390
01391 this->normalize_route_numbers();
01392
01393
01394 for(i=1;i<=this->total_number_of_routes;i++)
01395 {
01396 if(route[i].load > max_veh_capacity)
01397 is_feasible=false;
01398 if(route[i].length > max_route_length)
01399 is_feasible=false;
01400 }
01401
01402 if(is_feasible == false)
01403 {
01404
01405
01406
01407
01408 for(i=1;i<=n;i++)
01409 {
01410 routed[i]=false;
01411
01412
01413
01414 if(route[i].load > max_veh_capacity)
01415 {
01416
01417
01418 printf("Default routes load violation: %d > %d\n",route[i].load,
01419 max_veh_capacity);
01420
01421 if( (route[i].load - max_veh_capacity) >
01422 violation.capacity_violation)
01423 {
01424 violation.capacity_violation =
01425 (route[i].load - max_veh_capacity);
01426 }
01427
01428 }
01429 if(route[i].length > max_route_length)
01430 {
01431
01432
01433 printf("Default routes length violation: %f > %f\n",route[i].length,
01434 max_route_length);
01435
01436 if( (route[i].length -max_route_length ) >
01437 violation.length_violation)
01438 {
01439 violation.length_violation =
01440 (route[i].length -max_route_length );
01441 }
01442 }
01443 }
01444 return false;
01445 }
01446 else
01447 return true;
01448
01449 }
01450
01451
01452 int VRP::count_num_routes()
01453 {
01457
01458 int current,next;
01459 int num=0;
01460
01461 current=VRPH_DEPOT;
01462 next=-1;
01463
01464 while(next!=VRPH_DEPOT)
01465 {
01466 next= next_array[current];
01467 if(next<0)
01468 {
01469
01470 num++;
01471 current=-next;
01472 }
01473 else
01474 current=next;
01475 }
01476
01477 return num;
01478
01479
01480 }
01481
01482 bool VRP::perturb()
01483 {
01487
01488 int i, n, pre, post;
01489 int current,next;
01490 n= num_nodes;
01491 VRPNeighborElement *v;
01492 VRPMove M;
01493
01494 Postsert Postsert;
01495 Presert Presert;
01496
01497 v=new VRPNeighborElement[n];
01498
01499 current=VRPH_ABS(next_array[VRPH_DEPOT]);
01500 i=0;
01501 while(current!=VRPH_DEPOT)
01502 {
01503 pre= VRPH_MAX(VRPH_DEPOT,pred_array[current]);
01504 post= VRPH_MAX(VRPH_DEPOT, next_array[current]);
01505
01506
01507 v[i].val=((double) nodes[current].demand) /
01508 (VRPH_EPSILON+d[pre][current]+d[current][post]-d[pre][post]);
01509 v[i].position=current;
01510 i++;
01511 next=VRPH_ABS(next_array[current]);
01512 current=next;
01513 }
01514
01515
01516
01517
01518 qsort(v,n,sizeof(VRPNeighborElement), VRPNeighborCompare);
01519
01520 int m=(int) (VRPH_MIN(30,n/5));
01521
01522
01523
01524
01525 int a,b,c,j,k,node1=0, node2=0,b_route,b_load,k_demand, jj, ll;
01526 double best_savings,b_len,jk,kl,jl;
01527
01528
01529 for(j=0;j<m;j++)
01530 {
01531 best_savings=VRP_INFINITY;
01532
01533 k=v[j].position;
01534 k_demand= nodes[k].demand;
01535
01536
01537
01538
01539 jj=VRPH_MAX(pred_array[k],0);
01540 ll=VRPH_MAX(next_array[k],0);
01541
01542 jk= d[jj][k];
01543 kl= d[k][ll];
01544 jl= d[jj][ll];
01545
01546
01547
01548
01549 b=VRPH_ABS(next_array[VRPH_DEPOT]);
01550 while(b!=VRPH_DEPOT)
01551 {
01552 b_route= route_num[b];
01553 b_load= route[b].load;
01554 b_len= route[b].length;
01555
01556 a=VRPH_MAX(pred_array[b],0);
01557 c=VRPH_MAX(next_array[b],0);
01558
01559 if(a!=k && b!=k && c!=k)
01560 {
01561
01562 if(a!=VRPH_DEPOT)
01563 {
01564 if(Postsert.evaluate(this,k,a,&M)==true && M.savings<best_savings)
01565 {
01566 best_savings=M.savings;
01567 node1=a;node2=b;
01568 }
01569 }
01570 else
01571 {
01572
01573 if(Presert.evaluate(this,k,b,&M)==true && M.savings<best_savings)
01574 {
01575 best_savings=M.savings;
01576 node1=a;node2=b;
01577 }
01578
01579
01580
01581 }
01582
01583
01584 if(b!=VRPH_DEPOT)
01585 {
01586 if(Postsert.evaluate(this,k,b,&M)==true && M.savings<best_savings)
01587 {
01588 best_savings=M.savings;
01589 node1=b;node2=c;
01590 }
01591 }
01592 else
01593 {
01594
01595 if(Presert.evaluate(this,k,c,&M)==true && M.savings<best_savings)
01596 {
01597 best_savings=M.savings;
01598 node1=b;node2=c;
01599 }
01600
01601
01602 }
01603
01604 }
01605
01606 b=VRPH_ABS(next_array[b]);
01607
01608 }
01609 if(best_savings!=VRP_INFINITY)
01610 {
01611
01612 if(node1!=VRPH_DEPOT)
01613 Postsert.move(this,k,node1);
01614 else
01615 Presert.move(this,k,node2);
01616 }
01617
01618 }
01619
01620
01621
01622 delete [] v;
01623 return true;
01624 }
01625
01626
01627 bool VRP::eject_node(int j)
01628 {
01633
01634 if(j<=VRPH_DEPOT || routed[j]==false )
01635 {
01636 fprintf(stderr,"Tried to eject index %d\n",j);
01637 report_error("%s: VRPH_DEPOT or unrouted index used in eject node\n",__FUNCTION__);
01638 }
01639
01640
01641 int k=j;
01642
01643 int c, e, k_route, c_route, e_route, k_start, k_end,flag;
01644 double change,ck,ke,ce;
01645
01646 flag=0;
01647
01648
01649 routed[k]=false;
01650
01651
01652
01653 c= pred_array[k];
01654 e= next_array[k];
01655
01656 k_route= route_num[k];
01657 c_route= route_num[VRPH_ABS(c)];
01658 e_route= route_num[VRPH_ABS(e)];
01659
01660 k_start= route[k_route].start;
01661 k_end= route[k_route].end;
01662
01663 if(k_start==k && k_end!=k)
01664 {
01665
01666 route[k_route].start=VRPH_ABS(e);
01667 next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01668 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01669
01670 }
01671
01672 if(k_end==k && k_start!=k)
01673 {
01674
01675 route[k_route].end=c;
01676 next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01677 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01678 }
01679
01680 if(k_end==k && k_start==k)
01681 {
01682
01683 next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01684 pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01685 flag=1;
01686 }
01687
01688 if(k_start!=k && k_end!=k)
01689 {
01690
01691 next_array[c]=e;
01692 pred_array[e]=c;
01693
01694 }
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704 if(e<0)
01705 e=VRPH_DEPOT;
01706 if(c<0)
01707 c=VRPH_DEPOT;
01708
01709 ce= this->d[c][e];
01710 ck= this->d[c][k];
01711 ke= this->d[k][e];
01712
01713 change=ce-(ck+ke);
01714 total_route_length+=change;
01715 route[k_route].length+=change;
01716 route[k_route].load-= nodes[k].demand;
01717
01718
01719 route[k_route].num_customers--;
01720 num_nodes--;
01721
01722 if(flag)
01723
01724 this->total_number_of_routes--;
01725
01726 route_num[k]=-1;
01727 normalize_route_numbers();
01728
01729 return true;
01730 }
01731
01732 bool VRP::eject_route(int r, int *route_buff)
01733 {
01738
01739 int current,cnt;
01740
01741
01742 current=this->route[r].start;
01743 cnt=0;
01744 while(current>0)
01745 {
01746 route_buff[cnt]=current;
01747 cnt++;
01748 current=this->next_array[current];
01749 }
01750 route_buff[cnt]=-1;
01751
01752
01753
01754
01755
01756 for(int i=0;i<cnt;i++)
01757 this->eject_node(route_buff[i]);
01758
01759
01760 return true;
01761
01762 }
01763 bool VRP::check_move(VRPMove *M, int rules)
01764 {
01765
01770
01771
01772 double savings;
01773
01774 savings=M->savings;
01775
01776
01777 if(this->forbid_tiny_moves)
01778 {
01779
01780 if(savings>-VRPH_EPSILON && savings < VRPH_EPSILON)
01781 return false;
01782 }
01783
01784
01785
01786 if( (rules & VRPH_FREE) == VRPH_FREE )
01787 {
01788
01789 return true;
01790 }
01791
01792
01793
01794
01795 if( (rules & VRPH_DOWNHILL) == VRPH_DOWNHILL)
01796 {
01797 if(savings<-VRPH_EPSILON )
01798 return true;
01799 else
01800 return false;
01801
01802 }
01803
01804
01805
01806
01807 if( rules & VRPH_RECORD_TO_RECORD )
01808 {
01809 if(savings<=-VRPH_EPSILON)
01810 return true;
01811
01812
01813
01814
01815
01816 if(!has_service_times)
01817 {
01818 if(total_route_length+savings<= (1+deviation)*record)
01819 return true;
01820 else
01821 return false;
01822
01823 }
01824 else
01825 {
01826
01827
01828
01829 if( (total_route_length - total_service_time) + savings <=
01830 ((1+deviation)*(record-total_service_time)))
01831 return true;
01832 else
01833 return false;
01834 }
01835 }
01836
01837 if( rules & VRPH_SIMULATED_ANNEALING )
01838 {
01839 if(M->evaluated_savings==true)
01840 return true;
01841
01842
01843 if( exp(- (M->savings / this->temperature)) > lcgrand(10) )
01844 return true;
01845 else
01846 return false;
01847
01848 }
01849
01850
01851
01852 report_error("%s: didn't return yet!\n",__FUNCTION__);
01853
01854 return false;
01855
01856
01857 }
01858
01859
01860 bool VRP::is_feasible(VRPMove *M, int rules)
01861 {
01867
01868 for(int i=0;i<M->num_affected_routes;i++)
01869 {
01870 if( (M->route_lens[i]>this->max_route_length) || (M->route_loads[i]>this->max_veh_capacity) )
01871 return false;
01872 }
01873
01874 return true;
01875
01876 }
01877 bool VRP::inject_node(int j)
01878 {
01883
01884 if(j==VRPH_DEPOT)
01885 report_error("%s: Can't inject VRPH_DEPOT!!\n",__FUNCTION__);
01886
01887 int edge[4];
01888 double costs[4];
01889 this->find_cheapest_insertion(j,edge,costs,VRPH_USE_NEIGHBOR_LIST);
01890
01891 this->insert_node(j,edge[0],edge[1]);
01892
01893 return true;
01894
01895 }
01896
01897 bool VRP::insert_node(int j, int i, int k)
01898 {
01903
01904 if(routed[j] || !routed[i] || !routed[k])
01905 {
01906 fprintf(stderr,"insert_node(%d,%d,%d)\n",j,i,k);
01907 report_error("%s: Improper nodes used in insert_node\n",__FUNCTION__);
01908 }
01909
01910 double increase;
01911 int r;
01912 routed[j]=true;
01913
01914 if(i==k && k==VRPH_DEPOT)
01915 {
01916
01917 int last_node=VRPH_ABS(pred_array[VRPH_DEPOT]);
01918 this->next_array[last_node]=-j;
01919 this->pred_array[j]=-last_node;
01920 this->next_array[j]=VRPH_DEPOT;
01921 this->pred_array[VRPH_DEPOT]=-j;
01922
01923 increase=this->d[0][j]+this->d[j][0];
01924 this->total_number_of_routes++;
01925 this->route_num[j]=this->total_number_of_routes;
01926 this->route[total_number_of_routes].length=increase;
01927 this->route[total_number_of_routes].load=nodes[j].demand;
01928 this->route[total_number_of_routes].num_customers=1;
01929 this->route[total_number_of_routes].start=j;
01930 this->route[total_number_of_routes].end=j;
01931 this->num_nodes++;
01932 this->total_route_length+=increase;
01933
01934 return true;
01935
01936
01937
01938 }
01939
01940 if(i!=VRPH_DEPOT && k!=VRPH_DEPOT)
01941 {
01942
01943
01944 if(VRPH_MAX(VRPH_DEPOT,next_array[i]) != k || VRPH_MAX(VRPH_DEPOT,pred_array[k])!=i)
01945 {
01946 fprintf(stderr,"edge doesn't exist: next(i(%d)) != k(%d)\n",i,k);
01947 report_error("%s\n",__FUNCTION__);
01948 }
01949
01950 increase=d[i][j]+d[j][k]-d[i][k];
01951 r=route_num[i];
01952
01953
01954 num_nodes++;
01955 next_array[i]=j;
01956 pred_array[j]=i;
01957 next_array[j]=k;
01958 pred_array[k]=j;
01959 route_num[j]=r;
01960 route[r].length+=increase;
01961 route[r].load+= nodes[j].demand;
01962 route[r].num_customers++;
01963 total_route_length+=increase;
01964
01965 return true;
01966 }
01967
01968 if(i==VRPH_DEPOT)
01969 {
01970
01971 increase=d[i][j]+d[j][k]-d[i][k];
01972
01973 r=route_num[k];
01974
01975
01976 int pre=VRPH_ABS(pred_array[k]);
01977
01978
01979 num_nodes++;
01980
01981 next_array[pre]=-j;
01982 pred_array[j]=-pre;
01983 next_array[j]=k;
01984 pred_array[k]=j;
01985 route_num[j]=r;
01986 route[r].start=j;
01987 route[r].length+=increase;
01988 route[r].load+= nodes[j].demand;
01989 route[r].num_customers++;
01990
01991 total_route_length+=increase;
01992 return true;
01993
01994 }
01995
01996 if(k==VRPH_DEPOT)
01997 {
01998
01999 increase=d[i][j]+d[j][k]-d[i][k];
02000 r=route_num[i];
02001
02002 int post=VRPH_ABS(next_array[i]);
02003
02004
02005 num_nodes++;
02006
02007 next_array[i]=j;
02008 pred_array[j]=i;
02009 next_array[j]=-post;
02010 pred_array[post]=-j;
02011 route_num[j]=r;
02012 route[r].length+=increase;
02013 route[r].load+= nodes[j].demand;
02014 route[r].end=j;
02015 route[r].num_customers++;
02016 total_route_length+=increase;
02017 return true;
02018
02019 }
02020
02021 return false;
02022
02023
02024 }
02025
02026 void VRP::find_cheapest_insertion(int j, int *edge, double *costs, int rules)
02027 {
02041
02042 int h,i,k,m, best_route, new_route, next_node;
02043 double ij,ik,jk,min_feasible_increase, increase, min_increase;
02044
02045 this->normalize_route_numbers();
02046 best_route = -1;
02047 k=-1;
02048
02049
02050
02051 min_feasible_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT];
02052 min_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT];
02053 edge[0]=edge[1]=edge[2]=edge[3]=VRPH_DEPOT;
02054
02055
02056 if(!(rules & VRPH_USE_NEIGHBOR_LIST))
02057 {
02058 i=VRPH_DEPOT;
02059 next_node=VRPH_ABS(next_array[i]);
02060 int cnt=0;
02061 while(next_node != VRPH_DEPOT)
02062 {
02063 cnt++;
02064 if(next_node > 0)
02065 {
02066 k=next_node;
02067
02068
02069
02070 ij= d[i][j];
02071 ik= d[i][k];
02072 jk= d[j][k];
02073 increase=ij+jk-ik;
02074
02075 if(increase<min_increase)
02076 {
02077 min_increase=increase;
02078 edge[2]=i;
02079 edge[3]=k;
02080 }
02081 if(increase<min_feasible_increase)
02082 {
02083
02084 new_route= route_num[k];
02085
02086 if( (route[new_route].length+increase <= max_route_length) &&
02087 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02088 {
02089 edge[0]=i;
02090 edge[1]=k;
02091 best_route=new_route;
02092 min_feasible_increase=increase;
02093
02094 }
02095
02096 }
02097
02098 i=k;
02099
02100 }
02101 else
02102 {
02103
02104
02105 k=VRPH_DEPOT;
02106
02107 ij= d[i][j];
02108 ik= d[i][k];
02109 jk= d[j][k];
02110 increase=ij+jk-ik;
02111
02112 if(increase<min_increase)
02113 {
02114 min_increase=increase;
02115 edge[2]=i;
02116 edge[3]=k;
02117 }
02118
02119 if(increase<min_feasible_increase)
02120 {
02121
02122 new_route= route_num[i];
02123
02124 if( (route[new_route].length+increase <= max_route_length) &&
02125 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02126 {
02127 edge[0]=i;
02128 edge[1]=k;
02129 best_route=new_route;
02130 min_feasible_increase=increase;
02131 }
02132 }
02133
02134
02135 i=VRPH_DEPOT;
02136 k=VRPH_ABS(next_node);
02137
02138 ij= d[i][j];
02139 ik= d[i][k];
02140 jk= d[j][k];
02141 increase=ij+jk-ik;
02142
02143 if(increase<min_increase)
02144 {
02145 min_increase=increase;
02146 edge[2]=i;
02147 edge[3]=k;
02148 }
02149
02150 if(increase<min_feasible_increase)
02151 {
02152 new_route= route_num[k];
02153
02154 if( (route[new_route].length+increase <= max_route_length )&&
02155 (route[new_route].load + nodes[j].demand <= max_veh_capacity ) )
02156 {
02157 edge[0]=i;
02158 edge[1]=k;
02159 best_route=new_route;
02160 min_feasible_increase=increase;
02161 }
02162
02163 }
02164 i=k;
02165 }
02166
02167
02168 next_node= next_array[i];
02169
02170 }
02171
02172 k=VRPH_DEPOT;
02173 ij= d[i][j];
02174 ik= d[i][k];
02175 jk= d[j][k];
02176 increase=ij+jk-ik;
02177
02178 if(increase<min_increase)
02179 {
02180 min_increase=increase;
02181 edge[2]=i;
02182 edge[3]=k;
02183 }
02184
02185 if(increase<min_feasible_increase)
02186 {
02187
02188 new_route= route_num[i];
02189
02190 if( (route[new_route].length+increase <= max_route_length) &&
02191 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02192 {
02193 edge[0]=i;
02194 edge[1]=k;
02195 best_route=new_route;
02196 min_feasible_increase=increase;
02197 }
02198 }
02199
02200 costs[0]=min_feasible_increase;
02201 costs[1]=min_increase;
02202 return;
02203
02204 }
02205 else
02206 {
02207
02208
02209 for(m=0;m<neighbor_list_size;m++)
02210 {
02211 i=nodes[j].neighbor_list[m].position;
02212 if(routed[i])
02213 {
02214 h=VRPH_MAX(VRPH_DEPOT,pred_array[i]);
02215 increase=d[h][j]+d[j][i]-d[h][i];
02216
02217 if(increase<min_increase)
02218 {
02219 min_increase=increase;
02220 edge[2]=h;
02221 edge[3]=i;
02222 }
02223
02224 if(increase<min_feasible_increase)
02225 {
02226
02227 if(i!=VRPH_DEPOT)
02228 new_route= route_num[i];
02229 else
02230 new_route= route_num[k];
02231
02232
02233 if( (route[new_route].length+increase <= max_route_length) &&
02234 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02235 {
02236 edge[0]=h;
02237 edge[1]=i;
02238 best_route=new_route;
02239 min_feasible_increase=increase;
02240 }
02241 }
02242
02243 k=VRPH_MAX(VRPH_DEPOT,next_array[i]);
02244 increase=d[i][j]+d[j][k]-d[i][k];
02245
02246 if(increase<min_increase)
02247 {
02248 min_increase=increase;
02249 edge[2]=i;
02250 edge[3]=k;
02251 }
02252
02253 if(increase<min_feasible_increase)
02254 {
02255
02256 if(i!=VRPH_DEPOT)
02257 new_route= route_num[i];
02258 else
02259 new_route= route_num[k];
02260
02261
02262 if( (route[new_route].length+increase <= max_route_length) &&
02263 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02264 {
02265 edge[0]=i;
02266 edge[1]=k;
02267 best_route=new_route;
02268 min_feasible_increase=increase;
02269 }
02270 }
02271 }
02272 }
02273
02274 costs[0]=min_feasible_increase;
02275 costs[1]=min_increase;
02276 return;
02277 }
02278 }
02279 int VRP::inject_set(int num, int *nodelist, int rules, int attempts)
02280 {
02287
02288 if(rules != VRPH_RANDOM_SEARCH && rules != VRPH_REGRET_SEARCH)
02289 report_error("%s: invalid rules\n",__FUNCTION__);
02290
02291 int i,j,k;
02292 double best_obj=VRP_INFINITY;
02293 int *best_sol, *orderings,*best_ordering, *start_sol;
02294 int best_index=0;
02295
02296 best_sol=orderings=best_ordering=start_sol=NULL;
02297
02298 for(i=0;i<num;i++)
02299 {
02300 if(nodelist[i]==VRPH_DEPOT)
02301 {
02302 fprintf(stderr,"nodelist[%d] of %d=VRPH_DEPOT\n",i,num);
02303 report_error("%s: Cannot inject VRPH_DEPOT!\n",__FUNCTION__);
02304 }
02305 }
02306
02307 best_sol=new int[3+(this->num_nodes)+num];
02308 start_sol=new int[3+(this->num_nodes)+num];
02309 this->export_solution_buff(start_sol);
02310 this->import_solution_buff(start_sol);
02311
02312 if(rules==VRPH_RANDOM_SEARCH)
02313 {
02314
02315 orderings=new int[num];
02316 best_ordering=new int[num];
02317 int edge[4];
02318 double costs[2];
02319 for(i=0;i<num;i++)
02320 orderings[i]=i;
02321
02322 for(i=0;i<attempts;i++)
02323 {
02324
02325 random_permutation(orderings,num);
02326 for(j=0;j<num;j++)
02327 {
02328 if(nodelist[orderings[j]]==VRPH_DEPOT)
02329 {
02330 fprintf(stderr,"inject_set: Bizarre nodelist[%d]:\n",j);
02331 for(k=0;k<num;k++)
02332 fprintf(stderr,"(%d,%d) ",orderings[k],nodelist[orderings[k]]);
02333
02334 report_error("%s: Found VRPH_DEPOT in nodelist!\n",__FUNCTION__);
02335 }
02336 this->inject_node(nodelist[orderings[j]]);
02337 }
02338
02339
02340
02341 if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON)
02342 {
02343 best_index=i;
02344 best_obj=this->total_route_length;
02345 this->export_solution_buff(best_sol);
02346 memcpy(best_ordering,orderings,num*sizeof(int));
02347 }
02348
02349
02350 for(j=0;j<num;j++)
02351 this->eject_node(nodelist[orderings[j]]);
02352 }
02353
02354 for(j=0;j<num;j++)
02355 {
02356 this->find_cheapest_insertion(nodelist[best_ordering[j]],edge,costs,VRPH_USE_NEIGHBOR_LIST);
02357 this->insert_node(nodelist[best_ordering[j]],edge[0],edge[1]);
02358 }
02359 }
02360
02361 if(rules==VRPH_REGRET_SEARCH)
02362 {
02363
02364 orderings=new int[num];
02365 best_ordering=new int[num];
02366
02367 int *current_list;
02368 current_list=new int[num];
02369 int edge[4];
02370 double costs[2];
02371 for(i=0;i<num;i++)
02372 orderings[i]=i;
02373
02374 for(i=0;i<attempts;i++)
02375 {
02376
02377 random_permutation(orderings,num);
02378 for(j=0;j<num;j++)
02379 current_list[j]=nodelist[orderings[j]];
02380
02381
02382
02383 int cycle_ctr=0;
02384 for(j=0;j<num;j++)
02385 {
02386 cycle_ctr++;
02387 if(cycle_ctr==VRPH_MAX_CYCLES)
02388 {
02389 fprintf(stderr,"Cycle encountered in REGRET SEARCH!\nReverting to original solution\n");
02390 if(best_obj<VRP_INFINITY)
02391 {
02392 this->import_solution_buff(best_sol);
02393 delete [] best_sol;
02394 delete [] orderings;
02395 delete [] best_ordering;
02396 delete [] start_sol;
02397 delete [] current_list;
02398 return 0;
02399 }
02400 else
02401 report_error("%s: Couldn't escape cycle in REGRET search!\n",__FUNCTION__);
02402 }
02403
02404
02405 this->find_cheapest_insertion(current_list[j],edge,costs,VRPH_USE_NEIGHBOR_LIST);
02406
02407 double max_ejection_cost=-VRP_INFINITY;
02408 double current_cost;
02409 int node_to_eject=-1;
02410 for(k=0;k<j;k++)
02411 {
02412
02413 if(current_list[k]!=edge[0] && current_list[k]!=edge[1])
02414 {
02415 if(( current_cost = this->ejection_cost(current_list[k])) > max_ejection_cost)
02416 {
02417 max_ejection_cost=current_cost;
02418 node_to_eject=current_list[k];
02419 if(!routed[node_to_eject])
02420 {
02421 fprintf(stderr,"%d NOT ROUTED!!\n",node_to_eject);
02422 report_error("%s: Error in inject_set\n",__FUNCTION__);
02423 }
02424 }
02425 }
02426 }
02427
02428 if(node_to_eject==-1 || costs[0]>= max_ejection_cost)
02429 {
02430
02431 this->insert_node(current_list[j],edge[0],edge[1]);
02432 }
02433 else
02434 {
02435
02436 this->eject_node(node_to_eject);
02437 this->insert_node(current_list[j],edge[0],edge[1]);
02438 current_list[j]=node_to_eject;
02439 j--;
02440 }
02441 }
02442
02443
02444 if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON)
02445 {
02446 best_index=i;
02447 best_obj=this->total_route_length;
02448 this->export_solution_buff(best_sol);
02449 memcpy(best_ordering,orderings,num*sizeof(int));
02450 }
02451
02452
02453 this->import_solution_buff(start_sol);
02454 }
02455 delete [] current_list;
02456 }
02457
02458 this->import_solution_buff(best_sol);
02459
02460 delete [] best_sol;
02461 delete [] orderings;
02462 delete [] best_ordering;
02463 delete [] start_sol;
02464
02465
02466 return best_index;
02467 }
02468
02469 void VRP::eject_neighborhood(int j, int num, int *nodelist)
02470 {
02477
02478 int i,k,cnt;
02479 int *ejected;
02480
02481 ejected=new int[this->num_nodes+1];
02482 memset(ejected,0,(this->num_nodes+1)*sizeof(int));
02483
02484 cnt=0;
02485
02486
02487 nodelist[cnt]=j;
02488 ejected[j]=1;
02489 cnt++;
02490
02491 i=0;
02492 while(cnt<num)
02493 {
02494
02495 i=(int)(lcgrand(17)*2*num);
02496
02497
02498 i=VRPH_MIN(MAX_NEIGHBORLIST_SIZE-1,i);
02499
02500
02501 k=this->nodes[j].neighbor_list[i].position;
02502 if(ejected[k]==0)
02503 {
02504 if(k!=VRPH_DEPOT)
02505 {
02506 nodelist[cnt]=k;
02507 ejected[k]=1;
02508 cnt++;
02509 }
02510 }
02511
02512 }
02513
02514
02515 for(i=0;i<cnt;i++)
02516 {
02517 if(nodelist[i]==VRPH_DEPOT)
02518 {
02519 fprintf(stderr,"Trying to eject VRPH_DEPOT??? j=%d\n",j);
02520 report_error("%s: Error in eject_neighborhood\n",__FUNCTION__);
02521 }
02522
02523 this->eject_node(nodelist[i]);
02524 }
02525
02526 #if 0
02527 this->verify_routes("After ejecting neighborhood\n");
02528 #endif
02529 delete [] ejected;
02530
02531 return;
02532 }
02533
02534
02535
02536
02537 void VRP::normalize_route_numbers()
02538 {
02544
02545 int current, end, n, R, i;
02546 int current_route;
02547 int *indices;
02548 R= count_num_routes();
02549 n= this->num_original_nodes;
02550
02551
02552
02553
02554
02555 indices=new int[n+1];
02556
02557 for(i=0;i<=n;i++)
02558 indices[i]=1;
02559
02560 int ctr=0;
02561 i=VRPH_ABS(next_array[VRPH_DEPOT]);
02562 while(i!=VRPH_DEPOT)
02563 {
02564 if(route_num[i] <= R && routed[i]==true)
02565 {
02566
02567
02568 indices[route_num[i]]=0;
02569 }
02570 else
02571 ctr++;
02572
02573 i=VRPH_ABS(next_array[i]);
02574 }
02575
02576 if(ctr==0)
02577 {
02578 delete [] indices;
02579 return;
02580 }
02581
02582
02583
02584
02585 int next_index=1;
02586 while(indices[next_index]==0)
02587 next_index++;
02588
02589
02590
02591
02592 current=VRPH_ABS(next_array[VRPH_DEPOT]);
02593
02594 while(current!=VRPH_DEPOT)
02595 {
02596 current_route= route_num[current];
02597 end= route[current_route].end;
02598
02599
02600 while(indices[next_index]==0)
02601 next_index++;
02602
02603
02604 if(current_route>R)
02605 {
02606
02607 i=VRPH_ABS(next_array[VRPH_DEPOT]);
02608 while(i!=VRPH_DEPOT)
02609 {
02610 if(route_num[i]==current_route)
02611 route_num[i]=next_index;
02612 i=VRPH_ABS(next_array[i]);
02613
02614 }
02615
02616 route[next_index].start = route[current_route].start;
02617 route[next_index].end = route[current_route].end;
02618 route[next_index].length = route[current_route].length;
02619 route[next_index].load = route[current_route].load;
02620 route[next_index].num_customers = route[current_route].num_customers;
02621
02622 next_index++;
02623
02624
02625 }
02626
02627
02628 current=VRPH_ABS(next_array[end]);
02629
02630
02631 }
02632
02633 delete [] indices;
02634
02635 return;
02636
02637
02638
02639 }
02640
02641
02642 bool VRP::create_search_neighborhood(int j, int rules)
02643 {
02648
02649 int i,k, cnt;
02650
02651
02652
02653 if( rules & VRPH_USE_NEIGHBOR_LIST )
02654 {
02655
02656
02657 search_size=0;
02658 cnt=0;
02659 for(i=0;i<neighbor_list_size;i++)
02660 {
02661
02662 k=nodes[j].neighbor_list[i].position;
02663 if( routed[k] == true)
02664 {
02665
02666 if(rules & VRPH_INTER_ROUTE_ONLY)
02667 {
02668
02669 if(route_num[k]!=route_num[j])
02670 {
02671 search_space[cnt] = k;
02672 cnt++;
02673
02674 }
02675 }
02676 else
02677 {
02678 if(rules & VRPH_INTRA_ROUTE_ONLY)
02679 {
02680
02681
02682 if(route_num[k]==route_num[j])
02683 {
02684
02685 search_space[cnt] = k;
02686 cnt++;
02687
02688 }
02689 }
02690 else
02691 {
02692
02693 search_space[cnt] = k;
02694 cnt++;
02695
02696
02697 }
02698 }
02699 }
02700 }
02701 search_size=cnt;
02702 goto randomize;
02703 }
02704
02705
02706
02707 if( (rules & VRPH_INTRA_ROUTE_ONLY) )
02708 {
02709
02710 search_space[0]=VRPH_DEPOT;
02711 search_space[1]=this->route[route_num[j]].start;
02712 for(i=2;i<=this->route[route_num[j]].num_customers;i++)
02713 search_space[i]=this->next_array[search_space[i-1]];
02714
02715 search_size=this->route[route_num[j]].num_customers+1;
02716
02717 goto randomize;
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733
02734
02735
02736
02737 }
02738
02739 if( (rules & VRPH_INTER_ROUTE_ONLY) )
02740 {
02741
02742 i=0;
02743 search_space[i]=VRPH_DEPOT;
02744 search_size=1;
02745 k=VRPH_ABS(next_array[search_space[i]]);
02746 i++;
02747
02748 for(;;)
02749 {
02750 if(k==VRPH_DEPOT)
02751
02752 goto randomize;
02753
02754 if(route_num[k]!=route_num[j])
02755 {
02756 search_space[i]=k;
02757 i++;
02758 search_size++;
02759 }
02760 k=VRPH_ABS(next_array[k]);
02761 }
02762
02763 }
02764
02765
02766
02767
02768
02769 if((rules & VRPH_FORWARD))
02770 {
02771 search_size=0;
02772
02773 i=0;
02774 search_space[i]=VRPH_ABS(next_array[route[route_num[j]].end]);
02775 search_size++;
02776
02777
02778 for(;;)
02779 {
02780 search_space[i+1]=VRPH_ABS(next_array[search_space[i]]);
02781 search_size++;
02782 if(search_space[i+1]==VRPH_DEPOT)
02783 goto randomize;
02784 i++;
02785
02786 }
02787
02788 }
02789
02790 if((rules & VRPH_BACKWARD))
02791 {
02792 search_size=0;
02793
02794 i=0;
02795 search_space[i]=VRPH_ABS(pred_array[route[route_num[j]].start]);
02796 search_size++;
02797
02798
02799 for(;;)
02800 {
02801 search_space[i+1]=VRPH_ABS(pred_array[search_space[i]]);
02802 search_size++;
02803 if(search_space[i+1]==VRPH_DEPOT)
02804 goto randomize;
02805 i++;
02806
02807 }
02808
02809 }
02810
02811
02812 search_size=0;
02813 i=0;
02814
02815 search_space[i]=VRPH_ABS(next_array[VRPH_DEPOT]);
02816 search_size++;
02817
02818 for(;;)
02819 {
02820 search_space[i+1]=VRPH_ABS(next_array[search_space[i]]);
02821 search_size++;
02822 if(search_space[i+1]==VRPH_DEPOT)
02823 goto randomize;
02824 i++;
02825 }
02826
02827 randomize:
02828
02829
02830 if(rules & VRPH_RANDOMIZED)
02831 random_permutation(search_space, search_size);
02832
02833 return true;
02834
02835
02836
02837 }
02838
02839
02840 double VRP::insertion_cost(int u, int a, int b)
02841 {
02846
02847
02848
02849
02850
02851
02852 if(a==b && b==VRPH_DEPOT)
02853 return d[VRPH_DEPOT][u]+d[u][VRPH_DEPOT];
02854
02855 if(a==u || u==b || u==VRPH_DEPOT)
02856 report_error("%s: overlap or VRPH_DEPOT found\n",__FUNCTION__);
02857
02858 if(routed[a]==false || routed[b]==false)
02859 {
02860 fprintf(stderr,"%d,%d: not routed!\n",a,b);
02861 report_error("%s: Unrouted nodes in insertion_cost\n",__FUNCTION__);
02862 }
02863
02864 if( (a!=VRPH_DEPOT) && VRPH_MAX(next_array[a],VRPH_DEPOT)!=b)
02865 {
02866 fprintf(stderr,"%d,%d: not an edge!\n",a,b);
02867 report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__);
02868 }
02869
02870 if( (a==VRPH_DEPOT) && VRPH_MAX(pred_array[b],VRPH_DEPOT)!=a)
02871 {
02872 fprintf(stderr,"%d,%d: not an edge!\n",a,b);
02873 report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__);
02874 }
02875
02876
02877 int new_route;
02878 if(a==VRPH_DEPOT)
02879 new_route=route_num[b];
02880 else
02881 new_route=route_num[a];
02882
02883 if(nodes[u].demand+route[new_route].load>max_veh_capacity)
02884 return VRP_INFEASIBLE;
02885
02886 double increase=d[a][u] + d[u][b] - d[a][b];
02887 if(route[new_route].length+increase>max_route_length)
02888 return VRP_INFEASIBLE;
02889
02890 return increase;
02891
02892 }
02893
02894
02895
02896 double VRP::ejection_cost(int u)
02897 {
02902
02903 if( u==VRPH_DEPOT)
02904 report_error("%s: Cannot eject the VRPH_DEPOT!\n",__FUNCTION__);
02905
02906 if(!routed[u])
02907 return -VRP_INFINITY;
02908
02909
02910
02911 return this->d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][u]+d[u][VRPH_MAX(next_array[u],VRPH_DEPOT)] -
02912 d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][VRPH_MAX(next_array[u],VRPH_DEPOT)];
02913
02914
02915 }
02916 void VRP::clean_route(int r, int heuristics)
02917 {
02922
02923 int i,j, r_start, r_end;
02924 double start_val, end_val, start_rlen, end_rlen;
02925
02926 OnePointMove OPM;
02927 TwoPointMove TPM;
02928 TwoOpt TO;
02929 ThreeOpt ThreeO;
02930 OrOpt OrO;
02931
02932 int rules= VRPH_INTRA_ROUTE_ONLY+VRPH_DOWNHILL+VRPH_FIRST_ACCEPT+VRPH_SAVINGS_ONLY;
02933
02934
02935 start_improving:
02936
02937 start_val = route[r].length;
02938
02939 #if CLEAN_DEBUG
02940 printf("CLEAN::start_val=%f\n",start_val);
02941 #endif
02942
02943 end_val = - VRP_INFINITY;
02944 r_start=route[r].start;
02945 r_end=route[r].end;
02946
02947 if((heuristics & ONE_POINT_MOVE)==ONE_POINT_MOVE)
02948 {
02949
02950 r_start=route[r].start;
02951
02952 opm_search:
02953 start_rlen=route[r].length;
02954 i=r_start;
02955 j=VRPH_MAX(next_array[i],0);
02956 while(i!=VRPH_DEPOT)
02957 {
02958 while(OPM.search(this,i,rules));
02959 i=j;
02960 j=VRPH_MAX(next_array[j],0);
02961 }
02962
02963 end_rlen = route[r].length;
02964 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
02965 goto opm_search;
02966
02967 #if CLEAN_DEBUG
02968 printf("CLEAN::OPM end_rlen=%f\n",end_rlen);
02969 #endif
02970
02971 }
02972
02973 if((heuristics & TWO_POINT_MOVE)==TWO_POINT_MOVE)
02974 {
02975
02976 tpm_search:
02977 r_start=route[r].start;
02978 start_rlen=route[r].length;
02979 i=r_start;
02980 j=VRPH_MAX(next_array[i],0);
02981 while(i!=VRPH_DEPOT)
02982 {
02983 while(TPM.search(this,i,rules));
02984 i=VRPH_MAX(next_array[i],VRPH_DEPOT);
02985 }
02986
02987 end_rlen = route[r].length;
02988 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
02989 goto tpm_search;
02990
02991 #if CLEAN_DEBUG
02992 printf("CLEAN::TPM end_rlen=%f\n",end_rlen);
02993 #endif
02994
02995
02996 }
02997
02998 if((heuristics & TWO_OPT)==TWO_OPT)
02999 {
03000
03001 to_search:
03002 r_start=route[r].start;
03003 start_rlen=route[r].length;
03004 i=r_start;
03005 j=VRPH_MAX(next_array[i],0);
03006 while(i!=VRPH_DEPOT)
03007 {
03008 while(TO.search(this,i,rules));
03009 i=VRPH_MAX(next_array[i],0);
03010 }
03011
03012 end_rlen = route[r].length;
03013 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
03014 goto to_search;
03015
03016 #if CLEAN_DEBUG
03017 printf("CLEAN::TO end_rlen=%f\n",end_rlen);
03018 #endif
03019
03020
03021 }
03022
03023 if((heuristics & OR_OPT)==OR_OPT)
03024 {
03025
03026 or_search:
03027 r_start=route[r].start;
03028 start_rlen=route[r].length;
03029 i=r_start;
03030 j=VRPH_MAX(next_array[i],0);
03031 while(i!=VRPH_DEPOT)
03032 {
03033 OrO.search(this,i,3,rules);
03034 i=VRPH_MAX(next_array[i],0);
03035 }
03036
03037 i=r_start;
03038 j=VRPH_MAX(next_array[i],0);
03039 while(i!=VRPH_DEPOT)
03040 {
03041 OrO.search(this,i,2,rules);
03042 i=VRPH_MAX(next_array[i],0);
03043 }
03044 end_rlen = route[r].length;
03045 if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
03046 goto or_search;
03047
03048 #if CLEAN_DEBUG
03049 printf("CLEAN::TO end_rlen=%f\n",end_rlen);
03050 #endif
03051
03052 }
03053
03054 if((heuristics & THREE_OPT)==THREE_OPT)
03055 {
03056
03057
03058 while(ThreeO.route_search(this,r,rules))
03059
03060 end_val = route[r].length;
03061 #if CLEAN_DEBUG
03062 printf("CLEAN::3O end_val=%f\n",end_val);
03063 #endif
03064
03065
03066 }
03067
03068 end_val= route[r].length;
03069
03070 #if CLEAN_DEBUG
03071 printf("CLEAN::final end_val=%f\n",end_val);
03072 #endif
03073
03074
03075 if(VRPH_ABS(start_val-end_val)>VRPH_EPSILON)
03076 goto start_improving;
03077 else
03078 return;
03079
03080 }
03081
03082 bool VRP::before(int a, int b)
03083 {
03089
03090 int i;
03091
03092 if(a==VRPH_DEPOT || b==VRPH_DEPOT)
03093 report_error("%s: before called with VRPH_DEPOT\n",__FUNCTION__);
03094
03095 if(route_num[a]!=route_num[b])
03096 {
03097 fprintf(stderr,"Ordering error: before called with %d and %d not in the same route!\n",a,b);
03098 report_error("%s: differnet routes\n",__FUNCTION__);
03099 }
03100
03101 if(next_array[a]==b)
03102 return true;
03103 if(next_array[b]==a)
03104 return false;
03105
03106 i=a;
03107 while(i>0 && i!=b)
03108 i=next_array[i];
03109
03110
03111
03112
03113
03114 if(i==b)
03115 return true;
03116 else
03117 return false;
03118
03119 }
03120
03121
03122 void VRP::update(VRPMove *M)
03123 {
03127
03128 if(M->num_affected_routes==0)
03129 return;
03130
03131 int i;
03132
03133 for(i=0; i < M->num_affected_routes; i++)
03134 {
03135
03136 route[M->route_nums[i]].length = M->route_lens[i];
03137
03138
03139 route[M->route_nums[i]].load = M->route_loads[i];
03140
03141
03142 route[M->route_nums[i]].num_customers = M->route_custs[i];
03143 }
03144
03145
03146 total_route_length = M->new_total_route_length;
03147
03148
03149 total_number_of_routes = M->total_number_of_routes;
03150
03151 return;
03152
03153 }
03154
03155
03156
03157
03158 void VRP::compute_route_center(int r)
03159 {
03164
03165 int current_node = route[r].start;
03166 double total_x=0;
03167 double total_y=0;
03168
03169 while(current_node != VRPH_DEPOT)
03170 {
03171 total_x += nodes[current_node].x;
03172 total_y += nodes[current_node].y;
03173 current_node=VRPH_MAX(VRPH_DEPOT, next_array[current_node]);
03174 }
03175
03176 route[r].x_center = total_x / ((double) route[r].num_customers);
03177 route[r].y_center = total_y / ((double) route[r].num_customers);
03178
03179 }
03180 void VRP::find_neighboring_routes()
03181 {
03187
03188 int i,j;
03189 VRPNeighborElement **rd;
03190
03191
03192 normalize_route_numbers();
03193 for(i=1;i<=total_number_of_routes;i++)
03194 {
03195 compute_route_center(i);
03196 }
03197
03198
03199 rd = new VRPNeighborElement *[total_number_of_routes + 1];
03200 rd[0]= new VRPNeighborElement[(total_number_of_routes +1) * (total_number_of_routes + 1)];
03201 for(i=1;i<total_number_of_routes+1;i++)
03202 rd[i]=rd[i-1] + (total_number_of_routes+1);
03203
03204
03205 for(i=1;i<=this->total_number_of_routes;i++)
03206 {
03207 rd[i][0].position=VRP_INFINITY;
03208 rd[i][0].val=VRP_INFINITY;
03209 for(j=1;j<=this->total_number_of_routes;j++)
03210 {
03211 rd[i][j].position = j;
03212 rd[i][j].val = VRPDistance(VRPH_EUC_2D,route[i].x_center, route[i].y_center,
03213 route[j].x_center,route[j].y_center);
03214 }
03215 }
03216
03217
03218
03219
03220 for(i=1;i<=total_number_of_routes;i++)
03221 qsort(rd[i],total_number_of_routes+1,sizeof(VRPNeighborElement),VRPNeighborCompare);
03222
03223
03224
03225 for(i=1;i<=total_number_of_routes;i++)
03226 {
03227 for(j=0;j<MAX_NEIGHBORING_ROUTES;j++)
03228 {
03229 route[i].neighboring_routes[j]=rd[i][j+1].position;
03230
03231 }
03232 }
03233
03234 delete [] rd[0];
03235 delete [] rd;
03236
03237
03238
03239 return;
03240 }
03241
03242 void VRP::capture_best_solution()
03243 {
03247
03248 if( (this->total_route_length < this->best_total_route_length) &&
03249 (VRPH_ABS(this->total_route_length - this->best_total_route_length) > VRPH_EPSILON) )
03250 {
03251 this->best_total_route_length=this->total_route_length;
03252 this->export_solution_buff(this->best_sol_buff);
03253
03254 }
03255
03256
03257 if(this->total_route_length < this->solution_wh->worst_obj ||
03258 this->solution_wh->num_sols < this->solution_wh->max_size)
03259 {
03260 VRPSolution this_sol(this->num_nodes);
03261
03262 this_sol.obj=this->total_route_length;
03263 this_sol.in_IP=false;
03264
03265
03266 this->export_canonical_solution_buff(this_sol.sol);
03267
03268 this->solution_wh->add_sol(&this_sol, 0);
03269
03270 }
03271
03272 return;
03273
03274
03275 }
03276
03277 void VRP::update_solution_wh()
03278 {
03282
03283 VRPSolution this_sol(this->num_nodes);
03284
03285 this_sol.obj=this->total_route_length;
03286 this_sol.in_IP=false;
03287
03288
03289 this->export_canonical_solution_buff(this_sol.sol);
03290
03291 this->solution_wh->add_sol(&this_sol, 0);
03292
03293
03294
03295 return ;
03296
03297
03298 }
03299
03300
03301 void VRP::update_route(int j, VRPRoute *R)
03302 {
03308
03309 int i, current;
03310 double st=0;
03311
03312 R->x[0]=this->nodes[VRPH_DEPOT].x;
03313 R->y[0]=this->nodes[VRPH_DEPOT].y;
03314 R->start=route[j].start;
03315 R->end=route[j].end;
03316 R->num_customers=route[j].num_customers;
03317 R->load=route[j].load;
03318 R->length=route[j].length;
03319 R->obj_val=this->total_route_length-this->total_service_time;
03320
03321
03322 if(R->start < R->end)
03323 {
03324 current=R->start;
03325
03326 R->ordering[0]=current;
03327 R->x[1]=this->nodes[current].x;
03328 R->y[1]=this->nodes[current].y;
03329 st+=this->nodes[current].service_time;
03330
03331 for(i=1; i<R->num_customers; i++)
03332 {
03333 current=this->next_array[current];
03334 st+=this->nodes[current].service_time;
03335 R->ordering[i]=current;
03336 R->x[i+1]=this->nodes[current].x;
03337 R->y[i+1]=this->nodes[current].y;
03338 }
03339
03340 R->total_service_time=st;
03341
03342 return;
03343 }
03344
03345
03346 current=R->end;
03347 R->ordering[0]=current;
03348 R->x[1]=this->nodes[current].x;
03349 R->y[1]=this->nodes[current].y;
03350 st+=this->nodes[current].service_time;
03351
03352
03353 for(i=1; i<R->num_customers; i++)
03354 {
03355 current=this->pred_array[current];
03356 st+=this->nodes[current].service_time;
03357 R->ordering[i]=current;
03358 R->x[i+1]=this->nodes[current].x;
03359 R->y[i+1]=this->nodes[current].y;
03360
03361 }
03362
03363 R->total_service_time=st;
03364
03365 return;
03366
03367 }
03368
03369
03370
03371 double VRP::split(double p)
03372 {
03380
03381 if(p>.5)
03382 report_error("%s: p must be less than .5\n");
03383
03384 int i,j,k;
03385 double beta;
03386 struct double_int *thetas;
03387
03388 thetas=new double_int[this->num_nodes];
03389
03390 for(i=1;i<=this->num_nodes;i++)
03391 {
03392 thetas[i-1].k=i;
03393 thetas[i-1].d=this->nodes[i].theta;
03394 }
03395
03396
03397 qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare);
03398
03399
03400
03401
03402
03403 for(;;)
03404 {
03405 beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta);
03406
03407
03408 k=0;
03409 for(j=0;j<this->num_nodes;j++)
03410 {
03411
03412 if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x))
03413 k++;
03414 }
03415
03416 if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes))))
03417 break;
03418 }
03419
03420
03421 for(i=1;i<=this->num_original_nodes;i++)
03422 {
03423 if(routed[i])
03424 {
03425 if(this->nodes[i].y >= tan(beta)*(this->nodes[i].x))
03426 this->eject_node(i);
03427 }
03428 }
03429
03430
03431 delete [] thetas;
03432 return beta;
03433 }
03434
03435
03436 int VRP::split_routes(double p, int **ejected_routes, double *t)
03437 {
03448
03449 if(p>.5)
03450 report_error("%s: p must be less than .5\n");
03451
03452 int i,j,k;
03453 double beta;
03454 struct double_int *thetas;
03455
03456 thetas=new double_int[this->num_nodes];
03457
03458 for(i=1;i<=this->num_nodes;i++)
03459 {
03460 thetas[i-1].k=i;
03461 thetas[i-1].d=this->nodes[i].theta;
03462 }
03463
03464
03465 qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare);
03466
03467
03468
03469
03470
03471 for(;;)
03472 {
03473 beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta);
03474
03475
03476 k=0;
03477 for(j=0;j<this->num_nodes;j++)
03478 {
03479
03480 if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x))
03481 k++;
03482 }
03483
03484 if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes))))
03485 {
03486 break;
03487 }
03488 }
03489
03490
03491
03492 int num_ejected=0;
03493 for(i=1;i<=this->total_number_of_routes;i++)
03494 {
03495 int start=this->route[i].start;
03496 int current=start;
03497 bool will_eject=true;
03498
03499 while(current!=VRPH_DEPOT)
03500 {
03501
03502 if(this->nodes[current].y <= tan(beta)*(this->nodes[current].x))
03503 {
03504
03505 will_eject=false;
03506 break;
03507 }
03508
03509 current=VRPH_MAX(VRPH_DEPOT,this->next_array[current]);
03510 }
03511 if(will_eject)
03512 {
03513
03514
03515 this->eject_route(i, ejected_routes[num_ejected]);
03516 num_ejected++;
03517
03518 }
03519
03520 }
03521
03522 delete [] thetas;
03523
03524 *t=beta;
03525 return num_ejected;
03526 }
03527
03528
03529 void VRP::fix_edge(int start, int end)
03530 {
03538
03539 this->fixed[start][end]=true;
03540 this->fixed[end][start]=true;
03541
03542
03543 if(start==VRPH_DEPOT)
03544 {
03545 this->fixed[dummy_index][start]=true;
03546 this->fixed[start][dummy_index]=true;
03547 }
03548
03549 if(end==VRPH_DEPOT)
03550 {
03551 this->fixed[dummy_index][end]=true;
03552 this->fixed[end][dummy_index]=true;
03553 }
03554
03555
03556 }
03557
03558 void VRP::unfix_edge(int start, int end)
03559 {
03563
03564 if(this->fixed[start][end])
03565 report_error("%s: Edge %d-%d is not already fixed!\n",__FUNCTION__,start,end);
03566
03567 this->fixed[start][end]=false;
03568 this->fixed[end][start]=false;
03569
03570
03571 if(start==VRPH_DEPOT)
03572 {
03573 this->fixed[dummy_index][start]=false;
03574 this->fixed[start][dummy_index]=false;
03575 }
03576
03577 if(end==VRPH_DEPOT)
03578 {
03579 this->fixed[dummy_index][end]=false;
03580 this->fixed[end][dummy_index]=false;
03581 }
03582
03583
03584 }
03585
03586 void VRP::unfix_all()
03587 {
03591
03592 for(int i=0;i<=matrix_size;i++)
03593 for(int j=0;j<=matrix_size;j++)
03594 fixed[i][j]=false;
03595 }
03596
03597 void VRP::fix_string(int *node_string, int k)
03598 {
03605
03606
03607 for(int i=0;i<k-1;i++)
03608 this->fix_edge(node_string[i], node_string[i+1]);
03609
03610
03611 }
03612
03613
03614
03615
03616
03617 int VRP::read_fixed_edges(const char *filename)
03618 {
03631
03632
03633 int i, a, b, k;
03634 FILE *in;
03635
03636 if( (in=fopen(filename,"r"))==NULL)
03637 {
03638 fprintf(stderr,"Error opening %s for reading\n",filename);
03639 report_error("%s: Bad file in read_fixed_edges\n",__FUNCTION__);
03640 }
03641
03642 fscanf(in,"%d\n",&k);
03643
03644 for(i=0;i<k;i++)
03645 {
03646 fscanf(in,"%d %d\n",&a, &b);
03647 if(a<0 || b<0 || a>this->num_original_nodes || b>this->num_original_nodes)
03648 {
03649 fprintf(stderr,"Tried to fix edge %d-%d\n",a,b);
03650 report_error("%s: Error in read_fixed_edges\n",__FUNCTION__);
03651 }
03652 this->fix_edge(a,b);
03653 }
03654
03655
03656 return k;
03657 }
03658 void VRP::list_fixed_edges(int *fixed_list)
03659 {
03666
03667 int pos=0;
03668 int current=VRPH_DEPOT;
03669 int next=VRPH_ABS(next_array[current]);
03670
03671 while(next!=VRPH_DEPOT)
03672 {
03673 if(fixed[current][next])
03674 {
03675 fixed_list[pos]=current;
03676 fixed_list[pos+1]=next;
03677 pos+=2;
03678 }
03679
03680 current=next;
03681 next=next_array[current];
03682
03683 if(next<0)
03684 {
03685 if(fixed[current][VRPH_DEPOT])
03686 {
03687 fixed_list[pos]=current;
03688 fixed_list[pos+1]=VRPH_DEPOT;
03689 pos+=2;
03690 }
03691
03692 if(fixed[VRPH_DEPOT][-next])
03693 {
03694 fixed_list[pos]=VRPH_DEPOT;
03695 fixed_list[pos+1]=-next;
03696 pos+=2;
03697 }
03698
03699 next=-next;
03700 }
03701 }
03702
03703
03704
03705 }
03706 void VRP::perturb_locations(double c)
03707 {
03714
03715 int i, pre, post;
03716 double v,theta;
03717
03718
03719 this->export_solution_buff(this->current_sol_buff);
03720
03721 for(i=1;i<=this->num_nodes;i++)
03722 {
03723 pre= VRPH_MAX(VRPH_DEPOT,this->pred_array[i]);
03724 post=VRPH_MAX(VRPH_DEPOT,this->next_array[i]);
03725
03726 v= (this->d[pre][i]+this->d[i][post]) - (this->nodes[i].service_time + this->nodes[post].service_time);
03727 v= v*c;
03728
03729
03730 theta= VRPH_PI*lcgrand(8);
03731
03732 this->nodes[i].x += v*cos(theta);
03733 this->nodes[i].y += v*sin(theta);
03734
03735 }
03736
03737
03738
03739 this->create_distance_matrix(this->edge_weight_type);
03740
03741
03742 this->import_solution_buff(this->current_sol_buff);
03743
03744 return;
03745
03746 }
03747 void VRP::add_route(int *route_buff)
03748 {
03754
03755
03756
03757
03758
03759
03760 int *temp_buff;
03761
03762 this->verify_routes("Before adding route\n");
03763
03764 temp_buff=new int[this->num_original_nodes+2];
03765 this->export_solution_buff(temp_buff);
03766
03767
03768 int old_num=this->num_nodes;
03769 int i=0;
03770 while(route_buff[i]!=-1)
03771 {
03772 temp_buff[old_num+1+i]=route_buff[i];
03773 if(i==0)
03774 temp_buff[old_num+1+i]=-temp_buff[old_num+1+i];
03775
03776 temp_buff[0]++;
03777 i++;
03778 }
03779
03780 temp_buff[old_num+1+i]=VRPH_DEPOT;
03781
03782
03783 this->import_solution_buff(temp_buff);
03784
03785 this->verify_routes("After adding route\n");
03786
03787 delete [] temp_buff;
03788
03789 }
03790 void VRP::append_route(int *sol_buff, int *route_buff)
03791 {
03798
03799 int i,j,current_num;
03800
03801 current_num=sol_buff[0];
03802
03803
03804
03805 j=0;
03806 while(route_buff[j]!=-1)
03807 j++;
03808
03809
03810
03811 sol_buff[0]+=j;
03812
03813 sol_buff[current_num+1]=-route_buff[0];
03814 for(i=1;i<j;i++)
03815 sol_buff[current_num+1+i]=route_buff[i];
03816
03817
03818
03819 sol_buff[current_num+1+j]=VRPH_DEPOT;
03820
03821 return;
03822
03823
03824 }
03825
03826 int VRP::intersect_solutions(int *new_sol, int **route_list, int *sol1, int *sol2, int min_routes)
03827 {
03838
03839 int i,j,k,m,num_routes;
03840 int *rnums;
03841
03842 rnums=new int[this->num_original_nodes+1];
03843
03844
03845 j = this->find_common_routes(sol1, sol2, rnums);
03846
03847 if(j==0)
03848 {
03849
03850 memcpy(new_sol,sol1, (this->num_original_nodes+2)*sizeof(int));
03851 delete[] rnums;
03852 return 0;
03853 }
03854
03855
03856
03857 this->import_solution_buff(sol1);
03858 num_routes=this->total_number_of_routes;
03859 k=0;
03860 for(i=0;i<j;i++)
03861 {
03862
03863 route_list[i][0]=this->route[rnums[i]].start;
03864
03865 m=0;
03866 while(route_list[i][m]!=this->route[rnums[i]].end)
03867 {
03868 m++;
03869 route_list[i][m]=this->next_array[route_list[i][m-1]];
03870 }
03871
03872 m++;
03873 route_list[i][m]=-1;
03874
03875 k++;
03876 num_routes--;
03877 if(num_routes==min_routes)
03878 break;
03879 }
03880
03881
03882
03883
03884
03885 int *junk;
03886 junk=new int[this->num_original_nodes];
03887
03888 for(i=0;i<k;i++)
03889 this->eject_route(rnums[i], junk);
03890
03891
03892 this->export_canonical_solution_buff(new_sol);
03893
03894
03895 this->import_solution_buff(new_sol);
03896
03897 delete [] rnums;
03898 delete [] junk;
03899
03900 return k;
03901
03902 }
03903 bool VRP::osman_insert(int k, double alpha)
03904 {
03912
03913 int i,j,l,ll,m,mm, best_l, best_m;
03914 double savings, best_savings, ik,kj,ij;
03915 Postsert postsert;
03916 Presert presert;
03917 VRPMove M;
03918
03919 if(k==VRPH_DEPOT)
03920 return false;
03921
03922 i=VRPH_MAX(this->pred_array[k],VRPH_DEPOT);
03923 j=VRPH_MAX(this->next_array[k],VRPH_DEPOT);
03924
03925
03926 ik=this->d[i][k];
03927 kj=this->d[k][j];
03928 ij=this->d[i][j];
03929
03930 best_savings=VRP_INFINITY;
03931 best_l=-1;
03932 best_m=-1;
03933
03934 l=VRPH_DEPOT;
03935 m=abs(this->next_array[l]);
03936 while(m!=VRPH_DEPOT)
03937 {
03938 if(m>0)
03939 {
03940 savings = (ik + kj - this->d[l][m]) - alpha*(this->d[l][k] + this->d[k][m] - ij);
03941 if(savings<best_savings && l!=i && m!=j)
03942 {
03943
03944 if(presert.evaluate(this,k,m,&M)==true)
03945 {
03946 best_savings=savings;
03947 best_l=l;
03948 best_m=m;
03949 }
03950 }
03951 }
03952 else
03953 {
03954
03955
03956
03957 mm=VRPH_DEPOT;
03958 savings = (ik + kj - this->d[l][mm]) - alpha*(this->d[l][k] + this->d[k][mm] - ij);
03959 if(savings<best_savings && l!=i && mm!=j)
03960 {
03961
03962 if(postsert.evaluate(this,k,l,&M)==true)
03963 {
03964 best_savings=savings;
03965 best_l=l;
03966 best_m=mm;
03967 }
03968 }
03969
03970
03971 ll=VRPH_DEPOT;
03972 m=abs(m);
03973 savings = (ik + kj - this->d[ll][m]) - alpha*(this->d[ll][k] + this->d[k][m] - ij);
03974 if(savings<best_savings && ll!=i && m!=j)
03975 {
03976
03977 if(presert.evaluate(this,k,m,&M)==true)
03978 {
03979 best_savings=savings;
03980 best_l=ll;
03981 best_m=m;
03982 }
03983 }
03984 }
03985
03986 l=m;
03987 m=this->next_array[l];
03988 }
03989
03990
03991
03992
03993 if(best_savings==VRP_INFINITY)
03994 {
03995
03996 return false;
03997 }
03998
03999 if(best_l!=VRPH_DEPOT)
04000 {
04001 if(postsert.move(this,k,best_l)==false)
04002 report_error("%s: postsert.move is false\n",__FUNCTION__);
04003 }
04004 else
04005 {
04006 if(presert.move(this,k,best_m)==false)
04007 report_error("%s: presert.move is false\n",__FUNCTION__);
04008 }
04009
04010 return true;
04011
04012 }
04013
04014 int VRP::osman_perturb(int num, double alpha)
04015 {
04021
04022 int tot=0;
04023 int k;
04024 int num_attempts=0;
04025
04026 while(tot<num)
04027 {
04028
04029 k=VRPH_MAX((int)(this->num_original_nodes * lcgrand(10)),1);
04030
04031 if(routed[k])
04032 {
04033
04034
04035 if(this->osman_insert(k,alpha)==true)
04036 tot++;
04037
04038 num_attempts++;
04039 if(num_attempts>2*this->num_original_nodes)
04040 return tot;
04041 }
04042 }
04043
04044
04045
04046 return num;
04047
04048
04049 }
04050
04051 bool VRP::check_fixed_edges(const char *message)
04052 {
04058
04059
04060
04061 int i,j;
04062
04063 for(i=0;i<=this->num_original_nodes;i++)
04064 {
04065 for(j=0;j<=this->num_original_nodes;j++)
04066 {
04067 if(this->fixed[i][j])
04068 {
04069
04070 if(i!=VRPH_DEPOT)
04071 {
04072 if(VRPH_MAX(this->next_array[i],VRPH_DEPOT)!=j && VRPH_MAX(this->pred_array[i],VRPH_DEPOT)!=j)
04073 {
04074 fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j);
04075 fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[i],VRPH_DEPOT),i,
04076 VRPH_MAX(this->next_array[i],VRPH_DEPOT));
04077 fprintf(stderr,"%s", message);
04078
04079 if(this->fixed[j][i])
04080 fprintf(stderr,"%d-%d also fixed\n",j,i);
04081 else
04082 fprintf(stderr,"%d-%d NOT fixed!\n",j,i);
04083
04084 return false;
04085 }
04086 }
04087
04088
04089 if(j!=VRPH_DEPOT)
04090 {
04091 if(VRPH_MAX(this->next_array[j],VRPH_DEPOT)!=i && VRPH_MAX(this->pred_array[j],VRPH_DEPOT)!=i)
04092 {
04093 fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j);
04094 fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[j],VRPH_DEPOT),j,
04095 VRPH_MAX(this->next_array[j],VRPH_DEPOT));
04096 fprintf(stderr,"%s", message);
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 return false;
04102
04103 }
04104 }
04105 }
04106
04107 }
04108
04109 }
04110
04111 return true;
04112 }
04113 int VRP::find_common_routes(int *sol1, int *sol2, int *route_nums)
04114 {
04120
04121 int *h11, *h12, *h21, *h22;
04122 double *L1, *L2;
04123 int i, j, cnt, r1, r2;
04124 VRPRoute rte(num_nodes);
04125
04126
04127 this->import_solution_buff(sol1);
04128
04129 r1=this->total_number_of_routes;
04130 h11=new int[r1+1]; h12=new int[r1+1];
04131 L1=new double[r1+1];
04132 for(i=1;i<=r1;i++)
04133 {
04134
04135 this->update_route(i,&rte);
04136
04137 h11[i]=rte.hash(SALT_1);
04138 h12[i]=rte.hash(SALT_2);
04139 L1[i]=rte.length;
04140 }
04141
04142
04143 this->import_solution_buff(sol2);
04144
04145 r2=this->total_number_of_routes;
04146 h21=new int[r2+1]; h22=new int[r2+1];
04147 L2=new double[r2+1];
04148 for(i=1;i<=r2;i++)
04149 {
04150 this->update_route(i,&rte);
04151 h21[i]=rte.hash(SALT_1);
04152 h22[i]=rte.hash(SALT_2);
04153 L2[i]=rte.length;
04154 }
04155
04156
04157
04158
04159 cnt=0;
04160 for(i=1;i<=r1;i++)
04161 {
04162 for(j=1;j<=r2;j++)
04163 {
04164 if(h11[i]==h21[j] && h12[i]==h22[j] && VRPH_ABS(L1[i]-L2[j])<VRPH_EPSILON)
04165 {
04166
04167 route_nums[cnt]=i;
04168 cnt++;
04169 }
04170 }
04171 }
04172
04173
04174 delete [] h11; delete [] h12; delete [] h21; delete [] h22;
04175 delete [] L1; delete [] L2;
04176 return cnt;
04177
04178 }
04179
04180 void VRP::set_daily_demands(int day)
04181 {
04187
04188 int i;
04189
04190 if(day>0)
04191 {
04192 for(i=0;i<=this->num_original_nodes;i++)
04193 {
04194 if(this->nodes[i].daily_demands[day]>=0)
04195 this->nodes[i].demand=this->nodes[i].daily_demands[day];
04196 else
04197 this->nodes[i].demand=-1;
04198
04199 }
04200 }
04201 else
04202 {
04203
04204 for(i=0;i<=this->num_original_nodes;i++)
04205 {
04206 int mean_demand=0;
04207 int k=0;
04208 for(int j=1;j<=this->num_days;j++)
04209 {
04210 if(this->nodes[i].daily_demands[j]>=0)
04211 {
04212 mean_demand+=this->nodes[i].daily_demands[j];
04213 k++;
04214 }
04215 }
04216 if(k>0)
04217 this->nodes[i].demand=(int)((double)mean_demand/(double)k);
04218 else
04219 this->nodes[i].demand=-1;
04220
04221 }
04222
04223 }
04224 }
04225
04226 void VRP::set_daily_service_times(int day)
04227 {
04232
04233 int i;
04234
04235 for(i=0;i<=this->num_original_nodes;i++)
04236 {
04237 this->nodes[i].service_time=this->nodes[i].daily_service_times[day];
04238 }
04239
04240
04241
04242 this->create_distance_matrix(this->edge_weight_type);
04243 }
04244
04245
04246 void VRP::update_arrival_times()
04247 {
04251
04252 int routenum, next, current;
04253 double t;
04254
04255
04256
04257 for(int i=1;i<=this->num_original_nodes;i++)
04258 this->nodes[i].arrival_time=-1;
04259
04260
04261 for(int i=1;i<=this->num_original_nodes;i++)
04262 {
04263
04264 if(routed[i])
04265 {
04266
04267 routenum=this->route_num[i];
04268
04269
04270 current=this->route[routenum].start;
04271 t=this->d[VRPH_DEPOT][current];
04272
04273 while(current!=i)
04274 {
04275 next=VRPH_MAX(VRPH_DEPOT, this->next_array[current]);
04276 t+=this->d[current][next];
04277 current=next;
04278
04279 }
04280
04281
04282
04283 t-=this->nodes[i].service_time;
04284 this->nodes[i].arrival_time=t;
04285 }
04286 }
04287
04288 }
04289
04290 bool VRP::check_tabu_status(VRPMove *M, int *old_sol)
04291 {
04299
04300
04301 #if VRPH_TABU_DEBUG
04302 printf("Checking tabu status\n");
04303 #endif
04304
04305
04306 int i,j;
04307
04308
04309 if(M->num_affected_routes>1)
04310 {
04311 for(i=0;i<M->num_affected_routes;i++)
04312 {
04313 if(M->route_custs[i]==0)
04314 return true;
04315 }
04316 }
04317
04318 VRPRoute r(this->num_nodes);
04319 int num_tabu_routes=0;
04320
04321 #if VRPH_TABU_DEBUG
04322 printf("Checking tabu status of %d routes\n",M->num_affected_routes);
04323 #endif
04324
04325 for(i=0;i<M->num_affected_routes;i++)
04326 {
04327 this->update_route(M->route_nums[i],&r);
04328 r.hash_val = r.hash(SALT_1);
04329 r.hash_val2 = r.hash(SALT_2);
04330
04331 for(j=0;j<this->tabu_list->num_entries;j++)
04332 {
04333 if(r.hash_val==this->tabu_list->hash_vals1[j] && r.hash_val2==this->tabu_list->hash_vals2[j])
04334
04335 num_tabu_routes++;
04336 }
04337
04338 }
04339
04340
04341 if(num_tabu_routes>0)
04342 {
04343
04344 this->import_solution_buff(old_sol);
04345 #if VRPH_TABU_DEBUG
04346 printf("Move is Tabu! Reverting to old solution\n");
04347 #endif
04348 return false;
04349 }
04350
04351 #if VRPH_TABU_DEBUG
04352 printf("Move is not Tabu. Updating list of %d routes\n",M->num_affected_routes);
04353 #endif
04354
04355
04356 for(i=0;i<M->num_affected_routes;i++)
04357 {
04358 this->update_route(M->route_nums[i],&r);
04359
04360 this->tabu_list->update_list(&r);
04361 }
04362
04363 return true;
04364
04365 }
04366
04367
04368 void VRP::print_stats()
04369 {
04374
04375 printf(" (Moves, Evaluations)\n");
04376 printf(" One Point Move: (%010d, %010d)\n",
04377 this->num_moves[ONE_POINT_MOVE_INDEX], this->num_evaluations[ONE_POINT_MOVE_INDEX]);
04378 printf(" Two Point Move: (%010d, %010d)\n",
04379 this->num_moves[TWO_POINT_MOVE_INDEX], this->num_evaluations[TWO_POINT_MOVE_INDEX]);
04380 printf(" Three Point Move: (%010d, %010d)\n",
04381 this->num_moves[THREE_POINT_MOVE_INDEX], this->num_evaluations[THREE_POINT_MOVE_INDEX]);
04382 printf(" Two-opt Move: (%010d, %010d)\n",
04383 this->num_moves[TWO_OPT_INDEX], this->num_evaluations[TWO_OPT_INDEX]);
04384 printf(" Three-opt Move: (%010d, %010d)\n",
04385 this->num_moves[THREE_OPT_INDEX], this->num_evaluations[THREE_OPT_INDEX]);
04386 printf(" Or-opt Move: (%010d, %010d)\n",
04387 this->num_moves[OR_OPT_INDEX], this->num_evaluations[OR_OPT_INDEX]);
04388 printf("Cross-Exchange Move: (%010d, %010d)\n\n",
04389 this->num_moves[CROSS_EXCHANGE_INDEX], this->num_evaluations[CROSS_EXCHANGE_INDEX]);
04390
04391 return;
04392
04393 }
04394
04395
04396 void VRP::reset()
04397 {
04401
04402 this->solution_wh->liquidate();
04403 for(int j=1;j<=this->num_original_nodes;j++)
04404 this->routed[j]=false;
04405 }
04406