Thanks Neil.

It's definitely a "one-off' thing --- nothing commercial at stake here.

I played a bit with heuristics akin to what you mention. My progress was so slow I felt the method might not converge. Perhaps I should go back to this mixed, user-guided, approach, and manually tweaks routes until I find a set that works.

Alternatively, there's implementing a published algorithm. Unfortunately papers that seem relevant (e.g. references 3 to 10 in the linked dissertation above) are behind a paywall...

Hi Ted,

Yes, I had seen that. This is solving a different problem because each route is a tour. Ie. some kids would be picked up when the bus is on its way back from the school to the depot.

Or perhaps there is a way to transform the problem so it can be solved using the demonstration's code? I need to try making all costs from stops to depots infinite, and costs from school/destination to bus depots all equal to zero.

This is useful information, Neil.

I ended up with this simple, non-optimal, algorithm: 1. Cluster all stops according to their nearest depot 2. Within each cluster, start a route taking the furthest stop and greedily adding the nearest stops until the route does not satisfy the requirements any more.

travelTimes[route_List] := 
  Reverse@Accumulate@Reverse[Append[
      distStops[[#[[1]], #[[2]]]] & /@ Partition[route, 2, 1], 
      distStopsDestination[[Last@route]]]];

routeOKQ[route_List] := With[{\[Alpha] = 1.5}, 
   And @@ MapThread[#1 < \[Alpha]*#2 &, {travelTimes@route, distStopsDestination[[route]]}]];

findRoutes[indices_List] := Module[
   {pool, route, routes, i0, r, c},
   pool = indices;
   routes = {};
   While[Length@pool > 0,
       route = {pool[[First@Ordering[distStopsDestination[[pool]], -1]]]};
        pool = DeleteCases[pool, Alternatives @@ route];
        While[Length@pool > 0,
             c = pool[[First@Ordering[distStops[[route[[-1]], pool]], 1]]];
             r = Append[route, c];
             If[routeOKQ@r,
                  route = r;
                  pool = DeleteCases[pool, Alternatives @@ route]
                  ,
                  Break[]
              ]
         ];
        AppendTo[routes, route];
        pool = DeleteCases[pool, Alternatives @@ route];
    ];
   routes
   ];

Example:

nstops = 50;
ndepots = 4;

destination = {0, 0};
depots = Table[ReIm[(1 + RandomReal[.1*{-1, 1}])*Exp[I*2*(j*Pi/ndepots + RandomReal[.1*{-1, 1}])]], 
                {j, 0, ndepots - 1}];
stops = RandomReal[{-1, 1}, {nstops, 2}];

distStops = Table[If[i == j, 0, EuclideanDistance[stops[[i]], stops[[j]]]*RandomReal[{0.9, 1.2}]],
 {i, nstops}, {j, nstops}];

distDepotsStops = Table[EuclideanDistance[depots[[i]], stops[[j]]]*RandomReal[{0.9, 1.2}], {i, ndepots}, {j, nstops}];

distStopsDestination = Table[EuclideanDistance[stops[[i]], destination]*RandomReal[{0.9, 1.2}], {i, nstops}];

nf = Nearest[Range@ndepots, DistanceFunction -> (distDepotsStops[[#2, #1]] &)];
c = Flatten@nf[Range@nstops, 1];
allroutes = Flatten[Table[findRoutes[Flatten@Position[c, k]], {k, ndepots}], 1];

ListPlot[stops, Prolog -> {
   Table[{Hue[i/Length@allroutes], 
     Line[Append[stops[[allroutes[[i]]]], destination]]}, {i, 
     Length@allroutes}],
   PointSize[Large], Red, Point[destination], Black, Point[depots]
   }, PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, AspectRatio -> 1]

Attachments: