As you said, each line can be written, in parametric form, as
$$
\left( {\matrix{
x \cr
y \cr
x \cr
} } \right) = \left( {\matrix{
{a_x } \cr
{a_y } \cr
{a_z } \cr
} } \right) + \lambda \left( {\matrix{
{v_x } \cr
{v_y } \cr
{v_z } \cr
} } \right)
$$
So we can represent each group of three lines in matricial notation as
$$
{\bf X} = \left( {\matrix{
{x_{\,1} } & {x_{\,2} } & {x_{\,3} } \cr
{y_{\,1} } & {y_{\,2} } & {y_{\,3} } \cr
{z_{\,1} } & {z_{\,2} } & {z_{\,3} } \cr
} } \right) = {\bf A} + {\bf V}\;{\bf \Lambda }\quad \quad {\bf X}' = \left( {\matrix{
{x'_{\,1} } & {x'_{\,2} } & {x'_{\,3} } \cr
{y'_{\,1} } & {y'_{\,2} } & {y'_{\,3} } \cr
{z'_{\,1} } & {z'_{\,2} } & {z'_{\,3} } \cr
} } \right) = {\bf A}' + {\bf V}'\;{\bf \Lambda }'
$$
where is clear the meaning of the matrices , and in particular that the ${\bf \Lambda }$ and ${\bf \Lambda }'$ matrices are diagonal.
Now if each line of the first group was to intersect the corresponding line in the second group, that means
that there are specific values of $\lambda_k$ and ${\lambda '}_j$ such that ${\bf X} ={\bf X}' $, i.e.
$$
{\bf A} + {\bf V}\;{\bf \Lambda } = {\bf A}' + {\bf V}'\;{\bf \Lambda }'
$$
or
$$
\;{\bf \Lambda } = {\bf V}^{\, - \,{\bf 1}} \left( {{\bf A}' - {\bf A}} \right) + {\bf V}^{\, - \,{\bf 1}} {\bf V}'\;{\bf \Lambda }' = diag
$$
That translates into a linear system of six equations in the three unknowns $\lambda'$, which then determine the $\lambda$.
If instead each line of the first group was to intercept one line of the 2nd group, different but not necessarily the corresponding one,
then we shall have ${\bf X} \, {\bf P} ={\bf X}' $, or v.v, where ${\bf P}$ is one of the $6$ Permutation matrices.
Also, if each line of the first group was to intercept one line of the 2nd group, not necessarily distinct,
then we shall have ${\bf X} ={\bf X}' \, {\bf Q} $, where ${\bf Q}$ is one of the $9$ "$0/1$" matrices, having only one $1$ in each column.
Finally, as per your question, we can introduce a Rotation matrix, depending on $3$ angles, which goes to multiply ${\bf X}$
$$
{\bf R}\,\left( {{\bf A} + \,{\bf V}\;{\bf \Lambda }} \right) = \left( {{\bf A}' + {\bf V}'\;{\bf \Lambda }'} \right)\;\left( {{\bf P}\,{\rm or}\;{\bf Q}} \right)
$$
to solve by imposing the diagonalization condition as above. i.e.
$$
{\bf \Lambda } = {\bf V}^{\, - \,{\bf 1}} \;{\bf R}^{\, - \,{\bf 1}} \;\left( {{\bf A}'{\bf Q} - {\bf R}\,{\bf A}} \right) + {\bf V}^{\, - \,{\bf 1}} \;{\bf R}^{\, - \,{\bf 1}} {\bf V}'\;{\bf \Lambda }'{\bf Q}\;
$$
By imposing that the "6" off-diagonal elements of the RHS matrix be null,
we get system of $6$ equations in the $3$ unknowns ${\lambda}'$, and in the $3$ rotation angles,
for each of the $9$ ${\bf Q}$ matrices we are willing to consider as "crossing conditions".
( let comprise in ${\bf Q}$ also the particular cases of unit and permutation matrices).
Now, the homogeneous system is linear in the ${\lambda}'$ but contains product of $\sin $ and $ \cos$ of the angles,
and you are asking for a way to simplify that.
The first way I can see at the moment is that, restarting from the identity
$$
{\bf R}\,\left( {{\bf A} + \,{\bf V}\;{\bf \Lambda }} \right) = \left( {{\bf A}' + {\bf V}'\;{\bf \Lambda }'} \right)\;{\bf Q}
$$
taking the transposed (indicated by the hat)
$$
\left( {\overline {\bf A} + \;{\bf \Lambda }\,\overline {\bf V} } \right){\bf R}^{\, - \,{\bf 1}} = \overline {\bf Q} \left( {\overline {{\bf A}'} + {\bf \Lambda }'\;\overline {{\bf V}'} } \right)\;
$$
we can get rid of the rotation, by multiplying the first by the second on the left
$$
\left( {\overline {\bf A} + \;{\bf \Lambda }\,\overline {\bf V} } \right)\left( {{\bf A} + \,{\bf V}\;{\bf \Lambda }} \right) = \overline {\bf Q} \left( {\overline {{\bf A}'} + {\bf \Lambda }'\;\overline {{\bf V}'} } \right)\left( {{\bf A}' + {\bf V}'\;{\bf \Lambda }'} \right)\;{\bf Q}
$$
however at the cost of introducing an equation between two quadratic forms, symmetric, so $6$ equations.
Yet this equation has the advantage of not requiring $\bf V$ to be invertible.
At this point might be profitable to express each line with the position vector ($\bf a$)
normal to the direction vector ($\bf v$) and to normalize the latter.
That will give that $\overline {\bf A} \, {\bf V}$ has the main diagonal null, while for
$\overline {\bf V} \, {\bf V}$ it is unitary.
So we might possibly start from equating the diagonals ...
Sorry for not being this actually an answer, but this is the best help I can offer.
--- P.S. -----
As per your last comment, the case of two groups of incident lines introduces many simplifications with respect to the general case of being all skew, and can be approached quite differently. It is worthy that you open another post. In any case it remains to clarify the intersection scheme (e.g. two lines 1st group could intersect the same line in 2nd group, and the third any one of the remaining?)
Best Answer
Use the law of cosines. If $\rho_i$ denotes the distance of $p_i$ from the origin, and you want the distance between the centres to be $r_1+r_2$, you need the angle bewteen the lines from the origin to the centres to form an angle $\theta$ that satisfies
$$ (r_1+r_2)^2=\rho_1^2+\rho_2^2-2\rho_1\rho_2\cos\theta\;, $$
and thus
$$ \theta=\arccos\frac{\rho_1^2+\rho_2^2-(r_1+r_2)^2}{2\rho_1\rho_2}\;. $$
Then you just need to calculate the angle originally formed by these lines and rotate by the difference of the two angles.