Consider system of non-linear equations $$ \begin{align} \ln \frac{x_1}{1 - x_1 - \ldots - x_n} + P_1(x_1, \ldots, x_n) &= 0 \\ &\ldots \\ \ln \frac{x_n}{1 - x_1 - \ldots - x_n} + P_n(x_1, \ldots, x_n) &= 0 \end{align} $$ with $P_i(x_1, \ldots, x_n)$ being polynomials with coefficients in $\mathbb{R}$ and $x_1, \ldots, x_n$ confined to the interior of the standard real symplex $x_1, \ldots, x_n \in \mathbb{R}$, $x_1, \ldots, x_n > 0$, $x_1 + \ldots + x_n < 1$.

Is there a general method for solving such a system analytically and/or numerically?

Can we compute its Groebner basis using standard algorithms for systems of polynomial equations and determine if it's positive- or zero-dimentional? In what ring we would compute it? $\mathbb{R}[x_1, \ldots, x_n, l_1, \ldots, l_n]$ with $$ l_i = \ln \frac{x_i}{1 - x_1 - \ldots - x_n}? $$

If the system proves to be zero-dimentional, can we somehow bound its number of solutions from above (perhaps using some generalization of Bézout's theorem?), so that we may be able to solve it numerically via homotopy continuation starting from some system of polynomial equations with suitable number of solutions?