VRPH  1.0
src/VRPSolvers.cpp
Go to the documentation of this file.
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 
00015 double VRP::RTR_solve(int heuristics, int intensity, int max_stuck, int max_perturbs,
00016                       double dev, int nlist_size, int perturb_type, int accept_type, bool verbose)
00017 {
00025 
00026     // Make sure accept_type is either VRPH_BEST_ACCEPT or VRPH_FIRST_ACCEPT - matters only
00027     // for the downhill phase as we use VRPH_LI_ACCEPT in the diversification phase
00028 
00029     if(accept_type!=VRPH_BEST_ACCEPT && accept_type!=VRPH_FIRST_ACCEPT)
00030         report_error("%s: accept_type must be VRPH_BEST_ACCEPT or VRPH_FIRST_ACCEPT\n");
00031 
00032     int ctr, n, j,  i,  R, random, fixed, neighbor_list, objective, tabu;
00033 
00034     random=fixed=neighbor_list=0;
00035 
00036     if(heuristics & VRPH_RANDOMIZED)
00037         random=VRPH_RANDOMIZED;
00038 
00039     if(heuristics & VRPH_FIXED_EDGES)
00040         fixed=VRPH_FIXED_EDGES;
00041 
00042     if(heuristics & VRPH_USE_NEIGHBOR_LIST)
00043         neighbor_list=VRPH_USE_NEIGHBOR_LIST;
00044 
00045     objective=VRPH_SAVINGS_ONLY;
00046     // default strategy
00047 
00048     if(heuristics & VRPH_MINIMIZE_NUM_ROUTES)
00049         objective=VRPH_MINIMIZE_NUM_ROUTES;
00050 
00051 
00052     if(heuristics & VRPH_TABU)
00053     {
00054         tabu=VRPH_TABU; // We will use a primitive Tabu Search in the uphill phase
00055         // Clear the tabu list
00056         this->tabu_list->empty();
00057     }
00058     else
00059         tabu=0;
00060 
00061     n=num_nodes;
00062 
00063     // Define the heuristics we will use
00064 
00065     OnePointMove OPM;
00066     TwoPointMove TPM;
00067     TwoOpt         TO;
00068     OrOpt         OR;
00069     ThreeOpt     ThreeO;
00070     CrossExchange    CE;
00071     ThreePointMove ThreePM;
00072 
00073     double start_val;
00074     int *perm;
00075     perm=new int[this->num_nodes];
00076 
00077 
00078     j=VRPH_ABS(this->next_array[VRPH_DEPOT]);
00079     for(i=0;i<this->num_nodes;i++)
00080     {
00081         perm[i]=j;
00082         if(!routed[j])
00083             report_error("%s: Unrouted node in solution!!\n");
00084 
00085         j=VRPH_ABS(this->next_array[j]);
00086     }
00087     if(j!=VRPH_DEPOT)
00088         report_error("%s: VRPH_DEPOT is not last node in solution!!\n");
00089 
00090 
00091     int rules;
00092 
00093     // Set the neighbor list size used in the improvement search
00094     neighbor_list_size=VRPH_MIN(nlist_size, this->num_nodes);
00095 
00096     // Set the deviation
00097     deviation=dev;
00098 
00099     int num_perturbs=0;
00100 
00101     record=this->total_route_length;
00102     this->best_total_route_length=this->total_route_length;
00103     this->export_solution_buff(this->current_sol_buff);
00104     this->export_solution_buff(this->best_sol_buff);
00105 
00106     normalize_route_numbers();
00107 
00108     ctr=0;
00109 
00110 
00111 uphill:
00112     // Start an uphill phase using the following "rules":
00113     double beginning_best=this->best_total_route_length;
00114     rules=VRPH_LI_ACCEPT+VRPH_RECORD_TO_RECORD+objective+random+fixed+neighbor_list+tabu;
00115 
00116     if(verbose)
00117         printf("Uphill starting at %5.2f\n",this->total_route_length);
00118     
00119     for(int k=1;k<intensity;k++)
00120     {
00121         start_val=total_route_length;
00122 
00123         if(heuristics & ONE_POINT_MOVE)
00124         {
00125             if(random)
00126                 random_permutation(perm, this->num_nodes);
00127 
00128             for(i=1;i<=n;i++)
00129             {
00130 #if FIXED_DEBUG
00131                 if(fixed && !check_fixed_edges("Before 1PM\n"))
00132                     fprintf(stderr,"Error before OPM search(%d)\n",perm[i-1]);
00133 #endif
00134                 OPM.search(this,perm[i-1],rules);
00135 
00136 #if FIXED_DEBUG
00137                 if(fixed && !check_fixed_edges("After 1PM\n"))
00138                 {
00139                     fprintf(stderr,"Error after OPM search(%d)\n",perm[i-1]);
00140                     this->show_route(this->route_num[perm[i-1]]);
00141                 }
00142 #endif
00143             }
00144         }
00145 
00146 
00147         if(heuristics & TWO_POINT_MOVE)
00148         {
00149             if(random)
00150                 random_permutation(perm, this->num_nodes);
00151 
00152             for(i=1;i<=n;i++)    
00153                 TPM.search(this,perm[i-1],rules + VRPH_INTER_ROUTE_ONLY);
00154 
00155             //check_fixed_edges("After 2PM\n");
00156 
00157         }
00158 
00159 
00160         if(heuristics & THREE_POINT_MOVE)
00161         {
00162             if(random)
00163                 random_permutation(perm, this->num_nodes);
00164 
00165             for(i=1;i<=n;i++)    
00166                 ThreePM.search(this,perm[i-1],rules + VRPH_INTER_ROUTE_ONLY);
00167 
00168             //check_fixed_edges("After 3PM\n");
00169 
00170         }
00171 
00172 
00173 
00174         if(heuristics & TWO_OPT)
00175         {
00176             if(random)
00177                 random_permutation(perm, this->num_nodes);
00178 
00179             for(i=1;i<=n;i++)    
00180                 TO.search(this,perm[i-1],rules);
00181 
00182             //check_fixed_edges("After TO\n");
00183 
00184 
00185         }        
00186 
00187 
00188         if(heuristics & OR_OPT)
00189         {
00190             if(random)
00191                 random_permutation(perm, this->num_nodes);
00192 
00193             for(i=1;i<=n;i++)    
00194                 OR.search(this,perm[i-1],4,rules);
00195 
00196             for(i=1;i<=n;i++)    
00197                 OR.search(this,perm[i-1],3,rules);
00198 
00199             for(i=1;i<=n;i++)    
00200                 OR.search(this,perm[i-1],2,rules);
00201 
00202             //check_fixed_edges("After OR\n");
00203 
00204         }
00205 
00206         if(heuristics & THREE_OPT)
00207         {
00208             normalize_route_numbers();
00209             R=total_number_of_routes;
00210 
00211             for(i=1; i<=R; i++)    
00212                 ThreeO.route_search(this,i,rules-neighbor_list);
00213 
00214             //check_fixed_edges("After 3O\n");
00215 
00216         }
00217 
00218         if(heuristics & CROSS_EXCHANGE)
00219         {
00220             normalize_route_numbers();
00221             this->find_neighboring_routes();
00222             R=total_number_of_routes;
00223 
00224             for(i=1; i<=R-1; i++)    
00225             {
00226                 for(j=0;j<1;j++)
00227                     CE.route_search(this,i, route[i].neighboring_routes[j],rules-neighbor_list); 
00228             }
00229 
00230             //check_fixed_edges("After CE\n");
00231         }
00232     }
00233 
00234     if(total_route_length<record)
00235         record = total_route_length;
00236 
00237     if(verbose)
00238     {
00239         printf("Uphill complete\t(%d,%5.2f,%5.2f)\n",count_num_routes(),total_route_length, record);
00240         printf("# of recorded routes: %d[%d]\n",total_number_of_routes,count_num_routes());
00241 
00242     }
00243 
00244     if(this->best_total_route_length<beginning_best-VRPH_EPSILON)
00245     {
00246         if(verbose)
00247             printf("New best found in uphill!\n");
00248         // We found a new best solution during the uphill phase that might
00249         // now be "forgotten"!! I have seen this happen where it is never recovered
00250         // again, so we just import it and start the downhill phase with this solution...
00251         //this->import_solution_buff(this->best_sol_buff);
00252 
00253     }
00254 
00255 downhill:
00256 
00257     // Now enter a downhill phase
00258     double orig_val=total_route_length;
00259     if(verbose)
00260         printf("Downhill starting at %f (best=%f)\n",orig_val,this->best_total_route_length);
00261 
00262 
00263     if((heuristics & ONE_POINT_MOVE)|| (heuristics & KITCHEN_SINK) )
00264     {
00265         rules=VRPH_DOWNHILL+objective+random+fixed+neighbor_list+accept_type;
00266         for(;;)
00267         {
00268             // One Point Move
00269             start_val=total_route_length;
00270 
00271             if(random)
00272                 random_permutation(perm, this->num_nodes);
00273 
00274             for(i=1;i<=n;i++)
00275                 OPM.search(this,perm[i-1],rules );
00276 
00277 
00278             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00279                 break; 
00280 
00281         }
00282 
00283     }
00284 
00285 
00286 
00287     if((heuristics & TWO_POINT_MOVE) || (heuristics & KITCHEN_SINK) )
00288     {
00289         rules=VRPH_DOWNHILL+VRPH_INTER_ROUTE_ONLY+objective+random+fixed+neighbor_list+accept_type;
00290         for(;;)
00291         {
00292             // Two Point Move
00293             start_val=total_route_length;
00294 
00295             if(random)
00296                 random_permutation(perm, this->num_nodes);
00297 
00298             for(i=1;i<=n;i++)    
00299                 TPM.search(this,perm[i-1],rules);
00300 
00301             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00302                 break; 
00303 
00304         }
00305 
00306     }
00307 
00308 
00309 
00310     if((heuristics & TWO_OPT)|| (heuristics & KITCHEN_SINK) )
00311     {
00312         // Do inter-route first a la Li
00313         rules=VRPH_DOWNHILL+VRPH_INTER_ROUTE_ONLY+objective+random+fixed+neighbor_list+accept_type;
00314         for(;;)
00315         {
00316 
00317             start_val=total_route_length;
00318 
00319             if(random)
00320                 random_permutation(perm, this->num_nodes);
00321 
00322             for(i=1;i<=n;i++)    
00323                 TO.search(this,perm[i-1],rules);
00324 
00325             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00326                 break; 
00327         }
00328 
00329         // Now do both intra and inter
00330         rules=VRPH_DOWNHILL+objective+random+fixed+neighbor_list+accept_type;
00331 
00332         for(;;)
00333         {
00334 
00335             start_val=total_route_length;
00336 
00337             if(random)
00338                 random_permutation(perm, this->num_nodes);
00339 
00340             for(i=1;i<=n;i++)    
00341                 TO.search(this,perm[i-1],rules);
00342 
00343             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00344                 break; 
00345         }
00346     }
00347 
00348     if((heuristics & THREE_POINT_MOVE) || (heuristics & KITCHEN_SINK) )
00349     {
00350         rules=VRPH_DOWNHILL+VRPH_INTER_ROUTE_ONLY+objective+random+fixed+accept_type+neighbor_list;
00351         for(;;)
00352         {
00353             // Three Point Move
00354             start_val=total_route_length;
00355 
00356             if(random)
00357                 random_permutation(perm, this->num_nodes);
00358 
00359             for(i=1;i<=n;i++)    
00360                 ThreePM.search(this,perm[i-1],rules);
00361 
00362             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00363                 break; 
00364 
00365         }
00366     }
00367 
00368 
00369     if((heuristics & OR_OPT) || (heuristics & KITCHEN_SINK))
00370     {
00371 
00372         rules=VRPH_DOWNHILL+ objective +random +fixed + accept_type + neighbor_list;
00373 
00374         for(;;)
00375         {
00376             // OrOpt
00377             start_val=total_route_length;
00378             if(random)
00379                 random_permutation(perm, this->num_nodes);
00380 
00381             for(i=1;i<=n;i++)    
00382                 OR.search(this,perm[i-1],4,rules);
00383             for(i=1;i<=n;i++)
00384                 OR.search(this,perm[i-1],3,rules);
00385             for(i=1;i<=n;i++)
00386                 OR.search(this,perm[i-1],2,rules);
00387 
00388 
00389             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00390                 break; 
00391         }
00392     }
00393 
00394     if((heuristics & THREE_OPT) || (heuristics & KITCHEN_SINK) )
00395     {
00396         normalize_route_numbers();
00397         R= total_number_of_routes;
00398         rules=VRPH_DOWNHILL+objective+VRPH_INTRA_ROUTE_ONLY+ random +fixed + accept_type;
00399         for(;;)
00400         {
00401             // 3OPT
00402             start_val=total_route_length;
00403 
00404             for(i=1;i<=R;i++)    
00405                 ThreeO.route_search(this,i,rules);
00406 
00407             if(VRPH_ABS(total_route_length-start_val)<VRPH_EPSILON)
00408                 break; 
00409         }
00410     }
00411 
00412 
00413     if( (heuristics & CROSS_EXCHANGE) )
00414     {
00415         normalize_route_numbers();
00416         this->find_neighboring_routes();
00417         R=total_number_of_routes;
00418 
00419         rules=VRPH_DOWNHILL+objective+VRPH_INTRA_ROUTE_ONLY+ random +fixed + accept_type;
00420 
00421         for(i=1; i<=R-1; i++)    
00422         {
00423             for(j=0;j<=1;j++)
00424                 CE.route_search(this,i, route[i].neighboring_routes[j], rules); 
00425         }
00426     }
00427 
00428 
00429     // Repeat the downhill phase until we find no more improvements
00430     if(total_route_length<orig_val-VRPH_EPSILON)
00431         goto downhill;
00432     if(verbose)
00433         printf("Downhill complete: %5.2f[downhill started at %f] (%5.2f)\n",total_route_length,orig_val,
00434         this->best_total_route_length);
00435 
00436 
00437     if(total_route_length < record-VRPH_EPSILON)    
00438     {
00439         // New record - reset ctr
00440         ctr=1;
00441         record=total_route_length;
00442     }
00443     else
00444         ctr++;
00445 
00446     if(ctr<max_stuck)
00447         goto uphill;
00448 
00449     if(ctr==max_stuck)
00450     {
00451         if(num_perturbs<max_perturbs)
00452         {
00453             if(verbose)
00454                 printf("perturbing\n");
00455             if(perturb_type==VRPH_LI_PERTURB)
00456                 perturb();
00457             else
00458                 osman_perturb(VRPH_MAX(20,num_nodes/10),.5+lcgrand(20));
00459 
00460             // Reset record
00461             this->record=this->total_route_length;
00462             if(tabu)
00463                 this->tabu_list->empty();
00464 
00465             ctr=1;
00466             num_perturbs++;
00467             goto uphill;
00468         }
00469     }
00470 
00471 
00472     if(verbose)
00473     {
00474         if(has_service_times==false)
00475             printf("BEST OBJ:  %f\n",best_total_route_length);
00476         else
00477             printf("BEST OBJ:  %f\n",best_total_route_length-total_service_time);
00478     }
00479 
00480     delete [] perm; 
00481 
00482     // Import the best solution found
00483     this->import_solution_buff(best_sol_buff);
00484 
00485     if(has_service_times==false)
00486         return best_total_route_length;
00487     else
00488         return best_total_route_length-total_service_time;
00489 
00490 
00491 }
00492 
00493 double VRP::SA_solve(int heuristics, double start_temp, double cool_ratio,
00494                      int iters_per_loop, int num_loops, int nlist_size, bool verbose)    
00495 {
00501 
00502     this->temperature = start_temp;
00503     this->cooling_ratio = cool_ratio;
00504 
00505     int ctr, n, j,  i,  R, rules, random, fixed, neighbor_list, objective;
00506 
00507     if(heuristics & VRPH_RANDOMIZED)
00508         random=VRPH_RANDOMIZED;
00509     else
00510         random=0;
00511 
00512     if(heuristics & VRPH_FIXED_EDGES)
00513         fixed=VRPH_FIXED_EDGES;
00514     else
00515         fixed=0;
00516 
00517     if(heuristics & VRPH_USE_NEIGHBOR_LIST)
00518         neighbor_list=VRPH_USE_NEIGHBOR_LIST;
00519     else
00520         neighbor_list=0;
00521 
00522     objective=VRPH_SAVINGS_ONLY;
00523     // default strategy
00524 
00525     if(heuristics & VRPH_MINIMIZE_NUM_ROUTES)
00526         objective=VRPH_MINIMIZE_NUM_ROUTES;
00527     
00528     n=num_nodes;
00529 
00530     // The perm[] array will contain all the nodes in the current solution
00531     int *perm;
00532     perm=new int[this->num_nodes];
00533     j=VRPH_ABS(this->next_array[VRPH_DEPOT]);
00534     for(i=0;i<this->num_nodes;i++)
00535     {
00536         perm[i]=j;
00537         if(!routed[j])
00538             report_error("%s: Unrouted node in solution!!\n");
00539 
00540         j=VRPH_ABS(this->next_array[j]);
00541     }
00542     if(j!=VRPH_DEPOT)
00543         report_error("%s: VRPH_DEPOT is not last node in solution!!\n");
00544 
00545     // Define the heuristics we may use
00546 
00547     OnePointMove OPM;
00548     TwoPointMove TPM;
00549     TwoOpt         TO;
00550     OrOpt         OR;
00551     ThreeOpt     ThreeO;
00552     CrossExchange    CE;
00553     ThreePointMove ThreePM;
00554 
00555     double start_val;
00556 
00557     this->export_solution_buff(this->best_sol_buff);
00558     // We are assuming we have an existing solution
00559 
00560     // Set the neighbor list size used in the improvement search
00561     this->neighbor_list_size=VRPH_MIN(nlist_size, num_nodes);
00562 
00563     best_total_route_length=this->total_route_length;
00564     normalize_route_numbers();
00565 
00566     ctr=0;
00567 
00568     // The idea is to perform num_loops loops of num_iters iterations each.
00569     // For each iteration, run through the given heuristic operation at random
00570     // and perform a Simulated Annealing search.
00571 
00572     rules=VRPH_USE_NEIGHBOR_LIST+VRPH_FIRST_ACCEPT+VRPH_SIMULATED_ANNEALING+VRPH_SAVINGS_ONLY;
00573 
00574     double worst_obj=0;
00575     for(ctr=0;ctr<num_loops;ctr++)
00576     {
00577         if(verbose)
00578         {
00579             printf("\nctr=%d of %d, temp=%f, obj=%f (overall best=%f; worst=%f)\n",ctr,num_loops,
00580                 this->temperature, 
00581                 this->total_route_length,this->best_total_route_length,worst_obj);
00582             fflush(stdout);
00583         }
00584         // Reset worst_obj;
00585         worst_obj=0;
00586 
00587         // Cool it...
00588         this->temperature = this->cooling_ratio * this->temperature;
00589 
00590 
00591         for(int k=0; k < iters_per_loop; k++)
00592         {
00593             start_val=total_route_length;
00594             if(heuristics & THREE_OPT)
00595             {
00596                 rules=VRPH_SIMULATED_ANNEALING+VRPH_INTRA_ROUTE_ONLY+random+fixed+objective;
00597                 normalize_route_numbers();
00598                 R=total_number_of_routes;
00599                 for(i=1; i<=R; i++)    
00600                 {
00601                     ThreeO.route_search(this,i,rules);
00602                     if(this->total_route_length > worst_obj)
00603                         worst_obj=this->total_route_length;
00604                 }
00605 
00606             }
00607 
00608 
00609             if(heuristics & ONE_POINT_MOVE)
00610             {
00611                 rules=VRPH_SIMULATED_ANNEALING+neighbor_list+random+fixed+objective;
00612                 if(random)
00613                     random_permutation(perm, this->num_nodes);
00614 
00615                 for(i=1;i<=n;i++)    
00616                 {
00617 
00618                     OPM.search(this,perm[i-1],rules);
00619                     if(this->total_route_length > worst_obj)
00620                         worst_obj=this->total_route_length;
00621                 }
00622 
00623             }
00624 
00625 
00626 
00627             if(heuristics & TWO_POINT_MOVE)
00628             {
00629                 rules=VRPH_SIMULATED_ANNEALING+neighbor_list+random+fixed+objective;
00630                 if(random)
00631                     random_permutation(perm, this->num_nodes);
00632 
00633                 for(i=1;i<=n;i++)    
00634                 {
00635                     TPM.search(this,perm[i-1],rules);
00636                     if(this->total_route_length > worst_obj)
00637                         worst_obj=this->total_route_length;
00638                 }
00639 
00640 
00641             }
00642 
00643 
00644 
00645             if(heuristics & TWO_OPT)
00646             {
00647 
00648                 rules=VRPH_SIMULATED_ANNEALING+neighbor_list+random+fixed+objective;
00649                 if(random)
00650                     random_permutation(perm, this->num_nodes);
00651 
00652                 for(i=1;i<=n;i++)    
00653                 {
00654                     TO.search(this,perm[i-1],rules);
00655                     if(this->total_route_length > worst_obj)
00656                         worst_obj=this->total_route_length;
00657                 }
00658 
00659 
00660             }        
00661 
00662             if(heuristics & THREE_POINT_MOVE)
00663             {
00664                 rules=VRPH_SIMULATED_ANNEALING+VRPH_INTER_ROUTE_ONLY+neighbor_list+random+fixed+objective;
00665                 if(random)
00666                     random_permutation(perm, this->num_nodes);
00667 
00668 
00669                 for(i=1;i<=n;i++)    
00670                 {
00671                     ThreePM.search(this,perm[i-1],rules);
00672                     if(this->total_route_length > worst_obj)
00673                         worst_obj=this->total_route_length;
00674                 }
00675             }
00676 
00677             if(heuristics & OR_OPT)
00678             {
00679                 rules=VRPH_SIMULATED_ANNEALING+neighbor_list+random+fixed+objective;
00680                 if(random)
00681                     random_permutation(perm, this->num_nodes);
00682 
00683                 for(i=1;i<=n;i++)    
00684                 {
00685                     OR.search(this,perm[i-1],3,rules);
00686                     if(this->total_route_length > worst_obj)
00687                         worst_obj=this->total_route_length;
00688                 }
00689 
00690                 for(i=1;i<=n;i++)    
00691                 {
00692                     OR.search(this,perm[i-1],2,rules);
00693                     if(this->total_route_length > worst_obj)
00694                         worst_obj=this->total_route_length;
00695                 }
00696 
00697 
00698             }
00699 
00700             if(heuristics & CROSS_EXCHANGE)
00701             {
00702                 normalize_route_numbers();
00703                 this->find_neighboring_routes();
00704                 R=total_number_of_routes;
00705                 rules=VRPH_SIMULATED_ANNEALING+fixed+objective;
00706                 if(random)
00707                     random_permutation(perm, this->num_nodes);
00708 
00709 
00710                 for(i=1; i<=R-1; i++)    
00711                 {
00712                     for(j=0;j<=1;j++)
00713                     {
00714                         CE.route_search(this,i, route[i].neighboring_routes[j],rules);
00715                         if(this->total_route_length > worst_obj)
00716                             worst_obj=this->total_route_length;
00717                     }
00718                 }
00719             }            
00720         }
00721     }
00722 
00723     delete [] perm;
00724 
00725     // Restore the best sol
00726     this->import_solution_buff(this->best_sol_buff);
00727     // Now return the obj. function value
00728 
00729     if(has_service_times==false)
00730         return this->best_total_route_length;
00731     else
00732         return this->best_total_route_length-total_service_time;
00733 
00734 }
00735 
00736