Rob. I want exactly this. I tried the approach. It is giving me errors. See below. I have the code as well after this one. I think i am not setting the model to read it correctly.
256 endsource;
257 action optimization.runOptmodel / code=pgm groupBy='ModelWeekNumber' statusOut='statusOut';
258 quit;
NOTE: Active Session now MYSESS.
ModelWeekNumber=3
NOTE: There were 1 rows read from table 'VRPDATA' in caslib 'CASUSER(SRHY-Exp)'.
NOTE: Problem generation will use 16 threads.
ERROR: The array subscript 'UseNode[99999,1]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,2]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,3]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,4]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,5]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,6]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,7]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,8]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,9]' is invalid at line 51 column 23.
ERROR: The array subscript 'UseNode[99999,10]' is invalid at line 51 column 23.
NOTE: The maximum error limit was reached during OPTMODEL problem generation.
/*Variables*/
%let capacity = 3400;
%let num_vehicles=20;
%let dist=3000;
%let max_stops=5;
%let Depot='A';
%let DepotID=99999; /* Bfor Lebanon and Afor Lebec*/
/*%Let WeekNumber='2'; */
%Let Mode='Mode1'; /* Mode1 or Mode2*/
%let Scenario='Routing Alternates';
proc fedsql sessref=mysess;
create table CASUSER.vrpdata {options replace='TRUE'} AS
select ModelWeekNumber_2WK as ModelWeekNumber
,ModelNode
,ModelCity
,ModelLat
,ModelLong
,ModelZip
,ModelState
,Case_Count
,Unit_Count
,round(Volume)as Volume
,Weight
,Declared_Value
,ModelDepotState
,ModelDepotCity
,ModelDepotZip
,NodeID
,Mode
,WeekNumbers
,WeekNumber_EarliestStartDate
,WeekNumber_LatestEndDate
from Data
where ModelDepotCity=&Depot
and ModelWeekNumber_2WK in ('3','4')
and Mode=&Mode;
;
QUIT;
proc cas;
source pgm;
put _BY_LINE_;
/* read the node location and demand data */
set NODES;
str ModelNode {NODES};
num ModelLat {NODES};
num ModelLong {NODES};
num Volume {NODES};
num capacity = &capacity;
num num_vehicles = &num_vehicles;
str Mode {NODES};
read data CASUSER.vrpdata into NODES =[NodeID] ModelNode ModelLat ModelLong Volume Mode;
set ARCS = {i in NODES, j in NODES: i ne j};
set VEHICLES = 1..num_vehicles;
num depot = &DepotID;
/* define the arc dist as the rounded Euclidean distance */
num dist {<i,j> in ARCS} = (if i ne depot then GEODIST(ModelLat[i],ModelLong[i],ModelLat[j],ModelLong[j],'M') else 0);
/*print dist;*/
/* Creating an Sparse Matrix to remove any unwanted arcs*/
set LOGICAL_ARCS = {<i,j> in ARCS: dist[i,j] <= &dist};
/*print {<i,j> in LOGICAL_ARCS} dist; */
/* Flow[k] is the amount of demand carried 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, k in VEHICLES} dist[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; */
con VisitDepot {i in NODES diff {depot}, k in VEHICLES}:
UseNode[i,k] <= UseNode[depot,k]
suffixes=(block=k);
/*DISTANCE CONSTRAINT*/
CON Distance {k in VEHICLES}:
/* SUM {<i,j> in LOGICAL_ARCS: j not in {depot}} (dist[i, j] * UseArc[i, j, k]) <= &dist_limit; */
sum {<i,j> in LOGICAL_ARCS} dist[i,j] * UseArc[i,j,k] <= &dist;
/*STOPS CONSTRAINT*/
CON Stops {k in VEHICLES}:
sum {i in NODES diff {depot}} UseNode[i,k] <= &max_stops;
/* 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];
/* 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 Volume supplied by vehicle k to node i must equal Volume if UseNode[i,k] = 1; otherwise, it must equal 0 */
/*FLOW BALANCE CONSTRAINT*/
con FlowBalance {i in NODES diff {depot}, k in VEHICLES}:
sum {<(i),j> in LOGICAL_ARCS} Flow[i,j,k] - sum {<j,(i)> in LOGICAL_ARCS} Flow[j,i,k]
= Volume[i] * UseNode[i,k];
/* each vehicle is empty when it leaves the depot */
for {<(depot),j> in LOGICAL_ARCS, k in VEHICLES} fix Flow[depot,j,k] = 0;
/* 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];
/*Constraint VehicleCapacity[3,4,1]: Flow[3,4,1] - 3650*UseArc[3,4,1] <= 0 /*
/* 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 {k in VEHICLES} Distance[k].block = k;
for {k in VEHICLES} Stops[k].block = k;
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 / decomp varsel=ryanfoster presolver=basic maxtime=1200;
/* OUTPUT */
/* create solution data set */
create data solution_data from [i j k]=
{<i,j> in LOGICAL_ARCS, k in VEHICLES: UseArc[i,j,k].sol > 0.5 and i ne depot}
StartNode=ModelNode[i] StartCity=ModelNode[i] x1=ModelLat[i] y1=ModelLong[i]
EndNode=ModelNode[j] EndCity=ModelNode[j] x2=ModelLat[j] y2=ModelLong[j]
flow=Flow[i,j,k] dist=dist[i,j]
;
endsource;
action optimization.runOptmodel / code=pgm groupBy='ModelWeekNumber' statusOut='statusOut';
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;
DATA WORK.VRO_OUTPUT;
set WORK.solution_data;
Depot=&Depot;
WeekNumber=&WeekNumber;
VehicleByScenarioByWeekNumber= cat(Depot,'-',k,'-',WeekNumber);
Mode=&Mode;
Scenario=&Scenario;
run;
