@Carl Love Many thanks. I run your code with a slight modification (**dim** and **Model_3**) but the outcome does not seem right. What did I miss?

restart:
MakeMonoList:= proc(alpha::symbol, V::list(symbol), dim::nonnegint)
local i, k, T;
coeffs(expand((1+add(V))^dim), V, T);
unapply(
sort([T], rcurry(Groebner:-TestOrder, 'tdeg'(V[])))
*~ [seq](alpha[i,k], k= 0..binomial(nops(V)+dim, dim)-1),
i::posint,
'proc_options'= 'remember'
)
end proc
:
M:= MakeMonoList(alpha, [z,y,x], 2);
buildAllPossibleModels:= (model::list, M::procedure)->
[for local i,m in model do
map(subs, m=~ m+~ subs({`if`}(m::`+`, op(m), m)=~ (), M(i)), model)[]
od]
:
Models_3 := [[y*z*alpha[1, 8], x*z*alpha[2, 6], x*y*alpha[3, 5]], [y*z*alpha[1, 8], x*y*alpha[2, 5], y*z*alpha[3, 8]], [y*z*alpha[1, 8], x*y*alpha[2, 5], y^2*alpha[3, 7]], [y*z*alpha[1, 8], x*y*alpha[2, 5], x*y*alpha[3, 5]], [y*z*alpha[1, 8], x^2*alpha[2, 4], x*z*alpha[3, 6]], [y*z*alpha[1, 8], x^2*alpha[2, 4], x^2*alpha[3, 4]], [y*z*alpha[1, 8], z*alpha[2, 3], x*z*alpha[3, 6]], [y*z*alpha[1, 8], x*alpha[2, 1], x*z*alpha[3, 6]], [y*z*alpha[1, 8], x*alpha[2, 1], x*alpha[3, 1]], [y^2*alpha[1, 7], z^2*alpha[2, 9], x*z*alpha[3, 6]], [y^2*alpha[1, 7], y*z*alpha[2, 8], x*z*alpha[3, 6]], [y^2*alpha[1, 7], y*z*alpha[2, 8], x^2*alpha[3, 4]], [y^2*alpha[1, 7], y*z*alpha[2, 8], x*alpha[3, 1]], [y^2*alpha[1, 7], x*z*alpha[2, 6], x*z*alpha[3, 6]], [y^2*alpha[1, 7], x*z*alpha[2, 6], x*y*alpha[3, 5]], [y^2*alpha[1, 7], x*z*alpha[2, 6], y*alpha[3, 2]], [y^2*alpha[1, 7], z*alpha[2, 3], x*z*alpha[3, 6]], [y^2*alpha[1, 7], z*alpha[2, 3], x*y*alpha[3, 5]], [x*y*alpha[1, 5], z^2*alpha[2, 9], x*z*alpha[3, 6]], [x*y*alpha[1, 5], y*z*alpha[2, 8], x*z*alpha[3, 6]], [x*y*alpha[1, 5], y*z*alpha[2, 8], x*y*alpha[3, 5]], [x*y*alpha[1, 5], y*z*alpha[2, 8], x^2*alpha[3, 4]], [x*y*alpha[1, 5], y*z*alpha[2, 8], x*alpha[3, 1]], [x*y*alpha[1, 5], x*z*alpha[2, 6], y*z*alpha[3, 8]], [x*y*alpha[1, 5], x*z*alpha[2, 6], y^2*alpha[3, 7]], [x*y*alpha[1, 5], x*z*alpha[2, 6], x*z*alpha[3, 6]], [x*y*alpha[1, 5], x*z*alpha[2, 6], x*y*alpha[3, 5]], [x*y*alpha[1, 5], x*z*alpha[2, 6], x^2*alpha[3, 4]], [x*y*alpha[1, 5], x*z*alpha[2, 6], y*alpha[3, 2]], [x*y*alpha[1, 5], x*z*alpha[2, 6], x*alpha[3, 1]], [x*y*alpha[1, 5], z*alpha[2, 3], x*z*alpha[3, 6]], [x*y*alpha[1, 5], z*alpha[2, 3], x*y*alpha[3, 5]], [x*y*alpha[1, 5], z*alpha[2, 3], x^2*alpha[3, 4]], [x*y*alpha[1, 5], z*alpha[2, 3], x*alpha[3, 1]], [y*alpha[1, 2], z^2*alpha[2, 9], x*z*alpha[3, 6]], [y*alpha[1, 2], y*z*alpha[2, 8], x*z*alpha[3, 6]], [y*alpha[1, 2], y*z*alpha[2, 8], x*y*alpha[3, 5]], [y*alpha[1, 2], y*z*alpha[2, 8], x^2*alpha[3, 4]], [y*alpha[1, 2], y*z*alpha[2, 8], x*alpha[3, 1]], [y*alpha[1, 2], x*z*alpha[2, 6], y*z*alpha[3, 8]], [y*alpha[1, 2], x*z*alpha[2, 6], y^2*alpha[3, 7]], [y*alpha[1, 2], x*z*alpha[2, 6], x*z*alpha[3, 6]], [y*alpha[1, 2], x*z*alpha[2, 6], x*y*alpha[3, 5]], [y*alpha[1, 2], x*z*alpha[2, 6], x^2*alpha[3, 4]], [y*alpha[1, 2], x*z*alpha[2, 6], y*alpha[3, 2]], [y*alpha[1, 2], x*z*alpha[2, 6], x*alpha[3, 1]], [y*alpha[1, 2], z*alpha[2, 3], x*z*alpha[3, 6]], [y*alpha[1, 2], z*alpha[2, 3], x*y*alpha[3, 5]], [y*alpha[1, 2], z*alpha[2, 3], x^2*alpha[3, 4]]]:
Models_4:= CodeTools:-Usage(
Threads:-Map(buildAllPossibleModels, Models_3, M)
):
M := proc (i::posint) option remember; [alpha[i, 0], x*alpha[i,
1], y*alpha[i, 2], z*alpha[i, 3], x^2*alpha[i, 4],
y*x*alpha[i, 5], z*x*alpha[i, 6], y^2*alpha[i, 7],
z*y*alpha[i, 8], z^2*alpha[i, 9]] end proc
memory used=0.92MiB, alloc change=19.69MiB, cpu time=16.00ms, real time=4.00ms, gc time=0ns
nops(convert(Models_4,set));
49

The number of models should be much higher than 49, more like 1170.

How can I avoid the line where the result is converted to a set?