thx for your efforts Roman

Fantastic.
Does the method work when the solutions are probably not exact?
Example for n teams (n<9) the following generates the objective function to minimize and the constraints:
restart
with(Statistics)
Pl := [2, 6, 7, 8, 15, 16, 19, 21, 23, 25, 25, 27]
k := Count(Pl)
TeamNo := 1/4*k
q := [a, b, c, d, e, f, g, h]
T := seq(add(Pl[n]*seq(q[m][n], n = 1 .. k)[n], n = 1 .. k), m = 1 .. 8)
Obj := add((1/4*T[n]-add(1/4*T[n], n = 1 .. TeamNo)/TeamNo)^2, n = 1 .. TeamNo)
#Mind you I could try a variant on Obj to get rid of the squared function:
# Obj2 := add(abs(1/4*T[n]-add(1/4*T[n], n = 1 .. TeamNo)/TeamNo), n = 1 .. #TeamNo)
ConstCol := seq(add(seq(q[m][n], n = 1 .. k)[n], n = 1 .. k) = 4, m = 1 .. TeamNo)
ConstRow := seq(add(q[m][n], m = 1 .. TeamNo) = 1, n = 1 .. k)
ConstBin := seq(seq(seq(q[m][n], n = 1 .. k)[n]*(seq(q[m][n], n = 1 .. k)[n]-1) = 0, n = 1 .. k), m = 1 .. TeamNo)
F := expand(map(lhs-rhs, [ConstCol, ConstRow, ConstBin]))
#And I duplicated the below commands, modifying procedure M (please #check)
#Basically the prime decomposition step takes forever. Maybe my #computer needs upgrading… I gave up getting a solution
# factor the system of constraints
with(PolynomialIdeals):
infolevel[GroebnerBasis] := 4:
J := PolynomialIdeal(F):
P := Simplify(PrimeDecomposition(J));
# evaluate the function and sort by value
S := seq(solve(Generators(i)), i=P);
M := sort([seq([subs(i,Obj), i], i=S)], proc(a,b,c) evalb(a[1] < b[1]) end proc):
M[1]; # min
M[-1]; # max

Fantastic.
Does the method work when the solutions are probably not exact?
Example for n teams (n<9) the following generates the objective function to minimize and the constraints:
restart
with(Statistics)
Pl := [2, 6, 7, 8, 15, 16, 19, 21, 23, 25, 25, 27]
k := Count(Pl)
TeamNo := 1/4*k
q := [a, b, c, d, e, f, g, h]
T := seq(add(Pl[n]*seq(q[m][n], n = 1 .. k)[n], n = 1 .. k), m = 1 .. 8)
Obj := add((1/4*T[n]-add(1/4*T[n], n = 1 .. TeamNo)/TeamNo)^2, n = 1 .. TeamNo)
#Mind you I could try a variant on Obj to get rid of the #squared function:
# Obj2 := add(abs(1/4*T[n]-add(1/4*T[n], n = 1 .. TeamNo)/TeamNo), n = 1 .. TeamNo)
ConstCol := seq(add(seq(q[m][n], n = 1 .. k)[n], n = 1 .. k) = 4, m = 1 .. TeamNo)
ConstRow := seq(add(q[m][n], m = 1 .. TeamNo) = 1, n = 1 .. k)
ConstBin := seq(seq(seq(q[m][n], n = 1 .. k)[n]*(seq(q[m][n], n = 1 .. k)[n]-1) = 0, n = 1 .. k), m = 1 .. TeamNo)
F := expand(map(lhs-rhs, [ConstCol, ConstRow, ConstBin]))
#And I duplicated the below commands, modifying procedure M (please check)
#Basically the prime decomposition step takes forever. Maybe my #computer needs upgrading… I gave up getting a solution
# factor the system of constraints
with(PolynomialIdeals):
infolevel[GroebnerBasis] := 4:
J := PolynomialIdeal(F):
P := Simplify(PrimeDecomposition(J));
# evaluate the function and sort by value
S := seq(solve(Generators(i)), i=P);
M := sort([seq([subs(i,Obj), i], i=S)], proc(a,b,c) evalb(a[1] < b[1]) end proc):
M[1]; # min
M[-1]; # max

My apologies.
The objective function should be
> Obj := (1/4*a[1]+3/4*a[2]+7/8*a[3]+a[4]+15/8*a[5]+2*a[6]+19/8*a[7]+21/8*a[8]-1/4*b[1]-3/4*b[2]-7/8*b[3]-b[4]-15/8*b[5]-2*b[6]-19/8*b[7]-21/8*b[8])^2+(1/4*b[1]+3/4*b[2]+7/8*b[3]+b[4]+15/8*b[5]+2*b[6]+19/8*b[7]+21/8*b[8]-1/4*a[1]-3/4*a[2]-7/8*a[3]-a[4]-15/8*a[5]-2*a[6]-19/8*a[7]-21/8*a[8])^2;
Thanks to roman_pearce I have rewitten the binomal constraint 4. Also I forgot I needed another binomial constraint- containt 5
with constraints
> Constaint1:=add(a[n], n = 1 .. 8) = 4;
> Constaint2:=add(b[n], n = 1 .. 8) = 4;
> Constaint3:=seq(a[n]+b[n] = 1, n = 1 .. 8);
> Constaint4:=seq(a[n]*(a[n]-1) = 0, n = 1 .. 8);
> Constaint5:=seq(b[n]*(b[n]-1) = 0, n = 1 .. 8);
> with(Optimization)
> NLPSolve(Obj, {Constraint1, Constraint2, Constraint3, Constaint4, Constraint5}, assume = nonnegative)
Error, (in Optimization:-NLPSolve) constraints must be specified as a set or list of equalities and inequalities
Still get the same error. Constraint4 and 5 are the problems.
What I want is an output that is binomial, eg
> evalf(subs(a[1] = 1, a[2] = 0, a[3] = 0, a[4] = 1, a[5] = 0, a[6] = 1, a[7] = 0, a[8] = 1, b[1] = 0, b[2] = 1, b[3] = 1, b[4] = 0, b[5] = 1, b[6] = 0, b[7] = 1, b[8] = 0, Obj));
Believe it or not the above effort is to work out a better method to allocate table tennis teams from n players (n divisible by 4): it’s about minimizing the variance between the teams. Using the brute force and blind ignorance approach, I set up an excel spreadsheet with a whole stack of binomial indicator variables -either in the team (1) or not (0). This method works for teams up to 3. Beyond that the combinations increase according to C(n,4) and Excel churns away forever. Hence the above Maple approach.

thanks. a good resource that pdf file