Although you could find the connected components of the support graph of each intermediate solution by using OPTMODEL’s programming language capabilities as in Chapter 23, the following SAS macro instead uses the SOLVE WITH NETWORK statement together with the CONCOMP option:

%macro findConnectedComponents;
   if card(ARCS_SOL) > 0 then do;
      solve with NETWORK /
         links    = (include=ARCS_SOL)
         subgraph = (nodes=NODES_SOL)
         concomp
         out      = (concomp=component_id);
      COMPONENT_IDS = setof {i in NODES_SOL} component_id[i];
      for {c in COMPONENT_IDS} COMPONENT[c] = {};
      for {i in NODES_SOL} do;
         ci = component_id[i];
         COMPONENT[ci] = COMPONENT[ci] union {i};
      end;
   end;
   else COMPONENT_IDS = {};
%mend findConnectedComponents;

The following SAS macro contains a DO UNTIL loop that implements dynamic generation of subtour elimination constraints (“row generation”), as in Chapter 23:

%macro subtourEliminationLoop;
   /* loop until each vehicle's support graph is connected */
   do until (and {v in VEHICLES} num_components[v] <= 1);
      solve;
      /* find connected components for each vehicle */
      for {v in VEHICLES} do;
         NODES_SOL = {i in NODES: UseNode[i,v].sol > 0.5};
         ARCS_SOL = {<i,j> in ARCS: UseArc[i,j,v].sol > 0.5};
         %findConnectedComponents;
         num_components[v] = card(COMPONENT_IDS);
         /* create subtour from each component not containing depot node */
         for {k in COMPONENT_IDS: depot not in COMPONENT[k]} do;
            num_subtours = num_subtours + 1;
            SUBTOUR[num_subtours] = COMPONENT[k];
            put SUBTOUR[num_subtours]=;
         end;
      end;
      print UseVehicle TimeUsed num_components;
   end;
%mend subtourEliminationLoop;

The previous two macros are used within the main PROC OPTMODEL call. The first several PROC OPTMODEL statements declare index sets and parameters and then read the input data:

proc optmodel;
   num num_vehicles init &num_vehicles;
   set VEHICLES = 1..num_vehicles;
   str depot = "&depot";
   set <str> NODES;
   read data time_data into NODES=[location];
   set ARCS init NODES cross NODES;
   num travel_time {ARCS};
   read data time_data into [i=location]
      {j in NODES} <travel_time[i,j]=col(j)>;

The following statements make the travel times symmetric, except that travel time back to the depot is set to 0:

   for {<i,j> in ARCS: travel_time[i,j] = .}
      travel_time[i,j] = travel_time[j,i];
   /* ignore travel time back to depot */
   for {i in NODES}
      travel_time[i,depot] = 0;
   /* remove self-loops */
   ARCS = ARCS diff setof {i in NODES} <i,i>;
   print travel_time;

The PRINT statement results in the first section of output, shown in Figure 27.1.

The following model declaration statements correspond directly to the mathematical programming formulation that is described earlier:

   var UseNode {NODES, VEHICLES} binary;
   var UseArc {ARCS, VEHICLES} binary;
   var UseVehicle {VEHICLES} binary;
   var TimeUsed {v in VEHICLES} >= 0 <= &time_limit;

   min NumVehiclesUsed =
      sum {v in VEHICLES} UseVehicle[v];

   con NodeCover {i in NODES diff {depot}}:
      sum {v in VEHICLES} UseNode[i,v] = 1;

   con Outflow {i in NODES, v in VEHICLES}:
      sum {<(i),j> in ARCS} UseArc[i,j,v] = UseNode[i,v];

   con Inflow {j in NODES, v in VEHICLES}:
      sum {<i,(j)> in ARCS} UseArc[i,j,v] = UseNode[j,v];

   con UseVehicle_con1 {i in NODES, v in VEHICLES}:
      UseNode[i,v] <= UseVehicle[v];

   con UseVehicle_con2 {v in VEHICLES}:
      UseVehicle[v] <= UseNode[depot,v];

   con TimeUsed_con {v in VEHICLES}:
      TimeUsed[v] = sum {<i,j> in ARCS} travel_time[i,j] * UseArc[i,j,v];

The following statements declare optional symmetry-breaking constraints to reduce the number of essentially identical branch-and-bound nodes that are explored by the mixed integer linear programming solver:

   /* several alternatives for symmetry-breaking constraints */
   con Symmetry {v in VEHICLES diff {1}}:
      sum {i in NODES} UseNode[i,v] <= sum {i in NODES} UseNode[i,v-1];
*  con Symmetry {v in VEHICLES diff {1}}:
      UseVehicle[v] <= UseVehicle[v-1];
*  con Symmetry {v in VEHICLES diff {1}}:
      TimeUsed[v] <= TimeUsed[v-1];

Williams (2013) breaks symmetry by using the first alternative, but you could use one of the other alternatives instead.

In SAS/OR 13.1, the mixed integer linear programming solver automatically detects and exploits symmetry without you having to explicitly declare such symmetry-breaking constraints. You can control the aggressiveness of symmetry detection by using the SYMMETRY= option in the SOLVE WITH MILP statement.

The following statements declare the subtour elimination constraints:

   num num_subtours init 0;

   /* subset of nodes not containing depot node */
   set <str> SUBTOUR {1..num_subtours};

   /* if node k in SUBTOUR[s] is used by vehicle v, then
      must use at least two arcs across partition induced by SUBTOUR[s] */
   con Subtour_elimination
      {s in 1..num_subtours, k in SUBTOUR[s], v in VEHICLES}:
      sum {i in NODES diff SUBTOUR[s], j in SUBTOUR[s]: <i,j> in ARCS}
         UseArc[i,j,v]
    + sum {i in SUBTOUR[s], j in NODES diff SUBTOUR[s]: <i,j> in ARCS}
         UseArc[i,j,v]
   >= 2 * UseNode[k,v];

The following statements declare the index sets and parameters that are needed to detect violated subtour elimination constraints:

   num num_components {VEHICLES};
   set <str> NODES_SOL;
   set <str,str> ARCS_SOL;
   num component_id {NODES_SOL};
   set COMPONENT_IDS;
   set <str> COMPONENT {COMPONENT_IDS};
   num ci;

The following statements call the %subtourEliminationLoop macro to minimize the number of vehicles used and then use the .sol objective suffix to update the num_vehicles parameter to the resulting minimum value:

   %subtourEliminationLoop;
   num_vehicles = round(NumVehiclesUsed.sol);

Changing the value of num_vehicles automatically updates the VEHICLES index set and consequently all the model declarations that depend on VEHICLES.

The following statements declare the additional variables, objective, and constraints that are needed to minimize the makespan, given the minimum number of vehicles already found:

   var MaxTimeUsed >= 0 <= &time_limit;

   min Makespan = MaxTimeUsed;

   con MaxTimeUsed_con {v in VEHICLES}:
      MaxTimeUsed >= TimeUsed[v];

The following statements call the %subtourEliminationLoop macro again to minimize the makespan and then print the nodes and arcs that are used by each vehicle in the final solution:

   %subtourEliminationLoop;
   for {v in VEHICLES: UseVehicle[v].sol > 0.5} do;
      print v;
      print {<i,j> in ARCS: UseArc[i,j,v].sol > 0.5} travel_time[i,j];
   end;
quit;

Figure 27.2 through Figure 27.7 show the output from each iteration of subtour elimination.

Figure 27.8 shows the final problem and solution summaries from the mixed integer linear programming solver.

Figure 27.9 shows the travel times for the arcs that are used by each vehicle in the final solution.