@Carl Love thanks for your effort here.

 Sorry to keep bugging you, How does one choose an initial variance?

I ran your code for Var(P)>0.01 (any lower and the kernel crashed).

I got these results:

Var(P);

0.89246398631e-2

P; [{1236.0, 1292.0, 1404.0, 1829.0}, {1249.0, 1392.0, 1477.0, 1594.0}, {1277.0,1325.0, 1477.00, 1644.0}, {1297.0, 1313.0, 1520.0, 1576.0}]

add(abs(convert(P[j], `*`)-convert(P[j+1], `*`)),j=1..nops(P)-1)+abs(convert(P[nops(P)], `*`)-convert(P[1], `*`));

73085876073.82

 In constrast, Evolver got these results [after a couple of minutes running]....

P:=[{1829,1404,1292,1236},{1644,1477,1325,1277},{1576,1520,1313,1297},{1594,1477,1392,1249}];

add(abs(convert(P[j], `*`)-convert(P[j+1], `*`)),j=1..nops(P)-1)+abs(convert(P[nops(P)], `*`)-convert(P[1], `*`));

58155331960

(smaller)

@Carl Love ran this code:

restart:
S := [1829.0, 1644.0, 1594.0, 1576.0, 1520.0, 1477.0, 1477.0001,1404.0, 1392.0, 1325.0, 1313.0, 1297.0, 1292.0, 1277.0, 1249.0, 1236.0]:
SL:=[seq(A||i,i=1..16)]:
assign(Labels ~ (S) =~ SL); #Create remember table.
nops(%%);

ts := time():  AllP:= [seq(P, P= Iterator:-SetPartitions({op(S)},[[4,4]],compile= false))]:
lnp:= evalf(ln(`*`(S[])^(1/3))):
Var:= proc(P::list(set))
local r:= evalf(`+`(map(b-> abs(ln(`*`(b[]))-lnp), P)[]));
end proc:
ans:=sort(AllP, (x,y)-> Var(x) < Var(y))[1];
 time()-ts;

subsindets(ans, realcons, Labels);  

@Carl Love silly me, I meant 10^20

@Joe Riel thanks Joe and Acer, thats just what I wanted

@Carl Love 

When I try to order the set [1829.0, 1644.0, 1594.0, 1576.0, 1520.0, 1477.0, 1477.0, 1404.0, 1392.0,

   1325.0, 1313.0, 1297.0, 1292.0, 1277.0, 1249.0, 1236.0] into 4 groups of 4, the program takes (understandably) a long time. Rather than minimizing the variance, is it possible to reduce the time by setting the variance<= a predetermined value (like 0.021).

real_labels.mw
@Kitonum Could you kindly correct attached worksheet so it works out the 3 groups and attributes labels correctly, thanks.

[It needs to work with lists instead of sets, and real numbers instead of integers].

@Kitonum for those of us unenlightened can you show how you derived your formula?

@Carl Love 

@Carl Love thanks.

yes, there are 15 throws (out of 36 states)  which satisfy the given sums (=5/12).

if the number pairs can be generated using combinat:-permute(1,1,2,2,3,3,4,4,5,5,6,6)

would you kindly write some code that will count the pairs which sum to a given number?

eg sum =7-> n=6;  sum=8 -> n=5

play_rank.mw
rank.xls
@Carl Love Sorry to keep bugging you.

When I imported this file into Maple it reported 0's (maybe because the numbers are formula generated). but when i copy-paste-values into another column =, the numbers are impoted but with rabbits ears.

@Carl Love Thanks again.

tested well, only slight problem to report. if S contains a duplicate, then the labeling is confused. eg in the above list replace 30 by 28 and you end up with two H's and no G's.

@Carl Love 

@Carl Love thanks, thats just what I wanted.

An extra Q: If the original set S corresponds to A,B,C...I, (ie A=3,B=4....) is it possible to map the final result back to the corresponding letters. in this case, [{3, 9, 35}, {4, 8, 28}, {5, 6, 30}] = [{A,F,I},{B,.......]

@Kitonum thanks for your efforts.

What about this set :{3, 4, 5, 6, 8, 9, 28, 30, 35}?

Using Evolver I got {5,6,30},{4,8,28}and {3,9,35} as the three groups close to each other as possible

@Carl Love good-o.

its a pity one can't put the character font (eg Bold,Brown Aerial,size=10) into it