As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
Best Answer
Rewrite your first hyperboloid equation as $$ \sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2} - d_{12}= \sqrt{(x-x_2)^2+(y-y_2)^2+(z-z_2)^2}. $$ Square both sides and rearrange to get $$ \tag{1} 2x(x_2-x_1)+2y(y_2-y_1)+2z(z_2-z_1) +d_1^2 - d_2^2 + d_{12}^2= 2d_{12}\sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2}, $$ where $d_1^2=x_1^2+y_1^2+z_1^2$ and $d_2^2=x_2^2+y_2^2+z_2^2$. Do the same for the other hyperboloid: $$ \tag{2} 2x(x_3-x_1)+2y(y_3-y_1)+2z(z_3-z_1) +d_1^2 - d_3^2 + d_{13}^2= 2d_{13}\sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2}, $$ where $d_3^2=x_3^2+y_3^2+z_3^2$. If you now multiply $(1)$ by $d_{13}$, $(2)$ by $d_{12}$, and subtract $(2)$ from $(1)$, you end up with a linear equation: $$ \tag{3} 2x[d_{13}(x_2-x_1)-d_{12}(x_3-x_1)]+ 2y[d_{13}(y_2-y_1)-d_{12}(y_3-y_1)]+ 2z[d_{13}(z_2-z_1)-d_{12}(z_3-z_1)] +d_{13}(d_1^2 - d_2^2 + d_{12}^2) -d_{12}(d_1^2 - d_3^2 + d_{13}^2)= 0. $$ This is the equation of a plane: the locus you need is then the intersection between this plane and either hyperboloid, that is a conic section (hyperbola or ellipse)
As you can see the equations are getting more and more complicated, so I will just sketch a possible strategy to parameterize that conic.
Square equation $(1)$ (or $(2)$, it doesn't matter) and simplify the result.
Insert into that equation the expression for $z$ obtained from $(3)$ (if equation $(3)$ doesn't contain $z$ just solve for $x$ or $y$ instead).
The equation you thus get is that of a conic section in plane $xy$, projection onto that plane of the original conic. You can find its center and axes with the usual formulas and thus find parametric equations $x(t)$ and $y(t)$ for its coordinates.
Plug $x(t)$ and $y(t)$ into the expression for $z$ obtained above, to get also z(t).