First the equations for all vectors $x$ on line $g$ and all vectors $y$ on line $h$:
$$
\begin{align}
g: x &= a + \lambda b \quad \\
h: y &= c + \mu d
\end{align}
\quad (*)
$$
The difference vector between two of those vectors is
$$
D = y - x = c - a + \mu d - \lambda b
$$
The length of $D$ squared is:
\begin{align}
q(\lambda, \mu)
&= D \cdot D \\
&= (c - a)^2 + (\mu d-\lambda b)^2 + 2 (c-a)\cdot(\mu d - \lambda b) \\
&= (c - a)^2 + \mu^2 d^2 + \lambda^2 b^2 - 2 \mu \lambda (d\cdot b)
+ 2 \mu ((c-a)\cdot d) - 2 \lambda ((c-a)\cdot b) \\
\end{align}
The gradient of $q$ is:
$$
q_\lambda = 2\lambda b^2 - 2\mu (d\cdot b)-2((c-a)\cdot b) \\
q_\mu = 2\mu d^2 - 2\lambda (d\cdot b)+2((c-a)\cdot d)
$$
It should vanish for local extrema of $q$ which leads to the system
$$
A u = v
$$
with
$$
A =
\left(
\begin{matrix}
b^2 & - d\cdot b \\
- d \cdot b & d^2
\end{matrix}
\right)
\quad
u =
\left(
\begin{matrix}
\lambda \\
\mu
\end{matrix}
\right)
\quad
v =
\left(
\begin{matrix}
(c-a) \cdot b \\
-(c-a) \cdot d
\end{matrix}
\right)
$$
A unique solution exists for
$$
0 \ne \mbox{det A} = b^2 d^2 - (d \cdot b)^2
$$
Note that the dot operator stands for the scalar product.
That solution is
$$
u = A^{-1} v
$$
with
$$
A^{-1} = \frac{1}{\mbox{det } A}
\left(
\begin{matrix}
d^2 & d\cdot b \\
d \cdot b & b^2
\end{matrix}
\right)
$$
Inserting the found values for $\lambda$ and $\mu$ into equations $(*)$ will provide you the two vectors, whose points are closest to each other.
Example:
$$
a = (2,0,0), b=(1,1,1), c=(0,1,-1), d=(-1,0,-1)
$$
leads to
$$
\lambda = 1, \mu = -2.5, x_\min = (3,1,1), y_\min=(2.5,1,1.5), D=(-0.5,0,0.5)
$$
The image renders the $x$ values in green, the $y$ values in purple, and $D$ in red.
Since $AH||BG$ and $HF||DB$, we obtain: $(AHF)||(DBG)$.
Id est, our distance it's just the distance between planes: $(AHF)$ and $(DBG).$
But $\overrightarrow{EC}$ it's a normal to both planes.
Thus, the needed distance is equal to $$\frac{1}{3}EC=\frac{s\sqrt3}{3}=\frac{s}{\sqrt3}.$$
Best Answer
Here's an attempt at a purely geometric approach.
Label the two lines $\ell_1$ and $\ell_2$. Select any point $A$ on line $\ell_1$. Construct line $\ell_3$ through $A$ parallel to $\ell_2$. Then the plane $\pi_1$ containing the lines $\ell_1$ and $\ell_3$ is parallel to the line $\ell_2$.
Now find the perpendicular projection $\ell_2'$ of the line $\ell_2$ on the plane $\pi_1$. (One way to do this is to pick any two distinct points $M$ and $N$ on $\ell_2$, find the points $M'$ and $N'$ in plane $\pi_1$ closest to $M$ and $N$, respectively, and construct the line $\ell_2'$ through $M'$ and $N'$.) Let $P$ be the point where the lines $\ell_1$ and $\ell_2'$ intersect.
By similar methods, find the plane $\pi_2$ through line $\ell_2$ parallel to line $\ell_1$, find the projection $\ell_1'$ of the line $\ell_1$ on the plane $\pi_2$, and let $Q$ be the intersection point of the lines $\ell_1'$ and $\ell_2$.
Now planes $\pi_1$ and $\pi_2$ are parallel, lines $\ell_1$ and $\ell_1'$ are perpendicular projections of each other on planes $\pi_1$ and $\pi_2$, and lines $\ell_2$ and $\ell_2'$ are perpendicular projections of each other on planes $\pi_2$ and $\pi_1$, respectively.
In particular, $P$ (the intersection of lines $\ell_1$ and $\ell_2'$) is the perpendicular projection of $Q$ (the intersection of lines $\ell_1'$ and $\ell_2$) on plane $\pi_1$, and likewise $Q$ is the perpendicular projection of $P$ on plane $\pi_2$. The line $PQ$ is perpendicular to both planes $\pi_1$ and $\pi_2$ and to the lines $\ell_1$ and $\ell_2$ in those planes; that is, $PQ$ is the line that was to be found.