NSolve and domain specifications

I am using NSolve to find roots to a system of polynomials that describe some chemical reactions. However, I am seeing some odd behavior when I restrict the domain.

First, the basic setup:

c1 = 50.0;
c2 = 1.5;
c3 = 0.4;
c4 = 0.5;
c5 = 8.0;
c6 = 3.2;
k1 = 3.7*10^-18;
k2 = 8.3*10^-9;
k3 = 8.7*10^-13;
k4 = 3.2*10^-7;
k5 = 8.6*10^-3;
f1 = 1/2*(c1-c2+c3+2*c4-c5+c6-x1-2*x2-x3+x4+(x3*x4)/(k5*x1));
f2 = 1/(2*k5*x1)*(3*k5*c6*x1-k5*c1*x1+2*k5*c4*x1+k5*c3*x1-k5*c2*x1-k5*c5*x1+k5*x1^2-2*k5*x1*x2-k5*x1*x3+k5*x1*x4-x3*x4);
eqs={
f1*f2 - k1*x1^2 == 0,
f1*x3 - k2*x2*x1 == 0,
f1*(c6 - x3 - x2 - (x3*x4)/(k5*x1)) - k3*x1*x3 == 0,
f2*(c5 - x4 - (x3*x4)/(k5*x1)) - k4*x1*x4 == 0
};
NSolve[eqs,{x1, x2, x3, x4}]

This returns 32 complex roots. However, I know at the onset that there is at least one real-valued root, and exactly one real root with all four variables greater than zero. So, I try it again and use the specification suggested in the documentation to search for Reals only.

NSolve[eqs, {x1, x2, x3, x4}, Reals]

This, however, returns {} for no solutions. However, while playing around with the options I accidentally used "Real" instead of "Reals":

NSolve[eqs, {x1, x2, x3, x4}, Real]

Now, I receive a warning:

NSolve::bdomv: Warning: Real is not a valid domain specification. Assuming it is a variable to eliminate. >>

But I also get 8 real roots! Interestingly, these are the same roots I find if I run the original calculation at any WorkingPrecision that is not MachinePrecision...

NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->MachinePrecision]
NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->2]
NSolve[eqs, {x1, x2, x3, x4}, WorkingPrecision->20]

Lastly, I also tried restricting the domain to positive reals. This fails to return any result. If I combine this with a finite WorkingPrecision I get the correct root. If I instead specify "Reals" I get nothing, and I use "Real" I get warnings about the domain as well as infinite solutions, but the same root as above when I used Working Precision.

NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, WorkingPrecision -> 8]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, Reals]
NSolve[Catenate[{eqs, {x1 > 0, x2 > 0, x3 > 0, x4 > 0}}], {x1, x2, x3, x4}, Real]

{}
{{x1 -> 49.4605, x2 -> 0.00100605, x3 -> 0.511294, x4 -> 2.21332}, {x1 -> 53.0999, x2 -> 4.30992*10^-13, x3 -> 0.000114669, x4 -> 7.99797}}
{}
NSolve::bdomv: Warning: Real is not a valid domain specification. Assuming it is a variable to eliminate. >>
NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with -((41688 x1)/65167)-(153968 x2)/195501+(153196 x3)/195501+(185938 x4)/195501+(38650 SystemNSolveDumpY\$31197)/65167 == 1. >>
{{x1 -> 49.4605, x2 -> 0.00100591, x3 -> 0.511294, x4 -> 2.21332}}

Anyone have insight into these behaviors?

Seems to be a case of ill conditioning of the input system. If I redo using exact input and set NSolve to work on high precision then I get a plausible outcome.

c1 = 50;
c2 = 3/2;
c3 = 2/5;
c4 = 1/2;
c5 = 8;
c6 = 16/5;
k1 = 37*10^(-17);
k2 = 83*10^(-8);
k3 = 87*10^(-12);
k4 = 32*10^(-6);
k5 = 86*10^(-32);
eqs = {f1*f2 - k1*x1^2, f1*x3 - k2*x2*x1,
f1*(c6 - x3 - x2 - (x3*x4)/(k5*x1)) - k3*x1*x3,
f2*(c5 - x4 - (x3*x4)/(k5*x1)) - k4*x1*x4};
solns = N[
NSolve[eqs, {x1, x2, x3, x4}, Method -> "EndomorphismMatrix",
WorkingPrecision -> 200]]

(* Out= {{x1 -> 54.8004248857, x2 -> -0.0000977553043262,
x3 -> -1.70022938028,
x4 -> -2.2175062235*10^-28}, {x1 -> 53.0993625312,
x2 -> -8.87278874951*10^-12, x3 -> 0.000637468817923,
x4 -> 5.73084681075*10^-25}, {x1 -> 51.7000000017,
x2 -> 3.33952465735*10^-15, x3 -> -8.42954955495*10^-20,
x4 -> -1.68785294015*10^-9}, {x1 -> 46.9000236518,
x2 -> -4.79997634819, x3 -> -0.0000236518142965,
x4 -> -1.364259666*10^-23}, {x1 -> 53.0973683194,
x2 -> 2.32767391416*10^-37, x3 -> 3.04591909329*10^-29,
x4 -> 4.79736831936}, {x1 -> 50.0008781357,
x2 -> 2.05348705909*10^-33, x3 -> 8.0900808571*10^-29,
x4 -> 1.70087813781}, {x1 -> 54.8000177886, x2 -> 3.10212897782,
x3 -> -7.90424004255,
x4 -> -4.76989717362*10^-29}, {x1 -> 2.09855075596*10^-11,
x2 -> -51.7, x3 -> -8.85456356826*10^-18,
x4 -> -1.11898203258*10^-22}} *)

Check residuals:

eqs /. solns

(* Out= {{1.45782281955*10^-18, -5.83359947449*10^-15, \
-1.06351826451*10^-14, -7.54753185026*10^-19}, {1.25900746755*10^-15, \
-4.71927742874*10^-19,
3.55398404018*10^-15, -9.73773799727*10^-28}, \
{-7.83975673111*10^-15, 9.62964972194*10^-35, 3.79152709445*10^-28,
2.21357836972*10^-14}, {-5.8353686176*10^-15,
0., -1.72569257779*10^-15,
1.25720897678*10^-28}, {4.9846880603*10^-15,
4.90188927228*10^-44, -1.40705350532*10^-37, \
-1.61849700309*10^-15}, {6.95980079332*10^-19, 6.41188041986*10^-44,
4.67805560848*10^-25,
9.91090889912*10^-15}, {9.35073564361*10^-20, -2.37330771355*10^-14,
6.3636526536*10^-18,
1.10569373992*10^-22}, {4.59635768258*10^-13, -1.97215226305*10^-31,
1.61661234308*10^-38, -2.11965757437*10^-13}} *)

A couple of them are all positive.

Select[{x1, x2, x3, x4} /. solns,
Apply[And, Thread[# >= 0]] &]

(* Out= {{53.0973683194, 2.32767391416*10^-37, 3.04591909329*10^-29,
4.79736831936}, {50.0008781357, 2.05348705909*10^-33,
8.0900808571*10^-29, 1.70087813781}} *)

There may be other ways to do this. I had but little luck using FindMinimum on the sum of squares but there might be option settings that would make that work reasonably well.

This reminds me of NSolve finds real-valued results in version 9, but not in version 10, in which you can use any of the Method settings from Methods for NSolve

"EndomorphismMatrix"
"CompanionMatrix"
"Legacy"
"Aberth"
"JenkinsTraub"

or even a nonexistent method "Foo":

NSolve[eqs, {x1, x2, x3, x4}, Reals, Method -> "Foo"]
(*
{{x1 -> 23.2915, x2 -> -28.505, x3 -> 55.2247, x4 -> -0.0852683},
{x1 -> 46.3516, x2 -> -5.36977, x3 -> 9.01972, x4 -> -0.0198153},
{x1 -> 56.3754, x2 -> -0.000327403, x3 -> -3.16593, x4 -> -1.46634},
{x1 -> 54.9, x2 -> 3.2, x3 -> -1.80018*10^-7, x4 -> -6.26944*10^-10},
{x1 -> 52.6254, x2 -> 0.0121417, x3 -> -2.70055, x4 -> -0.976273},
{x1 -> 52.2298, x2 -> -0.0000285978, x3 -> 0.928352, x4 -> 2.62754},
{x1 -> 49.4605, x2 -> 0.00100605, x3 -> 0.511294, x4 -> 2.21332},
{x1 -> 53.0999, x2 -> 4.30992*10^-13, x3 -> 0.000114669, x4 -> 7.99797}}
*)