I tried something but not sure the result is right
n = 2;
varx = Array[x, n];
vary = Array[y, n];
contsX = And @@ (a < # < b & /@ varx);
constYcontsY = And @@ (c < # < d & /@ vary);
u = 0.4; a = 0; b = 1; c = 0; d = 1;
func = Sum[y[i]^u x[i]^(1 - u), {i, n}];
NMaximize[{func, Sum[x[i], {i, n}] == 1 && Sum[y[i], {i, n}] == 1 && contsX && constYcontsY}, Join[varx, vary]]
EDITED:
varx = Array[x, n];
vary = Array[y, n];
contsXconstX = And @@ MapThread[#1 <= #2 <= #3 &, {a, varx, b}];
constY = And @@ MapThread[#1 <= #2 <= #3 &, {c, vary, d}];
func = Sum[y[i]^uTotal[y[#]^u x[i]^x[#]^(1 - u), {i,& n}];/@ Range@n]
FindMaximum[{func, Total@varx == 1 && Total@vary == 1 && contsX && constYcontsY}, Join[varx, vary]]