查看: 1676|回复: 1

[问题求助] 一组空间点，求一条直线，使这些点到这个直线距离都相等，

字体大小: 正常放大

1 主题	1 听众	3 积分

升级 60%

该用户从未签到

电梯直达

1^#

发表于 2015-1-27 16:44 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

各位大侠，小弟初接触模拟退火算法，现在做一个程序，问题的描述是：
6个空间点，p1-p6.已知各点坐标x1，y1，z1。求一条直线，使该直线到上述点距离都相等。
空间直线可以用参数方程表达，x=at+px，y=bt+py，z=ct+pz（其中a，b，c是直线的方向向量。px，py，pz是直线上一点。t是直线参数。）
假设参数方程已知，我可以通过点到空间直线距离公式计算各点到直线距离。将计算得到的直线距离取平均值，然后用各距离减平均值，以该差值作为收敛目标（也就是各点到直线距离的差距越小越好，越接近目标）
我的问题是我的理论可以实现对直线方程的求解吗？以下是我的matlab程序，收敛不了，请大家帮小弟挑错
% 初始化过程计算初始状态
pp=[361.954086685118,222.958293409968,57.3284760947955,-282.138702545910,-357.764981326235,-321.219260549933;
12.7017055583522,283.682166162709,355.435259006668,221.691115425287,0.0772113586361366,-157.707415860041;
-471.015149591140,-236.049854262889,-150.230911572783,-306.421592445076,-479.622244334522,-514.184148298648;
1,1,1,1,1,1;];
x1=pp(1,

; y1=pp(2,

;z1=pp(3,

;
  a(1)=0;  b(1)=0;  c(1)=1;
  ja(1)=pi/2; jb(1)=pi/2; jc(1)=0;%根据a、b、c的值求得的对应余弦角
  px(1)=0;  py(1)=0;  pz(1)=-860;
for i=1:1:6%计算初始距离
   t(i)=(a(1)*(x1(i)-px(1))+b(1)*(y1(i)-py(1))+c(1)*(z1(i)-pz(1)))/(a(1)^2+b(1)^2+c(1)^2);%求垂足参数
   x(i)=a(1)*t(i)+px(1);%求垂足位置
   y(i)=b(1)*t(i)+py(1);
   z(i)=c(1)*t(i)+pz(1);
   rcy(i)=((x(i)-x1(i))^2+(y(i)-y1(i))^2+(z(i)-z1(i))^2)^0.5;%分别计算点到直线距离
end
  for i=1:1:6
   drcy(i)=rcy(i)-mean(rcy)%计算各距离偏差值
  end
   e(1)=max(abs(drcy));%计算最大误差，作为收敛目标，当e尽量小的时候收敛结束
   qar=0.5*3.14/180;  qbr=0.5*3.14/180;  qcr=0.5*3.14/180;%a，b，c对应余弦角取值范围
   qxr=5;  qyr=5;  qzr=5;%px，py，pz取值范围
   rotet=0.95%设置降温速率
   %计算初始温度
      for j=1:1:120%先进行120次试算，
         kja(j)=pi/2-qar+2*qar*rand;
         ka(j)=cos(kja(j));
         kjb(j)=pi/2-qbr+2*qbr*rand;
         kb(j)=cos(kjb(j));
         kc(j)=(1-ka(j)^2-kb(j)^2)^0.5
            kpx(j)=0-qxr+2*qxr*rand;
            kpy(j)=0-qyr+2*qyr*rand;
            kpz(j)=-860-qzr+2*qzr*rand;
            for i=1:1:6%计算初始半径
               t(i)=(ka(j)*(x1(i)-kpx(j))+kb(j)*(y1(i)-kpy(j))+kc(j)*(z1(i)-kpz(j)))/(ka(j)^2+kb(j)^2+kc(j)^2);
               x(i)=ka(j)*t(i)+kpx(j);
               y(i)=kb(j)*t(i)+kpy(j);
               z(i)=kc(j)*t(i)+kpz(j);
               rcy(i)=((x(i)-x1(i))^2+(y(i)-y1(i))^2+(z(i)-z1(i))^2)^0.5;
               drcy(i)=rcy(i)-360;
            end
         ee(j)=max(abs(drcy));%计算最大误差
      end
      maxee=max(ee); meanee=mean(ee); pp0=0.95;
   ctemperature(1)=(meanee-maxee)/log(pp0);%计算初始温度

%%%%%%%%%%%%%%%%%%%%%%开始进入循环
   j=2;%设置计数器，避免进入死循环
%    while j<10 & e(j-1)>5%迭代10次为实验用，实际上需要迭代不止10次，e（j）一定要小于5这是角度改变迭代后的目标
while j<10 & e(j-1)>0.5%怕计算不出答案，将程序循环过程设为10，收敛目标设的稍微大一点。
%       while j<10 & e(j-1)>3
      mja(j)=ja(j-1)-qar+2*qar*rand;ma(j)=cos(mja(j));
      mjb(j)=jb(j-1)-qbr+2*qbr*rand;mb(j)=cos(mjb(j));
      mc(j)=(1-ma(j)^2-mb(j)^2)^0.5;
      mpx(j)=px(j-1)-qxr+2*qxr*rand;
      mpy(j)=py(j-1)-qyr+2*qyr*rand;
      mpz(j)=pz(j-1)-qzr+2*qzr*rand;
      %赋值进入计算
         for i=1:1:6%计算初始半径
               t(i)=(ma(j)*(x1(i)-px(j))+mb(j)*(y1(i)-py(j))+mc(j)*(z1(i)-pz(j)))/(ma(j)^2+mb(j)^2+mc(j)^2);
               x(i)=ma(j)*t(i)+px(j);
               y(i)=mb(j)*t(i)+py(j);
               z(i)=mc(j)*t(i)+pz(j);
               rcy(i)=((x(i)-x1(i))^2+(y(i)-y1(i))^2+(z(i)-z1(i))^2)^0.5;
               drcy(i)=rcy(i)-mean(rcy)
         end
         e(j)=max(abs(drcy));%计算最大误差
      deltae=e(j)-e(j-1);%计算两次状态，最大误差值
      if deltae<=0 %如果小于0了，那就是状态收敛，所得数据可以接受
            ja(j)=mja(j);jb(j)=mjb(j);%用来下一步迭代，赋初值
            px(j)=mpx(j);py(j)=mpy(j);pz(j)=mpz(j);
            qar=(1-(j-1)*0.02)*qar;qbr=(1-(j-1)*0.02)*qbr;%收敛数据选取范围缩小
            qxr=(1-(j-1)*0.02)*qxr;qyr=(1-(j-1)*0.02)*qyr;qzr=(1-(j-1)*0.02)*qzr;
            ctemperature(j)=rotet*ctemperature(j-1);
            j=j+1;
      else%如果不小于0，未收敛，按照概率形式检查选取范围
               kexp(j)=exp(-deltae/ctemperature(j-1))
            if exp(-deltae/ctemperature(j-1))>rand
               ja(j)=mja(j);jb(j)=mjb(j);%用来下一步迭代，赋初值
               px(j)=mpx(j);py(j)=mpy(j);pz(j)=mpz(j);
               qar=(1-(j-1)*0.02)*qar;qbr=(1-(j-1)*0.02)*qbr;
               qxr=(1-(j-1)*0.02)*qxr;qyr=(1-(j-1)*0.02)*qyr;qzr=(1-(j-1)*0.02)*qzr;
               ctemperature(j)=rotet*ctemperature(j-1);
               j=j+1;
            else

            end
      end

end