问题描述:飞机在飞行过程中,能够收到地面上各个监控台发来的关于飞机当前位置的 信息,根据这些信息可以比较精确地确定飞机的位置。如图3所示,VOR 是高频多向导航设备的英文缩写,它能够得到飞机与该设备连线的角度信息;DME 是距离测量装 置的英文缩写,它能够得到飞机与该设备的举例信息。图中飞机接收到来自 3 个 VOR 给出的角度和 1 个 DME 给出的距离(括号内是测量误差限),并已知这 4 种设备的  x, y 坐标(假设飞机和这些设备在同一平面上)。如何根据这些信息精确地确定当前飞机的 位置? 
 

devops workplace导航 vor导航台_最小值

  (1)问题分析 

 记 4 种设备 VOR1、VOR2、VOR3、DME 的坐标为

devops workplace导航 vor导航台_最小值_02

 (以 km 为单位),i=1,2,3,4 ;VOR1、VOR2、VOR3 测量得到的角度为

devops workplace导航 vor导航台_devops workplace导航_03

(从图中可以看出,按照航空飞行管理的惯例,该角度是从北开始,沿顺时针方向的角度,取值在

devops workplace导航 vor导航台_devops workplace导航_04

 之间),角度的误差限为

devops workplace导航 vor导航台_ci_05

;DME 测量得到的距离为

devops workplace导航 vor导航台_ci_06

(单位:km),距离的误差限为 4 σ。设飞机当前位置的坐标为

devops workplace导航 vor导航台_devops workplace导航_07

,则问题就是在表 9 的已知数据下计算 

devops workplace导航 vor导航台_devops workplace导航_07

  。

 

devops workplace导航 vor导航台_devops workplace导航_09

(2)模型 1 及求解 

 图中角度 

devops workplace导航 vor导航台_devops workplace导航_03

 是点 

devops workplace导航 vor导航台_最小值_02

  和点

devops workplace导航 vor导航台_devops workplace导航_07

 的连线与 y 轴的夹角(以 y 轴正向为基准,顺时针方向夹角为正,而不考虑逆时针方向的夹角),于是角度 

devops workplace导航 vor导航台_devops workplace导航_03

 的正切 

devops workplace导航 vor导航台_最小二乘_14

    ( 1 )对 DME 测量得到的距离,显然有 

devops workplace导航 vor导航台_devops workplace导航_15

     ( 2 )

直接利用上面得到的 4 个等式确定飞机的坐标 y x, ,这是一个求解超定(非线性) 方程组的问题,在最小二乘准则下使计算值与测量值的误差平方和最小(越接近 0 越 好),则需要求解 

devops workplace导航 vor导航台_最小二乘_16

  ( 3 )

式(3)是一个非线性(无约束)最小二乘拟合问题。很容易写出其 LINGO 程序 如下:

MODEL: 
TITLE 飞机定位模型1; 
SETS: 
VOR/1..3/:x0,y0,cita,sigma; 
ENDSETS 
DATA: 
x0, y0, cita, sigma = 
746  1393 161.2   0.8 
629  375 45.1     0.6 
1571 259 309.0    1.3; 
x4 y4 d4 sigma4=155,987,864.3,2.0; 
ENDDATA 
calc: 
@for(VOR:cita=cita*3.14159/180;sigma=sigma*3.14159/180); 
endcalc 
min=@sum(VOR:@sqr((x-x0)/(y-y0)-@tan(cita)))+@sqr(d4-@sqrt(@sqr(x -x4)+@sqr(y-y4))); 
END

上述程序必须使用全局求解器进行求解,否则求得的是一个局部最优解。用 “LINGO|OPTIONS”菜单命令启动“Global Solver”选项卡上的“Use Global Solver”选项,然后求解,可以得到全局最优解 x=1019.306 ,y= 987.2909  ,对应的目标函
数值为 0.4729562,这里的解受π 的取值影响很大。 

(3)模型 2 及求解

注意到这个问题中角度和距离的单位是不一致的(角度为弧度,距离为公里),因 此将这 4 个误差平方和同等对待(相加)不是很合适。并且,4 种设备测量的精度(误差限)不同,而上面的方法根本没有考虑测量误差问题。如何利用测量设备的精度信息? 这就需要看对例中给出的设备精度如何理解。 一种可能的理解是:设备的测量误差是均匀分布的。以 VOR1 为例,目前测得的角度为

devops workplace导航 vor导航台_ci_17

 ,测量精度为

devops workplace导航 vor导航台_最小二乘_18

,所以实际的角度应该位于区间 

devops workplace导航 vor导航台_ci_19

  内。对其它设备也可以类似理解。由于

devops workplace导航 vor导航台_devops workplace导航_20

 很少,即测量精度很高,所以在相应区间内正切函数 tan 的单调性成立。于是可以得到一组不等式:

devops workplace导航 vor导航台_最小值_21

  也就是说,飞机坐标应该位于上述不等式组成的区域内。 由于这里假设设备的测量误差是均匀分布的,所以飞机坐标在这个区域内的每个 点上的可能性应该也是一样的,我们最好应该给出这个区域的 x和 y 坐标的最大值和最小值。于是我们可以分别以 min x ,  max x ,  min y,  max y为目标,以上面的区域限制条件为约束,求出x 和 y 坐标的最大值和最小值。 
 以 min x 为例,相应的 LINGO 程序为:  

MODEL: 
TITLE 飞机定位模型2; 
SETS: VOR/1..3/:x0,y0,cita,sigma; 
ENDSETS 
INIT: 
x=1000; y=900; 
ENDINIT 
DATA: 
x0, y0, cita, sigma = 
746  1393 161.2   0.8 
629  375 45.1     0.6 
1571 259 309.0    1.3; 
x4 y4 d4 sigma4=155,987,864.3,2.0; 
ENDDATA 
calc: 
@for(VOR:cita=cita*3.14159/180;sigma=sigma*3.14159/180); 
endcalc 
min=x; 
@for(VOR:(x-x0)/(y-y0)>@tan(cita-sigma)); 
@for(VOR:(x-x0)/(y-y0)<@tan(cita+sigma)); 
d4-sigma4 <((x-x4)^2+(y-y4)^2)^.5 ; 
d4+sigma4 >((x-x4)^2+(y-y4)^2)^.5 ; 
END

注意:用 LINGO9 求解非线性问题,必须对决策变量进行初始化,否则 LINGO 可 能找不到可行解。决策变量的初值也有范围限制,取的不合适也可能找不到可行解。 求得的 x的最小值为 974. 8433。类似地(只需要换目标函数就可以了),可得 到 x的最大值为 982.2005, y 的最小值为 717.1614, y 的最大值为 733.1582。 因此,最后得到的解是一个比较大的矩形区域,大致为

devops workplace导航 vor导航台_最小值_22


  (4)模型 3 及求解

模型 2 得到的只是一个很大的矩形区域,仍不能令人满意。实际上,模型 2 假设 设备的测量误差是均匀分布的,这是很不合理的。一般来说,在多次测量中,应该假设设备的测量误差是正态分布的,而且均值为 0。本例中给出的精度

devops workplace导航 vor导航台_devops workplace导航_20

 可以认为是测量

误差的标准差。 在这种理解下,用各自的误差限  

devops workplace导航 vor导航台_devops workplace导航_20

 对测量误差进行无量纲化(也可以看成是一种加权法)处理是合理的,即求解如下的无约束优化问题更合理。 

devops workplace导航 vor导航台_devops workplace导航_25

 由于目标函数是平方和的形式,因此这是一个非线性最小二乘拟合问题。相应的 LINGO 程序为: 

MODEL: 
TITLE 飞机定位模型3; 
SETS: 
VOR/1..3/:x0,y0,cita,sigma,alpha; 
ENDSETS 
INIT: 
x=1000; y=900; 
ENDINIT 
DATA: 
x0, y0, cita, sigma = 
746  1393 161.2   0.8 
629  375 45.1     0.6 
1571 259 309.0    1.3; 
x4 y4 d4 sigma4=155,987,864.3,2.0; 
ENDDATA 
calc: 
@for(VOR:cita=cita*3.14159/180;sigma=sigma*3.14159/180); 
endcalc 
min=@sum(VOR:((alpha-cita)/sigma)^2)+((d4-((x-x4)^2+(y-y4)^2)^.5 )/ sigma4 )^2; 
@for(VOR: @tan(alpha)=(x-x0)/(y-y0) ); 
END

 启动 LINGO 的全局最优求解程序求解,得到全局最优解 x=978.3071,y= 723.9841,对应的目标函数的值为 0.668035。 这里得到的误差比模型 1 的大,这是因为模型 1 中使用的是绝对误差,而这里使用的是相对于精度

devops workplace导航 vor导航台_devops workplace导航_20

 的误差。对角度而言,分母 

devops workplace导航 vor导航台_devops workplace导航_20

 很少,所以相对误差比绝对误差大,这是可以理解的。