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); constY = 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 && constY}, Join[varx, vary]] **EDITED:** varx = Array[x, n]; vary = Array[y, n]; contsX = And @@ MapThread[#1 <= #2 <= #3 &, {a, varx, b}]; constY = And @@ MapThread[#1 <= #2 <= #3 &, {c, vary, d}]; func = Sum[y[i]^u x[i]^(1 - u), {i, n}]; FindMaximum[{func, Total@varx == 1 && Total@vary == 1 && contsX && constY}, Join[varx, vary]]