Two spheres intersecate (if they do) in a circle that lies on a plane
that is normal to the line joining the centers, and that crosses that line
at a distance $d_{1,2}$ from point $C1$, that can be found by solving
the triangle with the given sides, as shown
that is
$$
d_{\,1,2} = {{r_{\,1} ^{\,2} - r_{\,2} ^{\,2} + \left| {{\bf v}_{\,12} } \right|^{\,2} } \over {2\left| {{\bf v}_{\,12} } \right|}}
$$
which can be either positive or negative.
The equation of the plane will then be:
$$
\left( {{\bf x} - {\bf c}_{\,1} } \right) \cdot {{{\bf v}_{\,12} } \over {\left| {{\bf v}_{\,12} } \right|}} = d_{\,1,2}
$$
which is the same as to write
$$
\eqalign{
& \left( {x - C_{\,1,\,x} } \right) \cdot \left( {C_{\,2,\,x} - C_{\,1,\,x} } \right) + \left( {y - C_{\,1,\,y} } \right) \cdot \left( {C_{\,2,\,y} - C_{\,1,\,y} } \right) + \left( {z - C_{\,1,\,z} } \right) \cdot \left( {C_{\,2,\,z} - C_{\,1,\,z} } \right) = \cr
& = {{r_{\,1} ^{\,2} - r_{\,2} ^{\,2} + \left| {{\bf v}_{\,12} } \right|^{\,2} } \over 2} \cr}
$$
Do the same for other two couples of points, choosen in such a way that the
three vectors $ {{\bf v}_{\,jk} }$ be oriented as much differently as possible.
Then find the point where the three planes cross, i.e. the solution to the system
of the three linear equations they define.
If the 4 spheres define a unique point, then it must be that, apart of course from
measurement errors/approximations.
Check the distances of the cross point from the four references.
In case of big discrepances, find some other planes and find a common solution
by the method of least squares, again checking vs. the known distances.
Note
If you are writing a "real applicable" program, the error checking and handling is fundamental. Then you need to know probability theory, statistics etc. (and GPS engineering, of course) to decide which method strategy to adopt.
One can always construct two parabolas passing through $p_1$, $p_2$, $p_3$, $p_4$ (green and pink in the figure below), each one possibly degenerating into a couple of parallel lines if two opposite sides of quadrilateral $p_1p_2p_3p_4$ are parallel. Point $p_5$ will determine an ellipse if it lies inside either parabola but not in their intersection.
This follows from the fact that five points always determine a conic section, and because the parabola is a limiting case between ellipse and hyperbola: each time $p_5$ crosses the boundary of a parabola, conic section $p_1p_2p_3p_4p_5$ switches from ellipse to hyperbola (or viceversa).
Best Answer
Here's a proof without coordinates. I will relabel $D_1=P_1P_4\cap P_2P_3$, $D_2=P_2P_4\cap P_1P_3$ and $D_3=P_3P_4\cap P_1P_2$.
Lemma 1. The centre of any rectangular hyperbola $\mathcal H$ through points $A$, $B$, $C$ lies on the nine-point circle of $\triangle ABC$.
Proof. Let $(ABC)$ meet $\mathcal H$ again at $D$, and let $H_A$, $H_B$, $H_C$, $H_D$ be the orthocentres of $\triangle BCD$, $\triangle CDA$, $\triangle DAB$, $\triangle ABC$, which all lie on $\mathcal H$. There is a $180^{\circ}$ rotation mapping $ABCD$ to $H_AH_BH_CH_D$, and the centre of this rotation is both the centre of $\mathcal H$ and the midpoint of $\overline{DH_D}$, so it lies on the nine-point circle of $\triangle ABC$. $\square$
Lemma 2. The incentre $I$ and excentres $I_1, I_2, I_3$ (opposite $D_1$, $D_2$, $D_3$) of $\triangle D_1D_2D_3$ lie on the hyperbola.
Proof. Take a projective transformation sending $P_1P_2P_3P_4$ to a square. Then $D_1$ and $D_3$ go to infinity, so $I_1I_2I_3I$ becomes a rectangle with sides parallel to $P_1P_2P_3P_4$. Then $D_2$ maps to the common centre of rectangle $I_1I_2I_3I$ and square $P_1P_2P_3P_4$, so these eight points lie on a common conic. $\square$
Now $(D_1D_2D_3)$ is the nine-point circle of $\triangle I_1I_2I_3$, so our two lemmas imply the result.