第26卷(1996)第1期 、模型假设及说明 1.不碰撞的标准为任意两架飞机的距离大于8公里 2.飞机飞行方向角调整的幅度不应超过30度 3.所有飞机飞行速度均为每小时800公里; 4.进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应在60公里以上; 5最多考虑6架飞机; 6、不必考虑飞机离开此区域后的状况; 7.飞杋进入控制区域后完全服从地面控制台的调度,飞机未接到指令时保持飞行状 态不变。 8.计算机从记录新进入飞机数据到给各飞机发指令间隔为t1,t1小于0.5分。 9.飞机接到指令后可立即转到所需角度,即不考虑转弯半径的影响。 说明 1.假设3假定所有飞机速度均为800km/h,是出于对问题的简化。我们将在模型的 推广中给出飞机速度各不相同时的对策。 2.假设4是必要的,否则可以给出无解的例证 如图所示。对该假设可作如下解释:飞机在区域D外 靠机上雷达自动保持与其他飞机距离大于60km,进 入区域冂后由地面控制台进行统一控制,保证飞机距 离大于8km JJ↓」J 3.假设3中6架飞机的假设是足够多的。以世界 最繁忙的国际航空港之一希恩罗机场邻近区域为例, 假设飞机在区域D作水平飞行,即知该区域内无机场。设在希思罗机场起降的飞机有 半穿过该区域,希思罗机场1992年起降总架次为225万次(文献6),则平均每小时有 15架飞机穿过该区域。而一架飞机穿过该区域最多参160×9≈0.28小时,则任一时 刻该区域上空飞机架数的期望值不超过4.5架。另外,事实上不同飞机的飞行高度是不同 的,这就进一步减少了该区域同一水平面上飞机的数目。以上讨论虽然稍嫌粗略,但是足 以说明6架飞机的假设是合理的。 1.假设8是因为计算机从接到数据到发出命令间存在一个时滞,该时滞固然越小越 好,但受机器限制,一般不能忽略。我们取0.5分为此时滞的一个上限,以使结果具有实际 意义。当κ<1秒时,可以认为实现的是实时控制,时滞可以忽略不计 5.虽然假设2给出的调整范围为30度,但实践证明,10度的调整范围就已足够(从 后面的模型1也可看出,即使两机相向飞行,各自所须的调整也不超过8度),因而在今后 绝大多数讨论和程序的编制中都将搜索区间定为[-10,10度 2 01995-2004 Tsinghua Tongfang Optical Disc Co, Ltd. All rights reserved
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
数学的实践与认识 四、文中用到的符号及说明 时滞 第;架飞机的坐标 an第i架飞机初始方位角t 时间参数 第i架飞机方位角 1),(a;a,) 时刻t两机距离 cosa:-cosa minD(a,a,)ij两机预计最短跗离 sIna, - sina 飞机速度 △△x / 偏差平方和函数 求精次数 逐步求精搜索法中每次求精每层循环次数 gn(x)1,j飞机最短距离构成的不等式约束 h,(x)关于第i架飞机的等式约束 p(X,r)罚函数 权因子 五、模型的建立及求解 模型一(两架飞机的情形.略) (二)模型二 模型二直接将两机模型运用于多机间题,因为它等价于飞机之间两两不相撞。该模型 讨论区域中有6架以下飞机时的情形.利用了模型1判别相撞的函数 crash.同模型 样采取方向角调整幅度平方和最优为调整原则,从而导出目标函数和约束条件如下: f=∑ inD)(a,a,)≥61 (1=1.2.…6,≠j),>0; 或t<0 其中minD)(x,a,)=( +(-△yS+△xC C+s2S,+△y,)2 △xS-△yC,)2 △yS,+△x,C 为使 crash函数只考虑区域中的飞机相撞情况,我们可作如下修改:因为飞机飞过 该区域时间不超过0.28小时(即飞正方形区域对角线时间),我们可认为仅当minD<8 且θ<κ<0.28小时的时候,飞机在区域中相撞;否则不相撞。在实际计算时,我们更把上 限加大到0.3小时·实际上是使不相撞的条件更苛刻了些,相当于对飞机飞离此区域后 的情况也作了部分考虑,提高了全局控制的安全性。 2 01995-2004 Tsinghua Tongfang Optical Disc Co, Ltd. All rights reserved
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
第26卷(1996)第1期 对该模型我们采用直接搜索法讨论了一般情况下对原问题的求解,可在规定时间内 得到一个近似解,如果放宽时限,则可得一个符合精度要求的最优解 直接搜索法原理十分简单:构造多重循环,对所有可能解进行判断,直接得出在一定 精度范围内无可质疑的最优解。但如不使用任何技巧进行直接搜索,必将耗费大量时间。 以本题为例,若在[-10,10]度范围内进行搜索,步长0.01度,共6层循环,需计算6.4 o次,在486DX66上计算一次循环内函数费时2.7×10秒,则此种算法算完显然是不 可能的 为在30秒内算出一个较精确的解,我们采用了逐步求精的方法。即每次用一定的步 长以较少的循环次数进行“粗选”,在“粗选”出的解附近以减小了的步长进行“精选”逐次 推进直到达到指定精度 设每次求精步长减小的倍数是相同的,则每次求精循环次数也相同,设为x,我们考 虑[-10、10度区间,精度要求达到0.01度,设进行了n次求精,则 0.01=2000 总循环次数L=nXx5=(g2000/(x/2))×x 由此式可知x减小则L减小,但若x太小则可能无法收敛到最优解。经验表明x在8 次以上时才能达到较好的搜索效果 n=(log2000/(log(8/2)=5.4≈5 此时共需搜索最多5×9°=265万次,需费时71.55秒 为将计算时间控制在30秒以内,我们又采取了一些优化方法 将底层循环内判别相撞的函数拆细分装在每层循环下,使在高层发现相撞后可提 前结束循环; b.每进入新一层循环把已积累偏差平方和与已得最小偏差平方和比较,若大则结束 该层循环 b.每进入新一层循环把已积累偏差平方和与已得最小偏差平方和比较,若大则结束 该层循环 这些措施大大减少了平均搜索次数,使得在多数情况下计算时间少于30秒,但程序 不能保证在30秒内结東运算,仍存在一些特异情况使计算时间接近最大耗时。我们又试 用偏差绝对值和来代替平方和,发现对最优解影响有限,未能明显缩短计算时间。 至此可知用直接搜索法在现有机器条件下难以满足题设要求,要利用该解法,地面控 制台必须满足以下条件之一 1)拥有速度至少为486DX66三倍以上的电脑; 2)降低精度要求至0.1度(4次求精步长为7,需时至多10.81秒)。 即使模型2不能直接用于飞行管理,它仍有以下作用: 1)可算岀符合精度要求的最优解供检验其他模型。 2)可在相当短的时间內算出具有一定精度的最优解作为其他算法的初值。当精度要 求为1度时,它算出最优解最多只需4.3秒(两次求精步长分别为7,6),大多情况下运算 2 01995-2004 Tsinghua Tongfang Optical Disc Co, Ltd. All rights reserved
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
数学的实践与认识 时间为2~3秒。 模型三 该模型将原问题归结为一个非线性规划冋题,弁用SUMT算法进行了求解 模型2给岀的解法虽然不能满足题设要求却能在较短时间內给出一个较接近最优 解的可行解。由此可行解岀发,用适当的非线性规划算法可较快得岀满足精度要求的最优 以α(i=1,2…6)为变量,在模型2中我们已经将问题归结为非线性规划问题,主要 约束条件为 (1) )≥64 1.2,…6,≠j) 由于minD(a,a)的形式复杂,求导有困难,我们对(1)作一些改变,将目标函数改为 Ca= cosa,- cosa, S 将变量改为C.,S(i=1,2,…6)共12个,增加等式约東 (C,+cos。i0)2+(S (=1,2,…,6) 使问题(1)变为一个12个变量,21个2次约束条件的目标函数求极小值的问题 min/(X)=>((+S%) s.t. g(X miD2(C1,n,C,n,S,,S,)-64≥0(i,j=1,2,…6i<j) 其中X=(C (2) min -(( 二((,a=C,+(a,)·△y +Sm)2+(C h, (X)=(Ctm,+ cosan)+(S, n +sina)2-1 Himmelblau在文献[1]中列举了三类12种算法,我们对这12种算法进行了比较,主 要是考虑其收敛速度的快慢,其次由于没有现成的软件包可供使用,必须自已编制有关程 序,程序准备时间的长短即算法变为程序的复杂程度也必须考虑。综合上述两方面因素 我们选择了序列无约束最小化方法( Sequential Unconstrained Minimization Technique 简称SUMT算法 SUMT法的基本思想是重复地求解一系列无约東问题,它们的解在极限情况下趋于 非线性规划问题的极小点。算法的详细內容参见文献匚1]2]-34].现概述如下 首先构造罚函数P(x,)=f(x)+nS、+nS(x),其中权因子 〃是一单调递减数列。对每个r值求X使P(X,n)取极小值。设精度要求ε,当X X:|<ε时结束运算,。即为符合所需精度要求的最优解。 2 01995-2004 Tsinghua Tongfang Optical Disc Co, Ltd. All rights reserved
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved
第26卷(1996)第1期 ·15 对此具体问题,我们置 P(X,)=f(X)+r∑h(X)+r(-∑切g,(X) r++,ro= mar110-, 1001, r:=10 E=0.5×102 为此问题编制程序SUMT1对题目中的数据进行运算,初始值用模型2以精度为1 度时算出的最优值代入,结果如下: 初始角 偏转角△α 飞机1243.00 0 飞机2236.00 0.00 飞机3220.54 飞机4159.00 0.42 飞机5230.0 0.00 飞机6 偏差平方和6.9 该程序运行时间为7.2秒。 迄今为止,我们尚未讨论时滞对解的影响。在模型3中, Planesumt的运行时间是稳定 的,约为4秒左右。用模型2的 Solve求粗略最优解耗时在2~-4.3秒,随最优解岀现位置 不同而改变,加上裕量后我们可以在10秒以内给岀最优解,故将时滞定为10秒,由假设 7只需将各架飞机的坐标按原速度方向向前移动10秒内的位移,在此基础上求解即可, 这在程序上实现起来相当简便。考虑时滞后得出的结果如下表所示: 初始角a 偏转角△a 飞机1243.00 飞机2236.00 0.00 飞机3220. 飞机4159.00 飞机5230.00 飞机652.00 偏差平方和11.3292 该模型运行结果在精度和耗时两方面都是令人满意的。我们建议地面控制台采用比 486DX66快8倍以上的机器,这样就可以进行实时控制了。对题目所给例子给出的第 个最优解,可以认为即是在实时控制的假设下给出的 六、模型的检验及误差分析 (一)模型检验 模型二用遍历所有可能解的方式得出最优解的,其解(在可能的精度内)的最优性是 2 01995-2004 Tsinghua Tongfang Optical Disc Co, Ltd. All rights reserved
© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved