It is well known, that given two skew lines \begin{align*} g: x &= a + \lambda b \quad \\ h: y &= c + \mu d \end{align*} the distance is given by $$D = \frac{|(a - c) \cdot b \times d|}{|b \times d|}$$ See e.g. here
I thought, it would be a nice exercise to derive the distance using an optimization method, but I get stuck. Short version: Why is \begin{equation*} \begin{aligned} \left\| \left(a - c\right) - \frac{|b|^2\cdot |d|^2}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right. & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{b}{|b|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{d}{|d|}\right)\right) \cdot \frac{b}{|b|} \right.\right. + \\ & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{d}{|d|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{b}{|b|}\right)\right) \cdot \frac{d}{|d|} \right] \right\| = \frac{|(a - c) \cdot b \times d|}{|b \times d|} \end{aligned} \end{equation*}
It seems clear to me, that if you subtract the orthogonal projections of $a - c$ in direction of $b$ and $d$, the result must be proportional to $b \times d$. On the other hand, there must be a clear line of arguments, why $\frac{|b|^2\cdot |d|^2}{|b \times d|^2}$ exactly does the job.
Long version. For notation and reference see e.g. here: The difference vector between two of those skew lines is $$D = |y - x| = |c - a + \mu d - \lambda b|$$ and the local extrema $(\lambda, \mu)$ is found to be at $$\frac{1}{\mbox{det } A} \left(\begin{matrix} d^2 & d\cdot b \\ d \cdot b & b^2 \end{matrix}\right) \cdot \left(\begin{matrix} (c-a) \cdot b \\ -(c-a) \cdot d \end{matrix}\right)$$
where $$\mbox{det A} = b^2 d^2 - (d \cdot b)^2 = |b \times d|^2$$
If you insert above solution into $D$, you get
\begin{equation*} \begin{aligned} D = |c - a + \mu d - \lambda b| \end{aligned} \end{equation*} \begin{equation*} \begin{aligned} D = \left\| \left(a - c\right) - \left[\vphantom{\frac{b}{b}}\right.\right. & \frac{|d|^2}{|b \times d|^2} \left.\left. \left(\left(a - c\right)\cdot \left(b - \frac{b \cdot d}{|b|\cdot|d|^2} \cdot d\right)\right) \cdot b \right.\right.+\\ & \frac{|b|^2}{|b \times d|^2} \left.\left. \left(\left(a - c\right)\cdot \left(d - \frac{b \cdot d}{|b|^2\cdot|d|} \cdot b\right)\right) \cdot d \right] \right\| \end{aligned} \end{equation*} \begin{equation*} \begin{aligned} D = \left\| \left(a - c\right) - \frac{|b|^2\cdot |d|^2}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right. & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{b}{|b|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{d}{|d|}\right)\right) \cdot \frac{b}{b} \right.\right.+\\ & \left.\left. \left(\left(a - c\right)\cdot \left(\frac{d}{|d|} - \frac{b \cdot d}{|b|\cdot|d|} \cdot \frac{b}{|b|}\right)\right) \cdot \frac{d}{d} \right] \right\| \end{aligned} \end{equation*}
Then I don't see how to proceed...
As pointed out by @ancientmathematician, we could set $u=a-c$, and assume $|b|=|d|=1$ without loss of generality, then \begin{equation*} \begin{aligned} D = \left\| u - \frac{1}{|b \times d|^2} \left[\vphantom{\frac{b}{b}}\right.\right. & \left.\left. \left(u\cdot \left(b - (b \cdot d) \cdot d\right)\right) \cdot b \right.\right.+\\ & \left.\left. \left(u\cdot \left(d - (b \cdot d) \cdot b\right)\right) \cdot d \vphantom{\frac{b}{b}} \right] \right\| \end{aligned} \end{equation*}
and $|b \times d|^{-2}= \sin(\alpha)^{-2}$. But the occurence of the angle between $b$ and $d$ at that place is quite puzzling to me. Do I miss something obvious?
Initially I thought \begin{equation*} \begin{aligned} D = \left\| u - \left[\right.\right. \left.\left. \left(u\cdot \left(b - (b \cdot d) \cdot d\right)\right) \cdot b \right.\right.+ \left.\left. \left(u\cdot \left(d - (b \cdot d) \cdot b\right)\right) \cdot d \right] \right\| \end{aligned} \end{equation*} which is wrong...
Without loss of generality, let $\ c=0,\ $ but we do not let $\ |b|=|d|=1\ $ because we would lose the valuable homogeneous property. We recall the standard scalar Triple product identity:
$$ \Delta := (a \cdot (b \times d))^2 = \begin{vmatrix} a\cdot a & a \cdot b & a\cdot d \\ b \cdot a & b\cdot b & b\cdot d \\ d \cdot a & d\cdot b & d\cdot d \end{vmatrix}. \tag{1} $$
$$ \textrm{Let }\quad b_v := \frac{b}{|b|} \!-\! \frac{b\cdot d}{|b|\,|d|} \frac{d}{|d|}, \quad d_v := \frac{d}{|d|} \!-\! \frac{b\cdot d}{|b|\,|d|} \frac{b}{|b|}, \quad \textrm{ and} \tag{2} $$ $$ a_v := a\ |b\times d|^2 - |b|^2|d|^2 \Big( a\cdot b_v\frac{b}{|b|} + a\cdot d_v\frac{d}{|d|}\Big), \quad D := a_v/|b\times d|^2. \tag{3} $$ A bit of simplification of equation $(2)$ gives us $$ |b|\ |d|^2\ b_v = b\ |d|^2 - d\ (b\cdot d), \quad |b|^2\ |d|\ d_v = d\ |b|^2 - b\ (b\cdot d). \tag{4} $$
Substituting equation $(3)$ in equation $(4)$ gives us
$$ a_v = (a_0)a + (b_0)b + (d_0)d \quad \textrm{ with } \tag{5} $$
$$ a_0 :=|b\!\times\!d|^2,\quad b_0 := -a\!\cdot\!\big(b\ |d|^2 \!-\! d\ (b\!\cdot\! d)\big),\quad d_0 := -a\!\cdot\!\big(d\ |b|^2 \!-\! b\ (b\!\cdot\! d)\big). \tag{6} $$
Taking the dot product of equation $(5)$ we get
$$ a_v\!\cdot\!a_v = (a\!\cdot\!a)a_0^2 \!+\! (b\!\cdot\!b)b_0^2 \!+\! (d\!\cdot\!d)d_0^2 \!+\! 2(a\!\cdot\!b)a_0b_0 \!+\! 2(a\!\cdot\!d)a_0d_0 \!+\! 2(b\!\cdot\!d)b_0d_0. \tag{7}$$
Using $\ |b\!\times\!d|^2 = |b|^2|d|^2\!-\!(b\!\cdot\!d)^2,\ $ and equations $(1)$ and $(6)$ and factoring equation $(7)$ we get
$$ |a_v|^2 = a_v\cdot a_v = |b\times d|^2 \Delta = |b\times d|^2 |a\cdot (b\times d)|^2 \tag{8}$$
Solving for $\ |a_v|\ $ and using equation $(3)$ we get
$$ |a_v| = |b\times d|\ |a\cdot b\times d|,\quad |D| = |a\cdot b\times d|/|b\times d| \tag{9} $$ which is what we wanted to prove.