The calculation had finished by the time I went to have another look.
The solutions were parameterized. I have attached a .mat with the values. (To get the conditions in the form I present them here, simplify() with 'steps', 100)
x1 -> z
x2 -> z1
x3 -> z2
x4 -> z3
y1 -> z4
y2 -> z5
Four independent solutions were created:
Set #1 (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)
0 <= z
0 <= z1
0 <= z2
0 <= z3
z4 == 0
0 <= z5
K22*epsilon2 + epsilon2*z5 == K22*beta12*f12*z + K22*beta22*f22*z1 + K22*beta32*f32*z2 + K22*beta42*f42*z3
s1 == z*(f12*z5 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f22*z5 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f32*z5 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f42*z5 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))
Set #2: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)
0 <= z
0 <= z1
0 <= z2
0 <= z3
0 <= z4
z5 == 0
K11*epsilon1 + epsilon1*z4 == K11*beta11*f11*z + K11*beta21*f21*z1 + K11*beta31*f31*z2 + K11*beta41*f41*z3
s1 == z*(f11*z4 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f21*z4 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f31*z4 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f41*z4 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))
Set #3: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)
N ~= 0
0 <= z
0 <= z1
0 <= z2
0 <= z3
z4 == 0
z5 == 0
N*s1 == r1*z*(eta11*z - N + eta12*z1 + eta13*z2 + eta14*z3)
N*s2 == r2*z1*(eta21*z - N + eta22*z1 + eta23*z2 + eta24*z3)
N*s3 == r3*z2*(eta31*z - N + eta32*z1 + eta33*z2 + eta34*z3)
N*s4 == r4*z3*(eta41*z - N + eta42*z1 + eta43*z2 + eta44*z3)
Set #4: (that is, x1 x2 x3 x4 y1 y2 are the tuples of values [z z1 z2 z3 z4 z5] such that the below conditions all hold)
0 <= z
0 <= z1
0 <= z2
0 <= z3
0 <= z4
0 <= z5
K11*epsilon1 + epsilon1*z4 + K11*g12*z5 + epsilon1*omega2*z5 == K11*beta11*f11*z + K11*beta21*f21*z1 + K11*beta31*f31*z2 + K11*beta41*f41*z3
K22*epsilon2 + epsilon2*z5 + epsilon2*omega1*z4 + K22*g12*gamma12*z4 == K22*beta12*f12*z + K22*beta22*f22*z1 + K22*beta32*f32*z2 + K22*beta42*f42*z3
s1 == z*(f11*z4 + f12*z5 + r1*((eta11*z + eta12*z1 + eta13*z2 + eta14*z3)/N - 1))
s2 == z1*(f21*z4 + f22*z5 + r2*((eta21*z + eta22*z1 + eta23*z2 + eta24*z3)/N - 1))
s3 == z2*(f31*z4 + f32*z5 + r3*((eta31*z + eta32*z1 + eta33*z2 + eta34*z3)/N - 1))
s4 == z3*(f41*z4 + f42*z5 + r4*((eta41*z + eta42*z1 + eta43*z2 + eta44*z3)/N - 1))
Code:
syms N x1 x2 x3 x4 y1 y2 real
assume(N >= 0 & x1 >= 0 & x2 >= 0 & x3 >= 0 & x4 >= 0 & y1 >= 0 & y2 >= 0)
syms r1 r2 r3 r4 positive
syms s1 s2 s3 s4 positive
syms epsilon1 epsilon2 omega1 positive
syms omega2 K11 K22 positive
syms beta11 beta21 beta31 beta41 positive
syms beta12 beta22 beta32 beta42 gamma12 positive
syms eta11 eta12 eta13 eta14 positive
syms eta21 eta22 eta23 eta24 positive
syms eta31 eta32 eta33 eta34 positive
syms eta41 eta42 eta43 eta44 positive
syms R c1 c2 c3 c4 positive
syms f11 f12 f21 f22 f31 f32 f41 f42 g12 positive
F1=R-c1*x1-c2*x2-c3*x3-c4*x4;
F2=x1*(r1*(1-(eta11*x1+eta12*x2+eta13*x3+eta14*x4)/N)-f11*y1-f12*y2)+s1==0;
F3=x2*(r2*(1-(eta21*x1+eta22*x2+eta23*x3+eta24*x4)/N)-f21*y1-f22*y2)+s2==0;
F4=x3*(r3*(1-(eta31*x1+eta32*x2+eta33*x3+eta34*x4)/N)-f31*y1-f32*y2)+s3==0;
F5=x4*(r4*(1-(eta41*x1+eta42*x2+eta43*x3+eta44*x4)/N)-f41*y1-f42*y2)+s4==0;
F6=y1*(-epsilon1*(1+(y1+omega2*y2)/K11)-g12*y2+beta11*f11*x1+beta21*f21*x2+beta31*f31*x3+beta41*f41*x4)==0;
F7=y2*(-epsilon2*(1+(omega1*y1+y2)/K22)-gamma12*g12*y1+beta12*f12*x1+beta22*f22*x2 +beta32*f32*x3+beta42*f42*x4)==0 ;
eqns=[F2,F3,F4,F5,F6,F7];
S = solve(eqns,[x1,x2,x3,x4,y1,y2], 'returnconditions', true);
