The answer for the shortest distance between two skew lines has already been answered here and is also found on Wikipedia.
However, I need to apply a further constraint. I need to find the shortest distance between two skew line segments. Or, if it is easier, two rays, whose starting points are $P_1$ and $P_2$. How can I modify the equations for $c_1$ and $c_2$ in the wiki link to account for this? Or is this a far more complicated problem?
Defining
$$ s_1(\lambda)=\lambda p_1+(1-\lambda)p_2\\ s_2(\mu) = \mu p_3+(1-\mu)p_4 $$
with $0\le \lambda\le 1$ and $0 \le \mu\le 1$ and with the help of four slack variables $\epsilon_1,\epsilon_2,\epsilon_3,\epsilon_4$ to handle the inequality restrictions, the problem can be stated using the Lagrange Multipliers as follows
$$ L(\lambda, \mu, \eta,\epsilon) = \frac{1}{2}||s_1(\lambda)-s_2(\mu)||^2+\eta_1 (\lambda-\epsilon_1^2) +\eta_2(\mu-\epsilon_2^2)+\eta_3(\lambda-1+\epsilon_3^2)+\eta_4(\mu-1+\epsilon_4^2) $$
with the stationary conditions
$$ \nabla L = 0 $$
or
$$ \left\{ \begin{array}{rcl} (p_2+\lambda(p_1-p_2))(p_1-p_2)+\eta_1 + \eta_3 & = & 0\\ (p_4+\mu(p_3-p_4))(p_3-p_4)+\eta_2 + \eta_4 & = & 0\\ \lambda-\epsilon_1^2 & = & 0\\ \mu-\epsilon_2^2 & = & 0\\ \lambda-1+\epsilon_3^2 & = & 0\\ \mu-1+\epsilon_4^2 & = & 0\\ \eta_1\epsilon_1 & = & 0\\ \eta_2\epsilon_2 & = & 0\\ \eta_3\epsilon_3 & = & 0\\ \eta_4\epsilon_4 & = & 0 \end{array} \right. $$
Solving those conditions we get the possible stationary points.