epsilon = 2;
eq1[y_, z_] :=
  1/(4 (1 - xx)^4) + (
   13 (1 - xx) (-1 + 3 xx))/(13 (1 - xx)^2 xx - 0.5 y - z)^2;
f[y_?NumericQ, z_?NumericQ] :=
 Reduce[eq1[y, z] == 0 && 0 < xx < 1, {xx}, Reals][[2]]
c1 = (1 - x)*0.5;
c2 = (1 - f[y, z])*0.5 + y;
g1 = 1.5*epsilon^epsilon*f[y, z]*(1 - f[y, z])^epsilon;
g2 = 1.5*epsilon^epsilon*f[y, z]*(1 - f[y, z])^epsilon - 0.5*y - z;
rho = 2;
Y2[y_?NumericQ, z_?NumericQ] := 1.125*epsilon^epsilon *(1 - f[y, z]);

NMaximize[{(c1^(1 - rho) - 1)/(1 - rho) + (g1^(1 - rho) - 1)/(
   1 - rho) + (c2^(1 - rho) - 1)/(1 - rho) + (g2^(1 - rho) - 1)/(
   1 - rho),
  1.5*epsilon^epsilon*f[y, z]*(1 - f[y, z])^epsilon - 0.5*y - z > 0 &&
    y + z <= Y2[y, z] && 0 < x < 1 && 0 <= y < Y2[y, z] &&
   0 <= z < Y2[y, z]}, {x, y, z}]

{-1.97619, {x -> 0., y -> 0.508291, z -> 0.}}

NMaximize[{(c1^(1 - rho) - 1)/(1 - rho) + (g1^(1 - rho) - 1)/(
   1 - rho) + (c2^(1 - rho) - 1)/(1 - rho) + (g2^(1 - rho) - 1)/(
   1 - rho),
  1.5*epsilon^epsilon*f[y, z]*(1 - f[y, z])^epsilon - 0.5*y - z > 0 &&
    y + z <= 4.5 && 0 < x < 1 && 0 <= y < 4.5 && 0 <= z < 4.5}, {x, y,
   z}]

{5.98763*10^11, {x -> 0.40521, y -> 0.982316, z -> 0.388413}}

NMaximize::cvdiv: Failed to converge to a solution. The function may be unbounded.