My problem setup is as shown below:
I know the location (x,y) of fixed points p1+, p1-, p2+, and p2-, and want to find the position marked "o" by trilateration. However, I do not know the individual distances r1+, r1-, r2+ and r2-. Rather, I have the following system of equations available:
$\begin{align} \frac{1}{r_{1+}}-\frac{1}{r_{1-}}&=V_1\\ \frac{1}{r_{2+}}-\frac{1}{r_{2-}}&=V_2\\ \frac{1}{r_{3+}}-\frac{1}{r_{3-}}&=V_3\\ \frac{1}{r_{4+}}-\frac{1}{r_{4-}}&=V_4\\ \frac{1}{r_{5+}}-\frac{1}{r_{5-}}&=V_5 \end{align}$
I couldn't figure out an analytical solution to this system of equations, even though there are really only 3 unknowns -- (x,y) of the the point marked o, and thought the solution needs to be found numerically...
Thus my question is, is my conclusion correct in that analytical solution isn't possible?
Edit: 11/15/2018
So after concluding analytical solution is most likely not feasible, I decided to use Newton's method:
$\mathbf{x}^{(k)}=\mathbf{x}^{(k+1)}-\mathbf{J}(\mathbf{x}^{(k+1)})^{-1}\mathbf{F}(\mathbf{x}^{(k-1)})$
where $\mathbf{x}^{(k)}$ is the estimate of the unknown point's position $[x^{(k)};y^{(k)};z^{(k)}]$ on the k-th iteration,
$\mathbf{F}$ is the value of my equations, where the i-th row is:
$V_i-\frac{1}{\sqrt{(x^{(k)}-x_{i+})^2+(y^{(k)}-y_{i+})^2+(z^{(k)}-z_{i+})^2}}+\frac{1}{\sqrt{(x^{(k)}-x_{i-})^2+(y^{(k)}-y_{i-})^2+(z^{(k)}-z_{i-})^2}}$,
where $[x_{i+},y_{i+},z_{i+}]$ and $[x_{i-}, y_{i-},z_{i-}]$ are the position of known points pi+ and pi- as shown in the picture.
$\mathbf{J}$ is the Jacobian matrix.
In the Newton update equation, getting $\mathbf{J}(\mathbf{x})^{-1}$ requires the Jacobian to be square. In my problem, there are more equations (5) than variables (3). Using the pseudo-inverse such that the update equation becomes
$\mathbf{x}^{(k)}=\mathbf{x}^{(k+1)}-(\mathbf{J}^T \mathbf{J})^{-1}\mathbf{J}^T\mathbf{F}(\mathbf{x}^{(k-1)})$
gives unstable behavior and yielded estimates of the unknown point falling outside the circle.
I also tried using a "boosting method" -- at each iteration k, pick three rows where F is the greatest...this also had convergence issues.
In the end, I decided on using gradient descent, with loss function being the sum-of-squared differences between estimated and actual $V$ values.

I realize now that I haven't taken into account the fact that the values $r1+, r1-$, etc. are only known through the system of equations you give. But nevertheless, I am going to show that even if we know these four values, we are practically unable to find an exact solution.
I recall first that an ellipse with focii $F$ and $F'$ is the locus of points $M$ such that the sum of distances $MF+MF'$ is a constant $k$.
Here, point $O$ belongs in particular to
the ellipse with focii $p1-$ and $p2-$ and constant $k=r1-+r2-$, and to
the ellipse with focii $p1+$ and $p2+$ and constant $k=r1++r2+$.
Otherwise said, point $O$ is an intersection point of these two ellipses. Knowing that an ellipse has a second degree equation, the coordinates of $O$ involve a $2^2=4$th degree equation (a rapid check : indeed, two ellipses can have up to 4 intersection points).
Thus, theoretically, you could get an (awfully complicated) analytical answer, but practically speaking, you have to look for a numerical method.
Remarks :
1) One could have chosen other ellipses such as the one with focii $p2+$ and $p1-$. There is no guarantee that the $\binom{4}{2}=6$ ellipses and the $\binom{6}{2}=15$ intersection points (at least) you can get in this way will intersect in a very same point due to measuring errors. Therefore, the numerical method you will use needs a "least squares" step.
2) one could have used hyperbolas instead of ellipses.
3) the fact that the 4 points $p1+, ...$ belong to a same circle hasn't been used.