I have about twenty five (multilinear) polynomials $f\_1(\mathbf{x}), f\_2(\mathbf{x}), \dots, f\_{25}(\mathbf{x})$ all in fifteen variables and I would like to decide if there is a $\mathbf{y} \in [0,1]^{15} \subset \mathbb{R}^{15}$ such that $f\_i(\mathbf{y}) \ge 0$ for all $1 \le i \le 25$.


I would like to know what the best way to approach this problem on a computer might be (which is why I specified the numbers 25 and 15 above). I have tried the FindInstance method in Mathematica but no luck there - <https://reference.wolfram.com/language/ref/FindInstance.html>


I am aware of the Positivstellensatz which says that the only obstacle to the existence of a solution is if $-1$ belongs to the "positive cone" of these polynomials (and the unit cube constraints) - <https://en.wikipedia.org/wiki/Stengle%27s_Positivstellensatz>


While looking for $-1$ in the positive cone sounds like an unbounded search, I would be happy if I could at least rule out making $-1$ using polynomials in the positive cone of degree bounded by at most 10, say (or find such a representation). This bounded degree search sounds like a (big) linear programming problem, and rather than reinvent the wheel, I was wondering if there might be a software package out there that already has the capability to do this. 



Solving polynomial inequalities -- efficient Positivstellensatz on a computer
You might try an optimization-based approach: give an optimizer the problem


maximize $z$ subject to $z - f\_i(y) \le 0$, $i = 1 \ldots 25$, and $0 \le y\_j \le 1$, $j = 1 \ldots 15$.


This might be too difficult for a global optimizer (but it might be worth a try).
If it's a local optimizer, you're not guaranteed to find a global maximum, so 
if the solutions found (perhaps with several randomized starting points) all have $z < 0$ you can't conclude that the answer to your question is no; but if you do find a solution with $z \ge 0$ you know the answer is yes. 


