Simply submitting the relevant query to msolve will do.

{x = 26730899*a[2]+1162213*a[5]+507520*a[9]+808519*a[4]+757969*a[8]+252655*a[3]+12167*a[1]+529*a[7], y = 8046000599*a[2]+748465172*a[5]+3822524979*a[9]+8821208032*a[4]+588083951*a[8]+5934310109*a[3]+8346050986*a[1]+7538374844*a[7]+4385029649*a[6]};

{a[5] = 0, a[2] = 0, a[3] = 0, a[4] = 0, a[8] = 0, a[1] = 0, a[6] = 6880584, a[9] = 4174, a[7] = -4004513};

Do you know your expression is annihilated by P:

P := 1+256*u*(256*u^3-135)*lambda^4+(16384*u^3+1728)*lambda^3+1536*lambda^2*u^2+64*lambda*u;
`@`(evala, Normal, subs)(lambda = expr, P);
 

I would guess you are meant to supply irreducible forms. You can always try this:

(degree, factor)(evala(Norm(lambda-bSol[1])));
Digits:=150;
lcm(op(map(p->PolynomialTools:-MinimalPolynomial(evalf(p), x, 36), bSol)));

Or start with this:

PA := evala(Algfield(bSol)):
PA[3];
 

If you map objective via exp function the following works naturally:

Optimization:-Maximize(x[1]*x[2]*x[3], {0 <= x[1], 0 <= x[2], 0 <= x[3], x[1]+x[2] <= 1, x[1]+x[3] <= 2});

I did not anticipate this limitation of Maple's solve. Given these responses, ?solve, parametric is the better choice.
Example:

solve({0 < 4*t1+4*r1, t1^2-4*t2 = r1^2, 0 <= r1}, {t2, r1}, parametric, parameters = {t1});
solve({0 < 4*t1+4*r1, t1^2-4*t2 = r1^2, 0 <= r1}, {t2, t1}, parametric, parameters = {r1});
 

Whenever dealing with radicals, use radnormal:

e1:= -(b*(b^2 - 4*a*c)^(1/2) - 4*a*c + b^2)/(2*a*(b^2 - 4*a*c)^(1/2))
e2:= -(b + (b^2 - 4*a*c)^(1/2))/(2*a);

radnormal(e1, rationalized);
radnormal(e1 - e2, rationalized);

From the help file: The roots function computes the exact roots of a polynomial over the rationals or an algebraic number field. No results means there are no rational roots in your example. I am quite certain roots works only with algebraic constants.

Examples drawn from the help file:
roots(x^4-4, x, {});
roots(x^4-4, x, {sqrt(2)});
roots(x^4-4, x, {sqrt(2), I});
roots((x-a)*(x^4-4), x, {sqrt(2), I});

The first will give you only rational roots, second rationals extended by sqrt(2), and third rationals extended by sqrt(2) and I.

ListTools:-Classify(degree, [x-y, y+z, x^3-x*y*z+z^3, y^3-x*z^2, y^6-x^4*y^2, x^5*y-z^6+x^2*y^2*z^2], {x, y, z})

You can also apply convert,list to the result.

This works, though it is not efficient:

rules := [g = d, e = a, a = b, f = d, c = g, b = b, d = c, c = f, a = e];
gs:=proc(S) local i; member(op(1, sort(S)), S, 'i'); [op(i..-1, S), op(1..i-1, S)]; end:
hr:=proc(n, rules)
local e:=n, S:=n, i, j, C;
do
e := subs(rules, e);
if member(e, [S], 'i') then
C := gs([op(i..-1, [S])]);
return seq(op(j, [S]) = ([op(j..i-1, [S])], op(i,[S]), C), j = 1..i-1), seq(op(j, [S]) = ([], op(j, [S]), C), j = i..nops([S]));
else S:=S, e; fi; 
od;
end:
temp:=`@`({op}, map)(op, rules); R:=NULL: 
while temp<>{} do
S := hr(op(1, temp), rules);
temp := temp minus map(lhs, {S});
R:=R, S;
od:
map(lprint, {R}):
cycle:=map(E->rhs(E)[3], {R});
entry:=map(E->rhs(E)[2], {R});
trail:=map(E->rhs(E)[1], {R});
 

The roots of this polynomial can be obtained as follows:

poly:=x^6-3*x+3;
polyroots:=map2(op, 1, roots(poly, proc () local result, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13; global _Z; k6 := RootOf(_Z^6-3*_Z+3); k2 := 6*k6; k5 := k6^2; k9 := k5^2; k3 := k6*k9; k4 := k6*k5; k11 := -15+6*k3+10*k4+12*k5+14*k9; k10 := k2-k11; k8 := _Z^2; k7 := 17*k8; k12 := 3*k3+7*k9+5*k4+6*k5-3*k6; k13 := k7-33+k12; k1 := RootOf(_Z*k13+k8*k10+k12+18); result := {RootOf(k10*_Z+(17*_Z+17*k1+k10)*k1+k13), RootOf(k7+(11*k6+k11)*_Z+11*k3+3*k9+7*k4+5*k5+k2-36)} end proc()));
 

Will this suffice?
 

'sum(
         (A-B)
         *(-1+binomial(N-2,i))
         *(A)^(i)
         *(B)^(N-2-i)
     ,i=0..N-2)';

simplify(evalindets('%',specfunc(anything, binomial), convert, GAMMA));
 

g & h can access z because you passed the variable to them in the parameters.

You are trying to solve EQ rather than eq:

eq := sqrt(x^2-10*x+1) = sqrt(-8*x^2+9*x-1);
EQ := (lhs, rhs)(map(Im=0, eq)), eq;
sol__eq := [solve]({eq}, {x});
sol__EQ := [solve]({EQ}, {x});
map((radnormal@subs), sol__eq, [EQ]);
map2(map, is, %);
 

`type/thatsum`:='`+`'(anything);
f:=x->mul([op](x));
applyrule((A::thatsum)='f(A)', [a+b,c+d,a+b+c,c+e+f+g+h]);
applyrule((A::'''`+`'''(anything))='f(A)', [a+b,c+d,a+b+c,c+e+f+g+h]);

minimize(expr, a=0..infinity, b=0..infinity, c=0..infinity, d=0..infinity, location);

Gives the same 4 solutions.