Hi Rob
this is directly related to the https://communities.sas.com/t5/Mathematical-Optimization/VRO-Question/m-p/879258#M3914
I do not have SAS Viya 4 so don't have the luxury of a special vrp solver. So settled with what I have but I am not able to get the model to solve. I have the code below. The model takes forever.. I stopped and I do not know where/why the model takes a long time. There are no nodes that are more than the fixed number (600 miles). That has been vetted out. Trying to have some sustainable and robust model so that I can run multiple VROs with probably 250 nodes. The demand range is single digits and goes as big as 34,000 ish. Is this something that is to be looked at? Not sure but am just thinking out loud.This model is not solvable, at least in my machine in my environment. Any help / leads is appreciated Rob.
/* number of vehicles available */
%let num_vehicles = 10;
/* capacity of each vehicle */
%let capacity = 40000;
%let cost_limit = 600;
/* node, x coordinate, y coordinate, demand */
data vrpdata;
input node nodename $ city $ x y demand;
datalines;
1 HAR HARRISBURG 39.7757 -86.1092 0
2 CLE C-44146 41.3889 -81.5272 6826
3 CLE C-44266 41.1677 -81.1638 109
4 CLE C-44707 40.7628 -81.3463 106
5 CMH C-25313 38.4143 -81.7652 5257
6 CMH C-26150 39.1584 -81.543 693
7 CMH C-43026 40.024 -83.1801 2942
8 CMH C-43078 40.1186 -83.7536 5510
9 CMH C-45601 39.3093 -82.9736 968
10 CVG C-41701 37.2923 -83.0836 7662
11 CVG C-45069 39.3377 -84.4056 9023
13 DTW C-43537 41.5745 -83.6734 4707
14 DTW C-48036 42.5987 -82.908 3049
15 DTW C-48120 42.3048 -83.1947 4119
16 DTW C-48165 42.4961 -83.6172 5440
17 DTW C-48170 42.3637 -83.5369 700
18 DTW C-48210 42.3377 -83.1278 838
19 DTW C-48601 43.4138 -83.8508 9565
20 DTW C-49508 42.8666 -85.6179 8559
21 IND C-46032 39.9727 -86.1671 4747
22 IND C-46052 40.0438 -86.4828 315
23 IND C-46075 40.0285 -86.3318 5830
24 IND C-46168 39.6824 -86.3897 34414
25 IND C-46226 39.8402 -86.0565 120
27 IND C-46242 39.7539 -86.1908 4385
28 IND C-46268 39.9009 -86.2329 21456
29 IND C-46733 40.8299 -84.9372 5
30 IND C-46808 41.1007 -85.1766 10296
32 IND C-47202 39.2023 -85.922 10
33 IND C-47203 39.2312 -85.8174 34
34 IND C-47265 39.0156 -85.6402 31
35 IND C-47274 38.9579 -85.9085 17085
41 MKE C-53005 43.0607 -88.0965 6213
42 MKE C-53154 42.8862 -87.8972 8839
43 MKE C-53190 42.8056 -88.7269 5
44 MKE C-53578 43.3146 -89.7764 99
45 MKE C-53589 42.9323 -89.2011 205
46 MKE C-53707 43.0866 -89.3597 642
proc optmodel;
/* read the node location and demand data */
set NODES;
str nodename {NODES};
str city {NODES};
num x {NODES};
num y {NODES};
num demand {NODES};
num capacity = &capacity;
num num_vehicles = &num_vehicles;
read data vrpdata into NODES=[node] nodename city x y demand;
set ARCS = {i in NODES, j in NODES: i ne j};
set VEHICLES = 1..num_vehicles;
/* define the RDC1 as node 1 */ /*how to add more than one depot?? */
num depot = 1;
/* set DEPOT= {1}; */
/* define the arc cost as the rounded Euclidean distance */
num cost {<i,j> in ARCS} = (if j ne depot then GEODIST(x[i],y[i],x[j],y[j],'M') else 0);
print cost;
set LOGICAL_ARCS = {<i,j> in ARCS: cost[i,j] <= &cost_limit};
print {<i,j> in LOGICAL_ARCS} cost;
/* Flow[i,j,k] is the amount of demand carried on arc (i,j) by vehicle k */
var Flow {LOGICAL_ARCS, VEHICLES} >= 0 <= capacity;
/* UseNode[i,k] = 1, if and only if node i is serviced by vehicle k */
var UseNode {NODES, VEHICLES} binary;
/* UseArc[i,j,k] = 1, if and only if arc (i,j) is traversed by vehicle k */
var UseArc {LOGICAL_ARCS, VEHICLES} binary;
/*-----------------------OBJECTIVE FUNCTION-------------------------------*/
/* minimize the total distance traversed */
min TotalCost = sum {<i,j> in LOGICAL_ARCS : j ne depot} sum{k in VEHICLES} cost[i,j] * UseArc[i,j,k];
/*----------------------------CONSTRAINTS-------------------------------------*/
/* each non-depot node must be serviced by at least one vehicle */
/*ASSIGNMENT CONSTRAINT*/
con Assignment {i in NODES diff {depot}}:
sum {k in VEHICLES} UseNode[i,k] >= 1;
for{k in VEHICLES} fix UseNode[depot,k] = 1;
/* each vehicle must start at the depot node */
/*DISTANCE CONSTRAINT*/
CON Distance {k in VEHICLES}:
sum {<i,j> in LOGICAL_ARCS} cost[i,j] * UseArc[i,j,k] <= &cost_limit;
/* expand Distance; */
/* some vehicle k traverses an arc that leaves node i if and only if UseNode[i,k] = 1 */
/* LEAVE NODE CONSTRAINT */
con LeaveNode {i in NODES, k in VEHICLES}:
sum {<(i),j> in LOGICAL_ARCS} UseArc[i,j,k] = UseNode[i,k];
/* expand LeaveNode; */
/* some vehicle k traverses an arc that enters node i if and only if UseNode[i,k] = 1 */
/* ENTER NODE CONSTRAINT */
con EnterNode {i in NODES, k in VEHICLES}:
sum {<j,(i)> in LOGICAL_ARCS} UseArc[j,i,k] = UseNode[i,k];
/* the amount of demand supplied by vehicle k to node i must equal demand if UseNode[i,k] = 1; otherwise, it must equal 0 */
/*FLOW BALANCE CONSTRAINT*/
con FlowBalance {i in NODES diff {depot}, k in VEHICLES}:
sum {<j,(i)> in LOGICAL_ARCS} Flow[j,i,k] - sum {<(i),j> in LOGICAL_ARCS} Flow[i,j,k]
= demand[i] * UseNode[i,k];
/* expand FlowBalance; */
/* if UseArc[i,j,k] = 1, then the flow on arc (i,j) must be at most capacity
if UseArc[i,j,k] = 0, then no flow is allowed on arc (i,j) */
/*VEHICLE CAPACITY CONSTRAINT*/
con VehicleCapacity {<i,j> in LOGICAL_ARCS, k in VEHICLES}:
Flow[i,j,k] <= Flow[i,j,k].ub * UseArc[i,j,k];
/* expand; */
/* decomp by vehicle *//*DECOMPOSITION*/
for {i in NODES, k in VEHICLES} do;
LeaveNode[i,k].block = k;
EnterNode[i,k].block = k;
end;
for {i in NODES diff {depot}, k in VEHICLES} FlowBalance[i,k].block = k;
for {<i,j> in LOGICAL_ARCS, k in VEHICLES} VehicleCapacity[i,j,k].block = k;
/* solve using decomp (aggregate formulation) */
solve with MILP / varsel=ryanfoster decomp=(logfreq=20 hybrid=TRUE relobjgap=0.5);
/* OUTPUT */
/* create solution data set */
create data solution_data (where=(x1<>x2 and y1<>y2)) from [i j k]=
{<i,j> in LOGICAL_ARCS, k in VEHICLES: UseArc[i,j,k].sol > 0.5}
StartNode=nodename[i] StartCity=city[i] x1=x[i] y1=y[i]
EndNode=nodename[j] EndCity=city[j] x2=x[j] y2=y[j]
flow=Flow[i,j,k] cost=cost[i,j]
function='line' drawspace='datavalue';
quit;
proc sgplot data=solution_data noautolegend;
scatter x=x1 y=y1 / datalabel=i;
vector x=x2 y=y2 / xorigin=x1 yorigin=y1 group=k noarrowheads;
xaxis display=none;
yaxis display=none;
run;
