First, I have a few corrections and suggestions for your code.
In PROC OPTMODEL, <> is the MAX operator, but here you want to express "not equal":
/* set ARCS = {i in NODES, j in NODES: i <> j};*/
set ARCS = {i in NODES, j in NODES: i ne j};
You want to avoid counting the (dummy) trip back to the depot, and you accounted for that in the Distance constraint, but you also need to do so in the objective:
/* min TotalCost = sum {<i,j> in LOGICAL_ARCS, k in VEHICLES} cost[i,j] * UseArc[i,j,k];*/
min TotalCost = sum {<i,j> in LOGICAL_ARCS: j ne depot} sum {k in VEHICLES} cost[i,j] * UseArc[i,j,k];
An alternative approach is to modify the cost definition:
* num cost {<i,j> in ARCS} = GEODIST(x[i],y[i],x[j],y[j],'M'); /*round(sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2))*/
num cost {<i,j> in ARCS} = (if j ne depot then GEODIST(x[i],y[i],x[j],y[j],'M') else 0);
And then you can just omit the logical conditions for the sums in both the objective and the Distance constraint. This idea is also illustrated in the Lost Baggage Distribution example from the mathematical programming examples book.
The Distance constraint can also be associated with the block for vehicle k, and you should not exclude the depot in the blocks for LeaveNode and EnterNode:
for {k in VEHICLES} Distance[k].block = k;
for {i in NODES /*diff {depot}*/, 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;
A best practice is to introduce a macro variable instead of hardcoding 600 in multiple places:
%let cost_limit = 600;
...
/* set LOGICAL_ARCS = {i in NODES, j in NODES: cost[i,j] <= 600};*/
set LOGICAL_ARCS = {<i,j> in ARCS: cost[i,j] <= &cost_limit};
...
con Distance {k in VEHICLES}:
/* sum {<i,j> in LOGICAL_ARCS: j not in {depot}} (Cost[i, j] * UseArc[i, j, k]) <= 600;*/
sum {<i,j> in LOGICAL_ARCS} cost[i,j] * UseArc[i,j,k] <= &cost_limit;
But even after all these changes, you will not be able to find a solution because the problem is infeasible. Notice that the distance from the depot to node 15 is 962.99269, so a cost_limit of 600 cannot be satisfied. After a few minutes, the decomposition algorithm reports this infeasibility:
7704 solve with MILP / varsel=ryanfoster decomp presolver=moderate;
NOTE: Problem generation will use 4 threads.
NOTE: The problem has 33300 variables (0 free, 10 fixed).
NOTE: The problem has 16860 binary and 0 integer variables.
NOTE: The problem has 17741 linear constraints (16450 LE, 1250 EQ, 41 GE, 0 range).
NOTE: The problem has 115710 linear constraint coefficients.
NOTE: The problem has 0 nonlinear constraints (0 LE, 0 EQ, 0 GE, 0 range).
NOTE: The initial MILP heuristics are applied.
NOTE: The MILP presolver value MODERATE is applied.
NOTE: The MILP presolver removed 7450 variables and 3720 constraints.
NOTE: The MILP presolver removed 26060 constraint coefficients.
NOTE: The MILP presolver modified 0 constraint coefficients.
NOTE: The presolved problem has 25850 variables, 14021 constraints, and 89650 constraint
coefficients.
NOTE: The MILP solver is called.
NOTE: The Decomposition algorithm is used.
NOTE: The Decomposition algorithm is executing in single-machine mode.
NOTE: The DECOMP method value USER is applied.
NOTE: All blocks are identical and the master model is set covering.
WARNING: The master model is not a set partitioning and VARSEL=RYANFOSTER. The objective function
must ensure that there exists at least one optimal solution that fulfills all of the
master constraints at equality.
NOTE: The Decomposition algorithm is using an aggregate formulation and Ryan-Foster branching.
NOTE: The number of block threads has been reduced to 1 threads.
NOTE: The problem has a decomposable structure with 10 blocks. The largest block covers 9.971% of
the constraints in the problem.
NOTE: The decomposition subproblems cover 25850 (100%) variables and 13980 (99.71%) constraints.
NOTE: The deterministic parallel mode is enabled.
NOTE: The Decomposition algorithm is using up to 4 threads.
Iter Best Master Best LP IP CPU Real
Bound Objective Integer Gap Gap Time Time
NOTE: Starting phase 1.
1 0.0000 40.0000 . 4.00e+01 . 35 25
8 20.0000 20.0000 . 0.00% . 575 205
9 20.0000 20.0000 . 0.00% . 815 277
NOTE: Phase 1 continuous lower bound has positive value. The original problem is Infeasible.
Node Active Sols Best Best Gap CPU Real
Integer Bound Time Time
0 1 0 . 2328.4549 . 815 277
NOTE: The Decomposition algorithm used 4 threads.
NOTE: The Decomposition algorithm time is 277.68 seconds.
NOTE: Infeasible.
If you change cost_limit to 963, the problem turns out to be feasible. But the decomposition algorithm documentation example is not the best way to solve a vehicle routing problem, and a note at the top of that example recommends instead using a specialized VRP solver that is available in SAS Optimization in SAS Viya. Support was recently added for time windows, and imposing a time window upper bound is a simple way to model the cost_limit. The following code solves your problem instantly:
%let cost_limit = 963;
proc optmodel sessref=mysess;
/* read the node location and demand data */
set NODES;
str nodename {NODES};
str city {NODES};
num x {NODES};
num y {NODES};
num demand {NODES};
read data vrpdata into NODES=[node] nodename city x y demand;
set ARCS = {i in NODES, j in NODES: i ne j};
set DEPOT = {1};
num cost {<i,j> in ARCS} = (if j not in DEPOT then GEODIST(x[i],y[i],x[j],y[j],'M') else 0);
set <num,num,num,num> VRPLINKS; /* route, order, from, to */
num twl{NODES} = 0;
num twu{NODES} = &cost_limit;
num route{NODES};
num route_order{NODES};
num delivery_time{NODES};
/* call network solver for vehicle routing with time windows */
solve with network /
direction = directed
links = (weight=cost)
vrp = (depot=DEPOT
capacity=&capacity
demand=demand
maxRoutes=&num_vehicles
timeWindowLower=twl
timeWindowUpper=twu
milp=FALSE)
out = (route=route
order=route_order
arrival=delivery_time
vrplinks=VRPLINKS);
set ROUTES init {};
num costUsed {ROUTES} init 0;
for {<r,o,i,j> in VRPLINKS} do;
ROUTES = ROUTES union {r};
costUsed[r] = costUsed[r] + cost[i,j];
end;
print costUsed;
/* create solution data sets */
create data RoutesNodes
from [node] demand twl twu route route_order delivery_time;
create data RoutesLinks(where=(cost ne 0))
from [r order i j]=VRPLINKS cost[i,j]
StartNode=nodename[i] StartCity=city[i] x1=x[i] y1=y[i]
EndNode=nodename[j] EndCity=city[j] x2=x[j] y2=y[j]
function='line' drawspace='datavalue';
quit;
proc sgplot data=RoutesLinks;
scatter x=x2 y=y2 / datalabel=j;
vector x=x2 y=y2 / xorigin=x1 yorigin=y1 group=r noarrowheads;
xaxis display=none;
yaxis display=none;
run;
... View more
