Hi Odenwald,   Actually 631 is much better than any value I have seen.  My experience with GA's is they tend to give you a quick approximate solution in the ballpark, but finding the globally optimal is not guaranteed.  So I think your instinct to explore the MILP route is excellent if seeking a robust way to find the global solution. Our team is still exchanging e-mails on the best way to support this problem in that respect.   So finding a way to make the GA find the global might be a needle in the haystack task to find the best heuristic.  That being said, I ran out of time; in the objective evaluation, there is a local search performed seeking to reduce the total distance.  That is PROC GA allows you to do a greedy search on any evaluation call and if a better point is found, you can simply overwrite the input point with your new point, and it will be used from that point on.  A very cool feature I wasn't fully cognizant until this project. do i = 2 to ncities; city1 = s[i]; inext = 1 + mod(i,ncities); if (inext=1) then inext=2; city2 = s[inext]; if i=1 then before = s[ncities]; else before = s[i-1]; after = s[1 + mod(inext,ncities)]; if (distances[before,city1]+distances[city2,after]) > (distances[before,city2]+distances[city1,after]) then do; s[i] = city2; s[inext] = city1; end; end; call WriteMember(selected,1,s);   If you notice, it will swap cities regardless of the cycle boundary.  So it might be improved if only swapping within a cycle.   The encoding I used was using integer sequences.  So for example (if you had 2 cars and 9 cities and evenly split), the sequence  1 3 9 4 2 6 8 7 5  could be interpreted as: Route 1: 1 ->3->9->4->2 -> 1 Route 2: 1->6->8->7->5->1 So I might want to reorder 3, 9, 4, 2  and 6,8,7,5 greedily in a way to reduce the cycle distance, not the total distance for a single car.  I left it there more as example of what could be done, though it is not clear if it actually helps/hurts at this point. Another assumption I made but did not optimize is to divide the number of cities as close to equally as possible. This part is hard-coded, so the optimal solution may in fact not be obtainable with this encoding.  For example, in the above, I could use the first 5 integers after the first one to define route 1 and the last 3 for the second: Route 1: 1 ->3->9->4->2->6 -> 1 Route 2: 1->8->7->5->1 In theory, you could greedily check which split is optimal in the evaluation or figure out a way to encode this choice so the GA may choose.  Note that vehmap in the code just hard-codes the bounds where the segments start-stop for each vehicle:   /* Populate vehicle map start/stop */ array vehmap[8,2] / nosym; cnt=2; inc=2; do i=1 to 8; vehmap[i,1]=cnt; vehmap[i,2]=cnt+inc; cnt=cnt+inc+1; if (i > 4) then inc=1; end; For example, you could add to this code to see what this code is doing:     do i=1 to 8; put "Car " i ": starts at loc " vehmap[i,1] " and stops at loc " vehmap[i,2]; end; To see the output     Car 1 : starts at loc 2 and stops at loc 4 Car 2 : starts at loc 5 and stops at loc 7 Car 3 : starts at loc 8 and stops at loc 10 Car 4 : starts at loc 11 and stops at loc 13 Car 5 : starts at loc 14 and stops at loc 16 Car 6 : starts at loc 17 and stops at loc 18 Car 7 : starts at loc 19 and stops at loc 20 Car 8 : starts at loc 21 and stops at loc 22 So when the GA passes in the selected array:   function CVRPSwap(selected[*],ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]); array s[1] /nosym; call dynamic_array(s,ncities); call ReadMember(selected,1,s); The array "s" is interpreted so that s[2], s[3], s[4] are the cities visited by vehicle one, s[5], s[6], s[7], by car 2, etc.  We really only need a sequence of size 21, since 1 is never permuted.  For simplicity I left it at 22 and always permute 1 back to s[1]:     /* Permut city 1 back to front */ do i = 1 to ncities; cityi=s[i]; if (cityi = 1) then do; s[i]=s[1]; s[1]=cityi; end; end; Otherwise, a lot of the TSP loop indices would need to change.  And because we call after the local optimization step mentioned above: call WriteMember(selected,1,s);   The first string entry is always 1 on exit from all evaluations.  So it is effectively the same.   Anyways, trying to make a little less of this code black box.  But in general, an efficient MILP type solution would be the way to go to reliably achieve global solutions.             ... View more

Hi Odenwald,   I am very sorry for the slow response.  I have been following up in the background.  I am very familiar with LSO, not so much PROC GA/MILP/CVRP.  I followed up with our MILP team and they are discussing the possibility of a dedicated VRP solver similar to our TSP.  So crossing fingers that project gets a green light.    Note that LSO has been renamed "Blackbox" for future releases.     Disclaimer, not being an expert in either PROC GA or CVRP I took a stab at transforming the PROC GA TSP problem into a CVRP example.  It seems to run correctly now, but I can't be certain.  So it is a sketch of how one would do it.  It gives results to the output like the below in cycle format.    Cycle for vehicle: 1 1 -> 5 (load = 1400 ) 5 -> 4 (load = 2200 ) 4 -> 8 (load = 3000 ) 8 -> 1 (load = 3000 ) Cycle for vehicle: 2 1 -> 12 (load = 1200 ) 12 -> 14 (load = 2500 ) 14 -> 15 (load = 2800 ) 15 -> 1 (load = 2800 ) Cycle for vehicle: 3 1 -> 22 (load = 700 ) 22 -> 18 (load = 1700 ) 18 -> 13 (load = 3000 ) 13 -> 1 (load = 3000 ) Cycle for vehicle: 4   data locations; input node x y demand; datalines; 1 145 215 0 2 151 264 1100 3 159 261 700 4 130 254 800 5 128 252 1400 6 163 247 2100 7 146 246 400 8 161 242 800 9 142 239 100 10 163 236 500 11 148 232 600 12 128 231 1200 13 156 217 1300 14 129 214 1300 15 146 208 300 16 164 208 900 17 141 206 2100 18 147 193 1000 19 164 193 900 20 129 189 2500 21 155 185 1800 22 139 182 700 ; run; proc ga data1 = locations seed = 5554 lastgen=lastgen; call SetEncoding('S22'); ncities = 22; ncar = 8; cap = 3000; /* v1 = 1 -> 2 -> 3 -> 4 -> 1 (indices of S owned by vehicial 1) v2 = 1 -> 5 -> 6 -> 7 -> 1 v3 = 1 -> 8 -> 9 -> 10 -> 1 v4 = 1 -> 11 -> 12 -> 13 -> 1 v5 = 1 -> 14 -> 15 -> 16 -> 1 v6 = 1 -> 17 -> 18 -> 1 v7 = 1 -> 19 -> 20 -> 1 v8 = 1 -> 21 -> 22 -> 1 */ /* Create distance matix once up front */ array distances[22,22] /nosym; do i = 1 to 22; do j = 1 to i; distances[i,j] = sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2); distances[j,i] = distances[i,j]; end; end; /* Populate vehicle map start/stop */ array vehmap[8,2] / nosym; cnt=2; inc=2; do i=1 to 8; vehmap[i,1]=cnt; vehmap[i,2]=cnt+inc; cnt=cnt+inc+1; if (i > 4) then inc=1; end; /* Objective function with local optimization */ function CVRPSwap(selected[*],ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]); array s[1] /nosym; call dynamic_array(s,ncities); call ReadMember(selected,1,s); /* Permut city 1 back to front */ do i = 1 to ncities; cityi=s[i]; if (cityi = 1) then do; s[i]=s[1]; s[1]=cityi; end; end; /* Permute first city /* First try to improve solution by swapping adjacent cities */ do i = 2 to ncities; city1 = s[i]; inext = 1 + mod(i,ncities); if (inext=1) then inext=2; city2 = s[inext]; if i=1 then before = s[ncities]; else before = s[i-1]; after = s[1 + mod(inext,ncities)]; if (distances[before,city1]+distances[city2,after]) > (distances[before,city2]+distances[city1,after]) then do; s[i] = city2; s[inext] = city1; end; end; call WriteMember(selected,1,s); /* Calculate cost and max load */ cost=0; loaderr=0; /* one norm of contraint violation */ do c = 1 to ncar; load=0; city_from=1; do loc = vehmap[c,1] to vehmap[c,2]; city_to = s[loc]; load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */ cost = cost + distances[city_from, city_to]; city_from=city_to; end; city_to=1; cost = cost + distances[city_from, city_to]; loaderr = loaderr+max(0, load-cap); /* need load <= 3000 */ end; merit = cost + 1e10*loaderr; /* put merit; */ return(merit); endsub; /* Called after each iteration */ subroutine logReport(itercnt, ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]); outargs itercnt; itercnt=itercnt + 1; array s[1,22] / nosym; call GetSolutions(s, 1, 1); /* Calculate cost and max load */ cost=0; loaderr=0; /* one norm of contraint violation */ do c = 1 to ncar; load=0; city_from=1; do loc = vehmap[c,1] to vehmap[c,2]; city_to = s[1,loc]; load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */ cost = cost + distances[city_from, city_to]; city_from=city_to; end; city_to=1; cost = cost + distances[city_from, city_to]; loaderr = max(0, load-cap); /* need load <= 3000 */ end; merit = cost + 1e10*loaderr; put "iter " itercnt ": " s merit; endsub; /* Called when finished to see solution */ subroutine printSol(itercnt, ncities,ncar, cap, distances[*,*], demand[*], vehmap[*,*]); outargs itercnt; itercnt=itercnt + 1; array s[1,22] / nosym; call GetSolutions(s, 1, 1); /* Calculate cost and max load */ cost=0; loaderr=0; /* one norm of contraint violation */ do c = 1 to ncar; load=0; put "Cycle for vehicle: " c; city_from=1; do loc = vehmap[c,1] to vehmap[c,2]; city_to = s[1,loc]; load = load + demand[city_to]; /* demand on 1 is 0, so okay to add on return as well */ cost = cost + distances[city_from, city_to]; put city_from " -> " city_to " (load = " load ")"; city_from=city_to; end; city_to=1; cost = cost + distances[city_from, city_to]; put city_from " -> " city_to " (load = " load ")"; loaderr = max(0, load-cap); /* need load <= 3000 */ end; merit = cost + 1e10*loaderr; put "iter " itercnt ": " s merit; endsub; itercnt=0; call SetObjFunc('CVRPSwap',0); call SetUpdateRoutine('logReport'); call SetFinalize('printSol'); call SetCross('Order'); call SetMut('Invert'); call SetMutProb(0.05); call SetCrossProb(0.8); call SetElite(1); call Initialize('DEFAULT',200); call ContinueFor(100); run; proc print data=lastgen; run;   Not sure if this is helpful as I realize it may not be easy to understand what is going on.  I am borrowing time from other projects, so I probably cannot develop this code more. In general, we do need a dedicated solver for this problem type (which would have heuristics); and so I hope we get clearance to work on in the near future.   ... View more

Wrt LSO:   Applying a GA to CVRP’s requires using sequence variable encoding and/or customized permutation/crossover functions.  The corresponding dimension of the MILP problem tends to have too many variables and constraints for a GA to be applicable in general-purpose form.   For example, consider the TSP formulation for MILP: https://go.documentation.sas.com/?docsetId=casmopt&docsetTarget=casmopt_milpsolver_examples04.htm&docsetVersion=8.5&locale=en versus the TSP formulation for a GA: https://go.documentation.sas.com/?docsetId=orlsoug&docsetTarget=orlsoug_ga_examples01.htm&docsetVersion=15.1&locale=en   Unfortunately, LSO does not support sequence variables or custom crossover/mutation functions, and would thus not work for this problem type.  In general, a GA that treats the problem as a Blackbox will have trouble handling integer variable constraints and/or variable dimensions exceeding 100; by working with sequence variables, the constraints would be implicitly satisfied for all permutations and crossovers while the variable dimension is dramatically reduced.  Similarly, custom crossover/mutation operators could simply jump from one feasible point to the next.   It might be possible to try something similar using PROC GA which does support both of these features.  Invoking the PROC OPTMODEL code directly may be challenging.   As of yet, we do not support Tabu search or Simulated Annealing (SA).  Adding support for both sequence variables and SA has been regularly discussed but is not yet on target for development, unfortunately. ... View more

This is something we could do for you in the next release.  But for now would need to be done manually:   One potential heuristic to overcome this hurdle if N "units" is greater than 2 is to create a binary LHS sampling (over which variables should be fixed at 0 and which should be free) that predetermines which variables will get fixed for a set of B "tests".  There is example code here: https://www.linkedin.com/pulse/latin-hypercube-sampling-emily-gao/    I.e. you could create an LHS sample over the unit-hypercube that is stored in data set "T" that is B x N and read into optmodel. For the ith row of T you would solve a nlp problem where Qj is fixed at 0 if Tij <= .5 and  So how many NLP problems you solve is controllable and a function of the number of "experiments" you decide to generate.    It is of course a heuristic ... ... View more

Thanks a lot for the two examples!  I was able to get similar solutions (using sledgehammer sized population:) if I used these two solve commands:   solve with lso / absfconv=1e-6 nabsfconv=10 maxgen=10 popsize=10000 feastol=1e-2; solve with nlp /algorithm=as;   But this could be a game of wack-a-mole.  Also it is much slower ... so likely a dead-end unless no 2 does not help.   Thanks a lot for the examples.  For better or worse, LSO treats all problems as black-boxes and avoids assumptions about derivatives existing.  So it may continue to optimize when other solvers get stuck--the price is when derivatives do exist it will still tick-tack to the solution.  Thus in general, if everything insight is smooth, the nlp solvers are more appropriate.  But adding binary variables will of course induce some (albeit structured) non-smoothness ...   ... View more

That is unfortunate.  It may be LSO is not a good fit.   I did try the first problem to verify LSO gave comparable results, but did not have time to run through the entire list.  If you had an instance in the loop where the objective Q.sol is much worse, that would help to diagnose.   Looking closer, I didn't realize the constraints were nonlinear.  The feastol is likely too tight as LSO is a DFO algorithm.   If you have time/interest, could you rerun but replace the solve command with:   solve with lso / absfconv=1e-6 nabsfconv=100 maxgen=100 popsize=500 feastol=1e-3; solve with nlp;   The second nlp solve is to increase the feasibility of the solution found.  When moving to binary we can can add a relaxint option that may work analogously for the mixed integer case.   Thanks.   ... View more

Could you please try the experiment with these options:   solve with lso / absfconv=1e-10 nabsfconv=100 maxgen=500 popsize=500 feastol=1e-7;   Increasing the popsize often improves the solution quality.   ... View more