Finding one non-negative solution to a system of linear equations

There is one and only solution in non-negative reals and that is probably why FindInstance has problems with it. FindInstance is good for cases where there are lots of solutions and you need few of them.

If you use Solve you get the solution in no time.

Solve[y[9] == 0 && y[8] == 1 - 2 y[12] && y[7] == 3 - 2 y[11] - 3 y[14] && y[6] == -2 y[10] - 3 y[13] - 4 y[15] && y[5] == 0 && y[4] == 1 && y[3] == 1 + y[12] && y[2] == y[11] + 2 y[14] && y[1] == 1 + y[10] + 2 y[13] + 3 y[15] && x[15] == 0 && x[14] == 1 + y[15] && x[13] == -y[15] && x[12] == y[14] - y[15] && x[11] == 2 + y[13] - y[14] + 3 y[15] && x[10] == -y[13] - 2 y[15] && x[9] == y[12] - y[14] && x[8] == y[11] - y[12] - y[13] + 3 y[14] - 2 y[15] && x[7] == 3 + y[10] - y[11] + 3 y[13] - 2 y[14] + 5 y[15] && x[6] == -y[10] - 2 y[13] - 3 y[15] && x[5] == -y[12] && x[4] == 1 - y[11] + y[12] - 2 y[14] && x[3] == -1 - y[10] + y[11] - 2 y[13] + 2 y[14] - 3 y[15] && x[2] == 1 + y[10] + 2 y[13] + 3 y[15] && x[1] == 0, NonNegativeReals] 

{{x[1] -> 0, x[2] -> 1, x[3] -> 0, x[4] -> 0, x[5] -> 0, x[6] -> 0, x[7] -> 2, x[8] -> 1, x[9] -> 0, x[10] -> 0, x[11] -> 2, x[12] -> 0, x[13] -> 0, x[14] -> 1, x[15] -> 0, y[1] -> 1, y[2] -> 1, y[3] -> 1, y[4] -> 1, y[5] -> 0, y[6] -> 0, y[7] -> 1, y[8] -> 1, y[9] -> 0, y[10] -> 0, y[11] -> 1, y[12] -> 0, y[13] -> 0, y[14] -> 0, y[15] -> 0}} 

Does this suffice?

eqs = List @@ (y[9] == 0 && y[8] == 1 - 2 y[12] && y[7] == 3 - 2 y[11] - 3 y[14] && y[6] == -2 y[10] - 3 y[13] - 4 y[15] && y[5] == 0 && y[4] == 1 && y[3] == 1 + y[12] && y[2] == y[11] + 2 y[14] && y[1] == 1 + y[10] + 2 y[13] + 3 y[15] && x[15] == 0 && x[14] == 1 + y[15] && x[13] == -y[15] && x[12] == y[14] - y[15] && x[11] == 2 + y[13] - y[14] + 3 y[15] && x[10] == -y[13] - 2 y[15] && x[9] == y[12] - y[14] && x[8] == y[11] - y[12] - y[13] + 3 y[14] - 2 y[15] && x[7] == 3 + y[10] - y[11] + 3 y[13] - 2 y[14] + 5 y[15] && x[6] == -y[10] - 2 y[13] - 3 y[15] && x[5] == -y[12] && x[4] == 1 - y[11] + y[12] - 2 y[14] && x[3] == -1 - y[10] + y[11] - 2 y[13] + 2 y[14] - 3 y[15] && x[2] == 1 + y[10] + 2 y[13] + 3 y[15] && x[1] == 0) 

variables:

yvars = Union@Cases[eqs, y[_], Infinity] xvars = Union@Cases[eqs, x[i_], Infinity] 

target:

eqsAndCons = Join[eqs, y[#] >= 0 & /@ Range[15], x[#] >= 0 & /@ Range[15]] 

instance:

fi=FindInstance[eqsAndCons, Join[yvars, xvars], Reals ] (*{{y[1] -> 1, y[2] -> 1, y[3] -> 1, y[4] -> 1, y[5] -> 0, y[6] -> 0, y[7] -> 1, y[8] -> 1, y[9] -> 0, y[10] -> 0, y[11] -> 1, y[12] -> 0, y[13] -> 0, y[14] -> 0, y[15] -> 0, x[1] -> 0, x[2] -> 1, x[3] -> 0, x[4] -> 0, x[5] -> 0, x[6] -> 0, x[7] -> 2, x[8] -> 1, x[9] -> 0, x[10] -> 0, x[11] -> 2, x[12] -> 0, x[13] -> 0, x[14] -> 1, x[15] -> 0}}*) 

Check:

Union[Flatten[eqs /. fi]] (*{True}*) 

Indeed, we can solve it as a linear optimization problem.

First convert the equations to a list, include the non-negativity condition, and call the LinearOptimization function.

cons1=y[9]==0&&y[8]==1-2 y[12]&&y[7]==3-2 y[11]-3 y[14]&&y[6]==-2 y[10]-3 y[13]-4 y[15]&&y[5]==0&&y[4]==1&&y[3]==1+y[12]&&y[2]==y[11]+2 y[14]&&y[1]==1+y[10]+2 y[13]+3 y[15]&&x[15]==0&&x[14]==1+y[15]&&x[13]==-y[15]&&x[12]==y[14]-y[15]&&x[11]==2+y[13]-y[14]+3 y[15]&&x[10]==-y[13]-2 y[15]&&x[9]==y[12]-y[14]&&x[8]==y[11]-y[12]-y[13]+3 y[14]-2 y[15]&&x[7]==3+y[10]-y[11]+3 y[13]-2 y[14]+5 y[15]&&x[6]==-y[10]-2 y[13]-3 y[15]&&x[5]==-y[12]&&x[4]==1-y[11]+y[12]-2 y[14]&&x[3]==-1-y[10]+y[11]-2 y[13]+2 y[14]-3 y[15]&&x[2]==1+y[10]+2 y[13]+3 y[15]&&x[1]==0; cons1=cons1/.And->List; vars=Flatten@Table[{x[i], y[i]}, {i, 1, 15}]; cons2=#>=0&/@vars; cons=Union[cons1,cons2]; LinearOptimization[0,cons,vars] 

The result is

{x[1]->0,y[1]->1,x[2]->1,y[2]->1,x[3]->0,y[3]->1,x[4]->0,y[4]->1,x[5]->0,y[5]->0,x[6]->0,y[6]->0,x[7]->2,y[7]->1,x[8]->1,y[8]->1,x[9]->0,y[9]->0,x[10]->0,y[10]->0,x[11]->2,y[11]->1,x[12]->0,y[12]->0,x[13]->0,y[13]->0,x[14]->1,y[14]->0,x[15]->0,y[15]->0}