Distance between two line segments in 3-space

You can easily see that if the nearest approach of the two infinite lines does not fall within both segments then you must calculate the nearest distance of all four end points to the other segment.

 pointsegdis[{seg_, pointlist_}] := 
       Module[{u = Subtract @@ seg, mean = Plus @@ seg/2},
       {#, (mean - u  Sign@# Min[1, Abs@#]/2) &@(-(2 ( # - mean).u )/u.u) } & 
        /@ pointlist ]
 segsegdis[s1_, s2_] := Module[{u, v, w, d, t1, t2},
   {u, v, w} = Subtract @@ # & /@ {Sequence @@ #, First /@ #} &@{s1, s2};
   If[(d = u.u v.v - (v.u)^2) != 0 &&
       Abs[t2 = 2 (u.u v.w - u.w v.u)/d + 1] <= 1 &&
       Abs[t1 = 2 (v.u v.w - v.v u.w)/d + 1] <= 1,
    Plus @@ (#[[1]] (1 + {1, -1} #[[2]] )) /2 & /@ {{s1, t1}, {s2, t2}} ,
    First@SortBy[Flatten[pointsegdis /@ {{s1, s2}, {s2, s1}}, 1], 
                       Total@((Subtract @@ #)^2) &]
       ] ] ;

 {s1, s2} = RandomReal[{-1, 1}, {2, 2, 3}];
 connect = segsegdis[s1, s2];
 Graphics3D[{Line /@ {s1, s2}, {Red, Line[connect]}}]

This is several orders of magnitude faster than Minimise or FindMinimum

 testcases = RandomReal[{-1, 1}, {100, 2, 2, 3}];
 (aa = Norm[ Subtract @@ segsegdis[#[[1]], #[[2]]] ] & /@ testcases) //Timing // First
 (bb = First@
      Minimize[{  Norm[(#[[1, 1]] (1 - t1)/2 + #[[1, 2]] (1 + t1)/2) -
         (#[[2, 1]] (1 - t2)/2 + #[[2, 2]] (1 + t2)/2) ]  , 
            {-1<=t1 <= 1 , -1 <= t2 <= 1}}, {t1, t2}] & /@ testcases) // Timing //First
 (cc = First@
       FindMinimum[{  
        Norm[(#[[1, 1]] (1 - t1)/2 + #[[1, 2]] (1 + t1)/2) -
         (#[[2, 1]] (1 - t2)/2 + #[[2, 2]] (1 + t2)/2) ]  , 
               {-1 <= t1 <= 1 , -1 <= t2 <= 1}}, {t1, t2}] & /@ testcases) // Timing // First

{0.015600, 12.620481, 2.714417}

verify all three give the same result to reasonable numerical precision

 Max[Abs[aa - bb]] -> 5.82532*10^-10
 Max[Abs[aa - cc]] -> 2.25578*10^-6

The link you gave shows how to find the distance between any two points on the lines. You can then use Mathematica's Minimize function to find the shortest distance. Like this.

minDist[{p1_, p2_}, {q1_, q2_}] :=
  Module[{P, Q, u, v, w},
         u = Normalize[p2 - p1];
         v = Normalize[q2 - q1];
         P[s_] := p1 + s u;
         Q[t_] := q1 + t v;
         w[s_, t_] := P[s] - Q[t];
         First[Minimize[Norm[w[s, t]], {s, t}]]
] 

Here, u and v are unit vectors pointing from p1 to p2 and q1 to q2, respectively. The functions P and Q represent the two lines. The function w is a vector going between P[s] and Q[t]. Minimize the length of that vector.

If you want a method that uses built-in functions, this was suggested to me by a colleague

distLineToLine[line1_, line2_] := Module[{distFunc},
  distFunc = RegionDistance[Line@line1];
  NMinValue[distFunc[{x, y, z}], {x, y, z} ∈ Line@line2]
  ]

But it isn't so fast,

lines = {{{24.134, -147.93, 193.37}, {-9.1854, -85.421, 
     96.284}}, {{-135.82, -104.67, 
     38.896}, {-177.11, -32.589, -74.348}}};  
distLineToLine @@ lines // RepeatedTiming

(* {0.17, 140.358} *)

The best answer here, in my opinion, is the code from the OP, which is about an order of magnitude faster than the accepted answer, although the code doesn't look very pretty. You can speed it up even more by compiling it (code pasted to gist to save space here)

distSegToSegComp = << "https://git.io/vDJMK"

Now compare the timing,

distSegToSeg @@ lines // RepeatedTiming
segsegdis @@ lines // Apply[Subtract] // Norm // RepeatedTiming
distSegToSegComp @@ lines // RepeatedTiming
(* {0.000054, 140.358} *)
(* {0.00012, 140.358} *)
(* {3.7*10^-6, 140.358} *)

So the compiled version wins hands down.