设P1P2线段 P1(P1.x,P1.y,P1.z),P2(P2.x,P2.y,P2.z),在P1P2上有一点为Q，

           Q=P1+t1*(P2-P1);

  设P3P4线段 P3(P3.x,P3.y,P3.z),P2(P4.x,P4.y,P4.z),在P3P4上有一点为Q，

           Q=P3+t2*(P4-P3);

则求AB，CD的距离就变成了QP的长度了，也就是f(t1,t2)的最小值了。

求f(t1,t2)的偏导数，解得极值点，然后检查范围0<t1<1,0<t2<1;

若t1,t2在（0，1）范围内，求f(t1,t2)。

否则求min(P1P3，P1P4，P2P3，P2P4)。

matlab代码如下：

function [ distance] = minDis( P1,P2,P3,P4)
%  计算空间两线段最小距离
%   P1,P2为第一个线段的两个端点
%   P3,P4为第二个线段的两个端点
%   distance为两个线段的最小距离
%计算过程临时变量
e1=P2-P1;
e2=P3-P4;
e3=P1-P3;
%计算二元一次方程组
K=[
    2*e1'*e1 2*e1'*e2
    2*e1'*e2 2*e2'*e2
    ];
B=[-e1'*e3 -e2'*e3]';
T=K\B;  %T=K^(-1)*B
%若两直线的公垂线的垂点在线段上，则直接求公垂线的长度
if (T(1)>0)&(T(2)>0)&(T(1)<1)&(T(2)<1)
    V=T(1)*e1+T(2)*e2+e3;
    distance=sqrt(V'*V);
%若两直线的公垂线的垂点不在线段上，则分别求两线段末端之间的距离，取最小
else
    d1=sqrt((P4-P1)'*(P4-P1));
    d2=sqrt((P3-P2)'*(P3-P2));
    d3=sqrt((P3-P1)'*(P3-P1));
    d4=sqrt((P4-P2)'*(P4-P2));
    distance=min(min(d1,d2),min(d3,d4));
end