Accepted Solutions
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;
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;Hi Rob
I wanted to try the solve with optwork ,the special vrp solver that you mentioned. However I get an error that it is not able to recognize the option VRP. Here is the log file and the SAS Viya release detail below. Can you please help in terms of what I need to do about this issue Rob?
Thanks much
%let cost_limit = 963;
173
174 proc optmodel sessref=mysess;
175 /* read the node location and demand data */
176 set NODES;
177 str nodename {NODES};
178 str city {NODES};
179 num x {NODES};
180 num y {NODES};
181 num demand {NODES};
182 read data vrpdata into NODES=[node] nodename city x y demand;
NOTE: There were 42 observations read from the data set WORK.VRPDATA.
183 set ARCS = {i in NODES, j in NODES: i ne j};
184 set DEPOT = {1};
185 num cost {<i,j> in ARCS} = (if j not in DEPOT then GEODIST(x[i],y[i],x[j],y[j],'M') else 0);
186
187 set <num,num,num,num> VRPLINKS;
187! /* route, order, from, to */
188 num twl{NODES} = 0;
189 num twu{NODES} = &cost_limit;
190 num route{NODES};
191 num route_order{NODES};
192 num delivery_time{NODES};
193
194 /* call network solver for vehicle routing with time windows */
195 solve with network /
196 direction = directed
197 links = (weight=cost)
198 vrp = (depot=DEPOT
---
690
22
ERROR 690-782: The option VRP is unrecognized.
ERROR 22-322: Expecting one of the following: BICONCOMP, BICONNECTEDCOMPONENTS, CLIQUE, CONCOMP, CONNECTEDCOMPONENTS, CYCLE,
GRAPH_DIRECTION, INCLUDE_SELFLINK, LAP, LINEARASSIGNMENT, LINKS, LOGFREQ, LOGLEVEL, MAXTIME, MCF, MINCOSTFLOW,
MINCUT, MINSPANTREE, MST, NODES, NTHREADS, OUT, SHORTESTPATH, SHORTPATH, SUBGRAPH, TIMETYPE, TRANSCL,
TRANSITIVECLOSURE, TSP.
199 capacity=&capacity
--------
22
76
ERROR 22-322: Syntax error, expecting one of the following: !, !!, &, (, ), *, **, +, -, .., /, <, <=, <>, =, >, ><, >=, AND, BY,
CROSS, DIFF, ELSE, IN, INTER, OR, SYMDIFF, TO, UNION, WITHIN, [, ^, ^=, |, ||, ~=.
ERROR 76-322: Syntax error, statement will be ignored.
200 demand=demand
201 maxRoutes=&num_vehicles
202 timeWindowLower=twl
203 timeWindowUpper=twu
204 milp=FALSE)
205 out = (route=route
206 order=route_order
207 arrival=delivery_time
208 vrplinks=VRPLINKS);
209
210 set ROUTES init {};
211 num costUsed {ROUTES} init 0;
212 for {<r,o,i,j> in VRPLINKS} do;
213 ROUTES = ROUTES union {r};
214 costUsed[r] = costUsed[r] + cost[i,j];
215 end;
ERROR: The symbol 'VRPLINKS' has no value at line 212 column 22.
216 print costUsed;
217
218 /* create solution data sets */
219 create data RoutesNodes
220 from [node] demand twl twu route route_order delivery_time;
NOTE: The data set WORK.ROUTESNODES has 42 observations and 7 variables.
221 create data RoutesLinks(where=(cost ne 0))
222 from [r order i j]=VRPLINKS cost[i,j]
223 StartNode=nodename[i] StartCity=city[i] x1=x[i] y1=y[i]
224 EndNode=nodename[j] EndCity=city[j] x2=x[j] y2=y[j]
225 function='line' drawspace='datavalue';
ERROR: The symbol 'VRPLINKS' has no value at line 222 column 26.
226 quit;
NOTE: The SAS System stopped processing this step because of errors.
NOTE: PROCEDURE OPTMODEL used (Total process time):
real time 0.01 seconds
cpu time 0.01 seconds
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.
-
8858 views
-
3 likes
-
-
