VRPH
1.0
|
00001 00002 // // 00003 // This file is part of the VRPH software package for // 00004 // generating solutions to vehicle routing problems. // 00005 // VRPH was developed by Chris Groer (cgroer@gmail.com). // 00006 // // 00007 // (c) Copyright 2010 Chris Groer. // 00008 // All Rights Reserved. VRPH is licensed under the // 00009 // Common Public License. See LICENSE file for details. // 00010 // // 00012 00013 // 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