VRPH  1.0
src/apps/vrp_glpk_sp.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 // Example of integrating VRPH with a Mixed Integer Programming solver to allow
00014 // generating solutions to the VRP by solving as a set partitioning problem.
00015 // We use the COIN-OSI interface in order to create "solver agnostic" code
00016 
00017 
00018 #include "OsiGlpkSolverInterface.hpp"
00019 #include "CoinPackedMatrix.hpp"
00020 #include "CoinPackedVector.hpp"
00021 #include "VRPH.h"
00022 
00023 #define MAX_ROUTES  50000
00024 // This is the maximum number of potentially different routes/columns
00025 // that we will bother to store
00026 
00027 // We allow only max_columns routes/columns in the IP
00028 // GLPK starts to struggle with more than 500 or so
00029 // Once we get more than this # of columns, we start fixing randomly
00030 // selected variables to 0 (GLPK column deletion routine seems buggy?)
00031 // CPLEX and other solvers can handle a couple thousand columns w/o difficulty
00032 int max_columns;
00033 int num_cols_to_delete; 
00034 // Global variable to control output
00035 bool verbose;
00036 time_t heur_time, mip_time;
00037 
00038 void OSI_recover_route(int id, int **orderings, VRPRoute *r)
00039 {
00044     
00045     // Now get the ordering using this id
00046     int i=0;
00047     r->num_customers=0;
00048     while(orderings[id][i]!=-1)
00049     {
00050         r->ordering[i]=orderings[id][i];
00051         r->num_customers++;
00052         i++;
00053     }
00054 
00055     // We now have the # of customers in the route and the ordering
00056     // Compute the two hash values 
00057     r->hash_val=r->hash(SALT_1);
00058     r->hash_val2=r->hash(SALT_2);
00059 
00060     return;
00061   
00062 }
00063 
00064 void OSI_recover_solution(OsiSolverInterface *si, int **orderings, int *sol_buff)
00065 {
00071  
00072     int i, j, k, id, nnodes, ncols;
00073     const double *x;
00074     std::string colname;
00075 
00076     nnodes= si->getNumRows();
00077     ncols= si->getNumCols();
00078     x = si->getColSolution();
00079 
00080     sol_buff[0]=nnodes;
00081 
00082     VRPRoute route(nnodes);
00083 
00084     k=1;
00085     for(i=0;i<ncols;i++)
00086     {
00087         if(VRPH_ABS(x[i]-1)<.01)
00088         {
00089             colname=si->getColName(i,VRP_INFINITY);
00090             id=atoi(colname.c_str());
00091 
00092             OSI_recover_route(id,orderings,&route);
00093         
00094             for(j=0;j<route.num_customers;j++)
00095                 sol_buff[k+j]=route.ordering[j];
00096 
00097             sol_buff[k]=-sol_buff[k];
00098 
00099             k += route.num_customers;
00100         }
00101     }
00102     sol_buff[nnodes+1]=0;
00103 
00104     return;
00105   
00106 }
00107 
00108 void OSI_add_route(OsiSolverInterface *si, VRP *V, VRPRoute *r, int id, int **orderings)
00109 {
00114 
00115     if(id>MAX_ROUTES)
00116     {   
00117         fprintf(stderr,"Too many routes added! Increase value of MAX_ROUTES=%d\n",MAX_ROUTES);
00118         exit(-1);
00119     }   
00120 
00121     int i, j, k, nvars, ncols;
00122     int *cols_to_delete;
00123     int *col_indices;
00124     std::string col_name;
00125     const double *U;
00126     time_t start, stop;
00127 
00128     // Create a new route object of the right size
00129     VRPRoute removed_route(V->get_num_nodes());
00130 
00131     // Create a CoinPackedVector for the new column
00132     CoinPackedVector col;
00133 
00134     for(i=0;i<r->num_customers;i++)
00135         col.insert(r->ordering[i]-1, 1);// 0-based
00136 
00137     // Check to see if we have exceeded max_columns
00138     // Since we have to fix variable values instead of deleting them in GLPK,
00139     // we have to calculate this
00140     U=si->getColUpper();
00141     ncols=si->getNumCols();
00142     nvars=0;
00143     for(i=0 ; i < ncols ; i++)
00144     {
00145         if(U[i]==1)
00146             nvars++;
00147     }
00148 
00149     if(nvars>=max_columns)
00150     {
00151         // Solve the problem and then delete num_cols_to_delete variables not in the current
00152         // optimal solution
00153         start=clock();
00154         si->branchAndBound();
00155         stop=clock();
00156         mip_time += (stop-start); 
00157         const double *y=si->getColSolution();
00158         // Select a set of num_cols_to_delete columns to delete   
00159         int nn=si->getNumCols();
00160         cols_to_delete = new int[nn];
00161         col_indices= new int[num_cols_to_delete];
00162         memset(cols_to_delete,0,nn*sizeof(int));
00163         memset(col_indices,0,num_cols_to_delete*sizeof(int));
00164         U=si->getColUpper();
00165         
00166         if(verbose)
00167             printf("Finding %d columns to delete \n",num_cols_to_delete);
00168         k=0;
00169         while(k<num_cols_to_delete)
00170         {
00171             while(true)
00172             {
00173                 j=(int)floor(lcgrand(10)*nn);
00174                 // Make sure the variable/column is not in the optimal solution 
00175                 // Since OSI/GLPK implementation of deleteCols doesn't seem right,
00176                 // Need to make sure that this column isn't already fixed!
00177                 if(y[j]<0.01 && cols_to_delete[j]==0 && U[j]==1)
00178                 {
00179                     cols_to_delete[j]=1;
00180                     col_indices[k]=j;
00181                     k++;
00182                     break;
00183                 }
00184             }    
00185         }
00186         
00187         // Delete the columns - this does not seem to work in GLPK - just fix the variable to 0
00188         // si->deleteCols(num_cols_to_delete,col_indices);
00189         
00190         for(k=0;k<num_cols_to_delete;k++)
00191         {
00192             col_name=si->getColName(col_indices[k],VRP_INFINITY);
00193             // Remove the route from VRPH's "Route Warehouse"
00194             int route_id=atoi(col_name.c_str());
00195             OSI_recover_route(route_id,orderings, &removed_route);
00196             
00197             V->route_wh->remove_route(removed_route.hash_val, removed_route.hash_val2);
00198             // Fix this variable/column to 0
00199             si->setColBounds(col_indices[k],0,0);
00200         }
00201         if(verbose)
00202             printf("Deleted/Fixed %d columns/variables\n",num_cols_to_delete);
00203     
00204         delete [] cols_to_delete;
00205         delete [] col_indices;
00206     }
00207 
00208     // Add the named column to the IP 
00209     sprintf(r->name,"%d",id);
00210     std::string new_name(r->name);
00211     si->addCol(col,0,1,r->length-r->total_service_time,new_name);
00212     ncols=si->getNumCols();
00213     // Set as integer (binary)
00214     si->setInteger(ncols-1); // 0-based
00215 
00216     // Copy the ordering of the route to the orderings[] array
00217     // Add a -1 to denote the end
00218     orderings[id]=new int[r->num_customers+1];
00219     memcpy(orderings[id],r->ordering,r->num_customers*sizeof(int));
00220     orderings[id][r->num_customers]=-1;
00221     return;
00222 }
00223 
00224 
00225 int main(int argc, char **argv)
00226 {
00227     VRPH_version();
00228 
00229     int i, j, k, n, status, num_attempts, *sol_buff, *IP_sol_buff;
00230     char in_file[200];
00231     double lambda, best_heur_sol=VRP_INFINITY;
00232     bool first_sol=false, bootstrap=false;;
00233     VRPSolution *fresh_solution;
00234     OsiSolverInterface *si;
00235     const double *x;
00236     int last_num_cols=0, route_id=0;
00237     time_t start, stop;
00238     int *orderings[MAX_ROUTES];
00239     for(i=0;i<MAX_ROUTES;i++)
00240         orderings[i]=NULL;
00241     
00242 
00243     // Set timing counters to 0
00244     heur_time=mip_time=0;
00245 
00246     // Check arguments
00247     if(argc<5)
00248     {
00249         fprintf(stderr,"Usage: %s -f input_file -n num_runs [-v,-b,-c max_columns -d cols_to_delete]\n",
00250             argv[0]);
00251         fprintf(stderr,"\t Will solve the problem num_solutions times and add the routes\n");
00252         fprintf(stderr,"\t to a set partitioning problem.\n");
00253         fprintf(stderr,"\t Other options:\n");
00254         fprintf(stderr,"\t -v runs in verbose mode\n");
00255         fprintf(stderr,"\t -b will use bootstrapping where we send the set partitioning\n"
00256                        "\t    solution back to the metaheuristic solver\n");
00257         fprintf(stderr,"\t -c max_columns will allow this many active columns/variables in the IP.\n");
00258         fprintf(stderr,"\t    Default value is max_columns=500\n");
00259         fprintf(stderr,"\t -d num_cols_to_delete will delete this many columns once we have too many\n");
00260         fprintf(stderr,"\t    in the IP. Default value is num_cols_to_delete=100\n");
00261         exit(-1);
00262     }
00263 
00264     // Set defaults
00265     verbose=false;
00266     max_columns=500;
00267     num_cols_to_delete=100;
00268 
00269     // Parse command line
00270     for(i=0;i<argc;i++)
00271     {
00272         if(strcmp(argv[i],"-f")==0)
00273             strcpy(in_file,argv[i+1]);
00274         if(strcmp(argv[i],"-n")==0)
00275             num_attempts=atoi(argv[i+1]);
00276         if(strcmp(argv[i],"-v")==0)
00277             verbose=true;
00278          if(strcmp(argv[i],"-b")==0)
00279             bootstrap=true;
00280          if(strcmp(argv[i],"-c")==0)
00281              max_columns=atoi(argv[i+1]);
00282          if(strcmp(argv[i],"-d")==0)
00283              num_cols_to_delete=atoi(argv[i+1]);
00284     }
00285 
00286     // This is the # of non-VRPH_DEPOT nodes
00287     n=VRPGetDimension(in_file);
00288     // This will be used to import/export solutions
00289     fresh_solution = new VRPSolution(n);
00290     // Create buffers for importing solutions
00291     sol_buff= new int[n+2];
00292     IP_sol_buff = new int[n+2];
00293 
00294     // Declare an OSI interface
00295     si=new OsiGlpkSolverInterface;
00296     si->setIntParam(OsiNameDiscipline,2);
00297 
00298     for(i=0;i<n;i++)
00299     {
00300         si->addRow(0,NULL,NULL,1,1);
00301     }
00302 
00303     // Declare a VRP of the right size and import the file
00304     VRP V(n);
00305     ClarkeWright CW(n);
00306     VRPRoute route(n);
00307     V.read_TSPLIB_file(in_file);
00308     // Set up a "route warehouse" to store the routes to be added to the IP
00309     V.route_wh=new VRPRouteWarehouse(HASH_TABLE_SIZE);
00310   
00311     // Set up a minimization problem
00312     si->setObjSense(1);
00313     // Set to error only output
00314     si->setHintParam(OsiDoReducePrint,true, OsiHintDo);
00315     // Unfortunately GLPK still prints out something regarding the conflict graph
00316 
00317     for(i=0;i<num_attempts;i++)
00318     {
00319         if(i==0 || !bootstrap) 
00320         {
00321             lambda=.5+1.5*lcgrand(0);
00322             // Start with a clean VRP object
00323             V.reset();
00324             CW.Construct(&V, lambda, false);
00325             if(verbose)
00326                 printf("CW solution %d[%5.3f]: %f\n",i,lambda,V.get_total_route_length()-V.get_total_service_time());
00327         }
00328         else
00329             // Use the solution from the IP
00330             V.import_solution_buff(IP_sol_buff);
00331         
00332         // Run VRPH's RTR algorithm to improve the solution
00333         start=clock();
00334         V.RTR_solve(ONE_POINT_MOVE | TWO_POINT_MOVE | TWO_OPT | VRPH_USE_NEIGHBOR_LIST,
00335             30, 5, 2, .01, 30, VRPH_LI_PERTURB, VRPH_FIRST_ACCEPT,false);
00336         stop=clock();
00337         heur_time += (stop-start);
00338 
00339         if(verbose)
00340             printf("RTR Metaheuristic found solution %5.3f\n",V.get_total_route_length()-V.get_total_service_time());
00341         
00342         // The RTR algorithm keeps a "warehouse" of the best solutions discovered during
00343         // the algorithm's search
00344         // Now go through the solutions in the solution warehouse and add the new routes
00345         // discovered to the IP
00346         for(j=0;j<V.solution_wh->num_sols;j++)
00347         {
00348             // Import solution j from the warehouse 
00349             V.import_solution_buff(V.solution_wh->sols[j].sol);
00350                 
00351             if(V.get_total_route_length()-V.get_total_service_time() < best_heur_sol)
00352                 best_heur_sol = V.get_total_route_length()-V.get_total_service_time() ;
00353 
00354             // Now add the routes from this solution to the IP
00355             for(k=1;k<=V.get_total_number_of_routes();k++)
00356             {
00357                 // Clean up the route by running INTRA_ROUTE optimizations only
00358                 // using the route_search method of the different local search
00359                 // heuristics, accepting improving moves only (VRPH_DOWNHILL)
00360                 OnePointMove OPM;
00361                 TwoOpt TO;
00362                 ThreeOpt ThO;
00363                 while(OPM.route_search(&V,k,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY )){}
00364                 while(TO.route_search(&V,k,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY  )){};
00365                 while(ThO.route_search(&V,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY )){};
00366 
00367                 // Copy route k from the solution to the VRPRoute R
00368                 V.update_route(k,&route);
00369                 route.create_name();
00370                 // Add it to the "route warehouse" - this uses a hash table to keep track
00371                 // of duplicate columns
00372                 status=V.route_wh->add_route(&route);
00373 
00374                 if(status!=DUPLICATE_ROUTE)
00375                 {
00376                     // This route is not currently in the WH and so it cannot be in the
00377                     // set partitioning problem
00378                     //OSI_add_route(si,&V,&route);
00379                     OSI_add_route(si,&V,&route,route_id,orderings);
00380                     route_id++;
00381                 }
00382             }
00383             
00384             // Set the row RHS's if we need to
00385             if(first_sol)
00386             {
00387                 first_sol=false;
00388                 for(int rownum=0;rownum<n;rownum++)
00389                     si->setRowBounds(rownum,1,1);
00390                 // Note that changing this to >= would be a set covering problem
00391                 // where each customer can be visited by more than one route
00392             }
00393         }
00394 
00395         // Now erase all the solutions from the WH
00396         V.solution_wh->liquidate();
00397 
00398         if(verbose)
00399         {
00400             printf("Attempt %02d:  Solving IP with %d columns\n",i,si->getNumCols());
00401             printf("%d routes in the WH\n",V.route_wh->num_unique_routes);
00402         }
00403 
00404         // Solve the current set partitioning problem using the MIP solver
00405         start=clock();
00406         si->branchAndBound();
00407         stop=clock();
00408         mip_time += (stop-start);
00409 
00410         double opt=si->getObjValue();
00411         x=si->getColSolution();
00412         last_num_cols=si->getNumCols();
00413         if(verbose)
00414             printf("Optimal solution (%d columns) is %f\n",last_num_cols,opt);
00415 
00416         
00417         // Now recover the solution from the IP solution
00418         OSI_recover_solution(si, orderings, IP_sol_buff);
00419        
00420         if(verbose)
00421             printf("IP solution has obj. function value: %5.2f\n"
00422             "Best heuristic  obj. function value: %5.2f\n",
00423             si->getObjValue(),best_heur_sol);     
00424     }
00425 
00426     if(verbose)
00427         printf(
00428         "\nResults\n"
00429         "--------\n"
00430         "After %d runs\n"
00431         "IP solution has obj. function value: %5.2f\n"
00432         "Best heuristic  obj. function value: %5.2f\n",
00433         num_attempts,si->getObjValue(),best_heur_sol);
00434 
00435     // print to stderr since GLPK prints "conflict graph" to stdout (a known bug...)
00436     // best_heur_sol best_mip_sol heur_time mip_time
00437     fprintf(stderr,"%5.3f %5.3f %5.3f %5.3f\n", best_heur_sol, si->getObjValue(),
00438         (double)(heur_time)/CLOCKS_PER_SEC, (double)(mip_time)/CLOCKS_PER_SEC);
00439 
00440     delete V.route_wh;
00441     delete fresh_solution;
00442     delete [] sol_buff;
00443     delete [] IP_sol_buff;
00444     delete si;
00445 
00446     for(i=0;i<MAX_ROUTES;i++)
00447         if(orderings[i])
00448             delete [] orderings[i];
00449 
00450      
00451     return 0;
00452 }
00453 
00454 
00455