Following @littleO's advice, I've set about to minimize $\sum_n ((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-d^2)^2$. Going using an exact Hessian (because the function is smooth definite) as follows:
$\textbf{H} = 4 * \begin{bmatrix} \sum_n [2(x-x_n)^2 + ((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-d_n^2)] && \sum_n 2(x-x_n)(y-y_n) && \sum_n 2(x-x_n)(z-z_n) \\ \sum_n 2(y-y_n)(x-x_n) && \sum_n [2(y-y_n)^2 + ((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-d_n^2)] && \sum_n 2(y-y_n)(z-z_n) \\ \sum_n 2(x-x_n)(z-z_n) && \sum_n 2(y-y_n)(z-z_n) && \sum_n [2(z-z_n)^2 + ((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-d_n^2)] \end{bmatrix}$
and gradient:
$\nabla F = 4 \begin{bmatrix} \sum_n (x-x_n)((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-{d_n}^2) \\ \sum_n (z-z_n)((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-{d_n}^2) \\ \sum_n (z-z_n)((x-x_n)^2+(y-y_n)^2+(z-z_n)^2-{d_n}^2) \end{bmatrix}$
I'm still taking the barycenter of points $(x_n, y_n, z_n)$ as initial guess, and calculating the next step with $\textbf{x}_{n+1} = \textbf{x}_n - H^{-1} \nabla F$.
The current code, which should run as-is, can be found on this gist.
I'm pretty sure something still went wrong, however. While building this, I added a maximum deviation from expected values (which is really $\max{ f_n(\textbf{x}_{guess})}$). I expect this to decrease as I iterate closer and closer to a minimum. However, this is not the case. Not only am I see seriously weird iterations (convergence within 5 iterations to a value that is 4200km away from one base station - theoretically impossible), but the converging value itself makes little to no sense.
I get the impression I have made a mistake somewhere.
Just recording my comment here so that this question receives an answer:
The suggestion was to minimize \begin{equation} \sum_n ((x - x_n)^2 + (y - y_n)^2 + (z - z_n)^2 - d_n^2)^2 \end{equation} using Newton's method.