2020年华为杯数学建模竞赛F题论文和代码

质心平衡供油策略设计：贪心和搜索的结合

在飞行器飞行过程中，其质心变化对飞行器的控制有着重要的影响，而各个油箱内油

量的分布和供油策略将对飞行器的质心变化产生重要影响。因此，为了最小化飞行过程中

的实际质心和理想质心的最大偏差，本文根据计划耗油速度和理想质心对供油策略和油量

分布进行求解。具体地，本文从四个角度对飞行器质心偏差算法进行说明：对飞行器存在

俯仰角的情况下的质心进行建模求解，对给定初始油量，耗油速度和理想质心对供油策略

进行建模求解，对给定剩余油量对初始和剩余油量分布进行建模求解，对固定理想质心为飞行器坐标系原点，且飞行器存在俯仰角的情况下对供油策略进行建模求解。

针对问题 1，对飞行器存在俯仰角的情况下的质心进行建模求解。本文依据数学几何知识，在每个油箱的中心建立坐标系，求解每个油箱各自的质心，然后进行坐标系转换得到飞行器坐标系下的油箱质心位置，再将油箱视为质点，求解得到飞行器的质心位置。

针对问题 2，对给定初始油量，耗油速度和理想质心对供油策略进行建模求解。本文利用贪心的思路设计逐时间节点的供油策略最优算法，并忽略连续供油的限制得到最优

解，再通过数据分析得到规律，利用聚类的方法修正最优解，迭代得到满足所有约束的最

优可行解。为了得到更少的质心偏差，本文允许飞行器耗油存在一定程度的浪费。在允许浪费20%计划耗油的基础上，本文通过逐时间节点的供油策略最优算法和利用聚类方法最优解修正，得到完整的供油策略，其中飞行过程中飞行器实际质心和理想质心的欧氏距离

最大值为

6.95978 10−5

米，其最大值在第 5365 秒产生，而飞行过程的整体耗油为

6441.746226 千克，相比于计划耗油 6441.524212 千克，仅浪费 0.222014 千克。而在相同条件和算法下，不允许浪费的供油策略得到的欧氏距离最大值为0.1246米，其最大值在第4159 秒产生，相较于允许浪费时的情况，增大了约 1790 倍。因此允许耗油的策略以少量的燃油浪费带来了质心偏移的大幅度降低，更值得被采用。

针对问题 3，不仅需要对给定剩余油量对初始和剩余油量分布进行建模求解，还需要根据耗油速度和理想质心进一步地对供油策略进行建模求解。本文从剩余油量的分布进行

求解，反向推导完整飞行过程的供油策略，并根据最大质心偏移的油量数据修正剩余油量

分布，直到得到最终的油量分布和供油策略。根据

t

=

7200 s

时刻数据，剩余油量的初始解

为(0,54,  137,  111,  120,  428)千克，通过迭代修正为(0.0,  47.18759,  131.5116,  144.8741, 130.9014, 395.5253)千克。通过反向推导得到油箱1-6初始化的油量为(147.50475, 1214.0329, 2019.6,  2254.2,  1560.29485,  1020)千克。此时飞行整体耗油 7298.005094kg，计划耗油6805.174669kg，合计浪费约7.2%。飞行器质心与理想质心距离的最大值为0.119036m，发生在第5372秒。

针对问题 4，对固定理想质心为飞行器坐标系原点，且飞行器存在俯仰角的情况下对供油策略进行建模求解。由于俯仰角的存在，质心计算公式从线性关系的计算中变成了分段式非线性计算，利用问题2的求解思路难以实现。对于不同的俯仰角，对于y轴和z轴的计算均使用平飞状态近似，引入最大的误差为1.8%，在误差允许的情况下能极大的简化计算，提高计算速度。另一方面，将完整的飞行过程按 60 秒为一个时间段划分，强制每一个时间段使用同样的油箱供油，其油箱使用方案有34种。求解采用不同方案60秒供油策略，选择质心偏差和最小的方案作为该 60 秒的最终供油策略，使之在可行解的基础上寻找最优解。最终得到的实际质心与理想质心的欧氏距离最大值为0.038594米，这一最大距离产生在第 3839 秒。其中飞行过程总消耗油量 7565.579365 千克，相比于计划耗油7035.545163千克，浪费7.53%。

2.模型假设

             假设1：可以将油箱简化为长方体且固定在飞行器内部，油箱的长、宽、高的三个方向与飞行器坐标系的x,y,z轴三个方向平行。

 	假设2：以飞行器不载油时的质心	0c	，为原点建立惯性坐标系、飞行器坐标系；
	假设3：第i个油箱的供油速度上限为			U   ( i	U  i	0	)，每个油箱一次供油的持续时间

不少于60秒。

             假设4：主油箱 2、3、4、5 可直接向发动机供油，油箱 1和油箱 6 作为备份油箱分别为油箱2和油箱5供油，不能直接向发动机供油。

             假设5：由于受到飞行器结构的限制，至多2个油箱可同时向发动机供油，至多3个油箱可同时供油。

             假设6：飞行器在执行任务过程中，各油箱联合供油的总量应至少满足发动机的对耗油量的需要，若某时刻供油量大于计划耗油量，多余的燃油可通过其它装置排出飞行器。

             假设7：飞行器的飞行过程中只有俯仰和平飞的情况。飞行器的俯仰将导致各油箱相对地面的姿态发生偏斜，在重力作用下，油箱的燃油分布随之变化。

             假设8：飞行器飞行过程是平稳的，加速度为0，即油箱内燃油液面始终与水平面平行。

             假设9：油箱中不会存在油珠飞溅、油珠挂壁等情况，且油箱内的燃油可以被完全使用。

     假设10：燃油是完全密度均匀的。

            假设 11：假设飞行器飞行过程中，每个油箱为发动机或其他油箱供油是瞬时完成的，且输送频率为每秒1次（或者不输送）。

4. 问题1模型的建立与求解

4.1问题1描述与分析
          问题1给出了某次任务中飞行器的6个油箱的供油速度及飞行器在飞行过程中的俯仰角变化数据，根据飞行器的简化模型和假设条件，建立物理模型对不同情况作分类讨论即可解决问题，如图4-1给出了根据题中所给数据得到的油箱相对位置图。

图4- 1 油箱三维位置示意图

为了分别求解每个油箱在时间 t 的质心位置，对每个油箱分别建立坐标系。以一个油

箱为例，以其在中心位置为原点，建立坐标系

O

−x y z ' ' '

，如图4-2所示。

图4- 2 三维油箱坐标系示意图

8

由于不考虑飞行器的左右偏转，故油箱内的燃油形状始终是柱体，在飞行器俯仰角不

同、油箱内燃油体积不同的情况下，其形态可以分为三棱柱、四棱柱、五棱柱三种情况。

由于燃油是密度均匀的，故在油箱坐标系内，无论燃油处于何种形态， 'y 方向上的坐标始

终为0，可以不用考虑。另一方面，燃油棱柱质心在 'x 方向和 'z 方向上的坐标与棱柱在 ' x Oz '平面上的截面的质心在 'x 方向和 'z 方向上的坐标相同，因此，我们可以将问题进一步简化为求解燃油截面图形的质心。

因此将单一油箱的三维质心问题降维成为xOz平面上求解多边形二维质心问题。根据油箱中所剩余油量不同，共存在三种不同的二维图形，分别为三角形，四边形，五边形。

以飞机向上仰冲时为例，示意图如图4-3所示。

z

x

O

图4- 3 所剩油量为三角形示意图

          当油箱中所剩余油量较少时，会在油箱的xOz平面形成如上图 4-3 所示的三角形，假设油是质量分布均匀且稳定的情况下。此时油箱内油的质心是三角形油面的质心。

z

x

O

图4- 4 所剩油量为四边形示意图

        当油箱中所剩余油量较少时，会在油箱的xOz平面形成如图4- 4所示的四边形，假设油是质量分布均匀且稳定的情况下。此时油箱内油的质心是四边形油面的质心。 同时若油箱内油量继续增加，则会在油箱内形成五边形的截面，其示意图如下图4-  5所示。

此次分别在matlab与Visual Studio 2019通过C++求解整体质心偏移模型。两种方式解得数据相同。解得质心在x轴、y轴、z轴的偏移情况如图4-16所示。

图4- 16 质心各方向偏移示意图

根据上图所示，发现质心在 x 轴向上的偏移量最大，x轴向的偏移对整体的偏移影响最大。y轴方向的质心偏移仅与各个油箱所剩油量相关，在单独计算各个油箱质心时，y轴

方向质心不变，所以y轴质心偏移量也较小。同时z轴向的质心偏移最为稳定，变化幅度

变化量最小。将三个方向的偏移合成一个方向，求解整体的质心偏移量示意图如图4-17所示。

图4- 17 质心整体偏移示意图

如上图所示，质心整体偏移呈现先上升，再下降，然后再上升，最后下降的趋势。基本与 x 轴向的偏移趋势相同。同时质心整体偏移的最大值发生在第 5177 秒，最大偏移量为1.1091。也与x轴的最大值出现时刻大致吻合。同样说明x轴向的偏移对整体的偏移影响最大。

          问题2要求规划满足任务需求的油箱供油策略，目标是使得飞行器每一时刻的质心位

置与理想质心位置

c t 2( )

的欧氏距离的最大值达到最小。题目给定了飞行器计划耗油速度及

飞行器在飞行器坐标系下的理想质心位置数据。与问题1不同的是，问题2限定了飞行器俯仰角始终为 0，简化了质心求解过程，可以认为在油箱坐标系下，油箱的质心位置只会存在z方向上的变化。由于飞行器每个油箱内的初始油量已经给定，我们需要确保油量始终足够飞行器完成本次飞行任务，在此基础上寻找问题的最优解。

问题2的目标在于使整体情况的质心偏移最小，但探索完整飞行过程中的质心偏移最大值最小需要求解大量的参数，例如每个时刻的油箱供油速度和油箱的选择，这使得求解

十分困难。另一方面，在实际飞行器飞行过程中，后续飞行过程的油耗和质心偏移难以估

测，因此通过当前油量和当前的油耗分析进行供油策略的设计是一个比较合理的可行方

案。基于此，我们对模型进行简化，每个时刻仅根据当前时刻剩余油量和需求油量求解供油速度和油箱选择，即用每个时刻的局部最优解代替全局最优解，使得问题得以解决。

5.2 模型建立

          （1）目标函数
          在问题2中，变量是每个油箱每一时刻的供油速度，目标是使得当前时刻的燃油被消耗后，质心位置与当前时刻质心理想位置的欧氏距离最小，依此建模，目标函数如下：

(

x real

min  f

= ||

c t ( )

−

c t ( ) ||

2

( ))

2

+

(

y real

( )

−

y ideal

( ))

2

+

(

z

real

( )

−

z ideal

( ))

t

(5-1)

3 模型实现与求解

5.3.1 算法设计
如果每一时刻都考虑各油箱供油策略要有利于飞行器完成一次飞行任务的全局耗油

需求，情况过于复杂。因此，我们考虑采用贪心算法的思想，逐个时间节点考虑当前情况的最优解。另一方面，约束 5 的引入使得需要对 60 秒的过程整体求解，极大的增加了求解器的压力，为了简化模型方便计算，同样引入贪心的思想，暂时搁置约束5的限制，先求出逐个时间节点的最优解，再对该解进行修正迭代以得到满足约束5的最终结果。

对于逐个时间节点的求解，我们使用了 LINGO 数学软件进行求解。LINGO 是 Linear Interactive and General Optimizer的缩写，即“交互式的线性和通用优化求解器”，由美国LINDO系统公司（Lindo System Inc.）推出的，为非线性规划和线性规划提供了很好的求解方案。问题二的求解模型已有前文给出，其模型为MINLP，即混合整数非线性规划，可以利用 LINGO 软件很好的求解。此外，题目给出的精度十分高，某些初始状态数据达到1e-9的量级，因此对LINGO的非线性求解器参数进行修改，其中Final Nonl Feasibility Tol设置为1e-16，Nonlinear Optimality Tol设置为1e-17，并且开启相关策略加速求解和提高求解精度，其设置可参考下图进行。

          在得到逐个时间节点的最优解方案后，需要对其进行修正以满足约束5的限制。其基于以下两个现象：
          1)  数据出现“集群”现象：如下表所示，第 95-100 秒的供油决策均倾向于 2号油箱供油，可以认为最终方案中该时间段均由2号油箱供油。同时存在少数数据不完全满足该现象，例如表中的 95 时刻的数据，猜想其可能是求解器的精度和结果的不稳定性导致，将其用2号油箱供油的方案相比于原方案带来的误差可以容忍（1e-4级别），并且几乎不会对整个飞行过程的实际质心与理想质心最大偏移产生影响。因此我们可以通过类似聚类

的方法，对逐个时间节点的最优解方案进行修正。

表格5- 1部分时间节点的供油策略

时刻	1号油箱供	2号油箱供	3号油箱供	4号油箱供	5号油箱供	6号油箱供
	油速度	油速度	油速度	油速度	油速度	油速度
95	0	0.009824	0.000232	0	0	0
96	0	0.01017	0	0	0	0
97	0	0.010286	0	0	0	0
98	0	0.010403	0	0	0	0
99	0	0.010521	0	0	0	0
100	0	0.01064	0	0	0	0

2)  修正带来的误差积累：虽然现象 1 表示其单个时间节点上的修正不会带来不可容忍的偏差，但随着时间的积累，其误差存在积累而使得后续数据不可信，另外修正使得油

箱油量数据偏差大，进一步加剧了误差。因此在完成对一定时间段供油策略的修正后，需要对后续数据进行重新求解，保证后续逐时间节点最优解的可靠性。

基于以上两个现象，我们对供油方案进行修正，其流程如下所示。其中初始化供油策略即对所有油箱不设置约束，然后对求解逐个时间节点的最优解方案。再利用聚类的思路，对数据进行修正，此处我们设置了 1000 秒的时间段，即修正 0-1000 秒的供油策略，然后将前1000秒的供油策略固定，重新求解逐时间节点的最优解，再修正下一个1000秒的时间段，直到全部数据修正完成，输出供油策略。

根据上述结果，我们可以得到六个油箱的供油策略和选择方案。其中 1 号油箱在

[387, 864]时间段供油，2号油箱在[95, 2056]时间段供油，3号油箱在[2802, 5251]时间段供油，4号油箱在[1422,  3328]时间段供油，5号油箱在[3959,  5987]时间段供油，6号油箱在[3959, 5365]时间段供油。

将4个主油箱的供油速度相加，得到了4个主油箱的总供油速度（每秒一组），供油曲线如图5-6所示。与如图5-7所示的计划耗油曲线相比，二者基本可以重合。通过计算，4个主油箱的总供油量为6441.746226千克，而发动机实际共耗油6441.524212千克，说明根据 5.3.1 所述的算法得到的供油策略，在保证实际质心与理想质心的偏移距离极小的前提下，导致的燃油浪费极少，浪费总量仅为0.222014千克。

图5-6问题二4个主油箱的总供油速度曲线

图5-7问题2计划耗油速度

图5- 8实际质心与理想质心的欧氏距离

计算每一秒钟飞行器的实际质心与理想质心的欧氏距离，作出曲线如图 5-8 所示。实际质心与理想质心的欧氏距离最大值为6.959780e-05米，这一最大距离产生在第5365秒。可以认定，5.3.1所述的算法求出了一种满足实际质心与理想质心的欧氏距离最大值很小的供油策略，完成了预期目标，且浪费的燃油极少。

为了进一步说明上述结论，我们设置不允许燃油浪费，即在5.2所述模型的约束2中，将设定为 1，得到了一种新的供油策略。使用新的供油策略计算实际质心与理想质心的欧氏距离，得到的曲线如图5-9所示。

图5- 9不允许浪费燃油得到的实际质心与理想质心的欧氏距离

          不允许浪费燃油得到的实际质心与理想质心的欧氏距离最大值为0.1246米，这一最大距离产生在第4159秒。这一最大距离相较于允许浪费时的情况，增大了约1790倍。

综上，使用所建立的模型和求解算法得到了一种供油策略，完成问题2所要求的任务后，4 个主油箱的总供油量为 6441.746226 千克，浪费燃油仅 0.222014 千克，实际质心与理想质心的欧氏距离最大值为6.959780e-05米。

代码

 import os   

4.  from glob import glob   

5.  import subprocess as sp   

6.  import sys   

7.  import openpyxl   

8.  import time   

9.  from tqdm import tqdm    

10. class PowerShell:   

11.     # from scapy   

12.     def __init__(self, coding, ):   

13.         cmd = [self._where('PowerShell.exe'),   

14.                "-NoLogo", "-NonInteractive",  # Do not print headers   

15.                "-Command", "-"]  # Listen commands from stdin   

16.         startupinfo = sp.STARTUPINFO()   

17.         startupinfo.dwFlags |= sp.STARTF_USESHOWWINDOW   

18.         self.popen = sp.Popen(cmd, stdout=sp.PIPE, stdin=sp.PIPE, stderr=sp.STDOUT, sta

rtupinfo=startupinfo)   

19.         self.coding = coding   

20.    

21.     def __enter__(self):   

22.         return self   

23.    

24.     def __exit__(self, a, b, c):   

25.         self.popen.kill()   

26.    

27.     def run(self, cmd, timeout=15):   

28.         b_cmd = cmd.encode(encoding=self.coding)   

29.         try:   

30.             b_outs, errs = self.popen.communicate(b_cmd, timeout=timeout)   

31.         except sp.TimeoutExpired:   

32.             self.popen.kill()   

33.             b_outs, errs = self.popen.communicate()   

34.         outs = b_outs.decode(encoding=self.coding)   

35.         return outs, errs   

36.   

37.     @staticmethod   

38.     def _where(filename, dirs=None, env="PATH"):   

39.         """Find file in current dir, in deep_lookup cache or in system path"""   

40.         if dirs is None:   

41.             dirs = []   

42.         if not isinstance(dirs, list):   

43.             dirs = [dirs]   

44.         if glob(filename):   

45.             return filename   

46.         paths = [os.curdir] + os.environ[env].split(os.path.pathsep) + dirs   

47.         try:   

48.             return next(os.path.normpath(match)   

49.                         for path in paths   

50.                         for match in glob(os.path.join(path, filename))   

51.                         if match)   

52.         except (StopIteration, RuntimeError):   

53.             raise IOError("File not found: %s" % filename)   

54. class Fly:   

55.     def __init__(self):   

56.         self.tanks_location = [[8.91304348,1.20652174,0.61669004],   

57.                                 [6.91304348,-1.39347826,0.21669004],   

58.                                 [-1.68695652,1.20652174,-0.28330996],   

54

59.                                 [3.11304348,0.60652174,-0.18330996],   

60.                                 [-5.28695652,-0.29347826,0.41669004],   

61.                                 [-2.08695652,-1.49347826,0.21669004]]   

62.         self.tanks_size     = [[1.5, 0.9, 0.3],   

63.                                 [2.2, 0.8, 1.1],   

64.                                 [2.4,1.1,0.9],   

65.                                 [1.7,1.3,1.2],   

66.                                 [2.4,1.2,1],   

67.                                 [2.4,1,0.5]]   

68.         self.tanks_rest     = [0.3,1.5,2.1,1.9,2.6,0.8]   

69.         self.speed = [1.1, 1.8, 1.7, 1.5, 1.6, 1.1]   

70.         self.mess = 3000   

71.         self.density = 850   

72.    

73.     def calEcomsume(self, comsume):   

74.         ecomsume = [0] * 6   

75.         for i in range(6):   

76.             ecomsume[i] = comsume[i]   

77.         ecomsume[1] -= comsume[0]   

78.         ecomsume[4] -= comsume[5]   

79.         return ecomsume   

80.    

81.     def UpdateComsume(self,ecomsume,time):   

82.         for i in range(len(ecomsume)):   

83.             self.tanks_rest[i] -= ecomsume[i]*time / self.density   

84.             if self.tanks_rest[i] < 0 :   

85.                 if abs(self.tanks_rest[i]) > 1e-5:   

86.                     print("something wrong")   

87.                 else:   

88.                     self.tanks_rest[i] = 0   

89.    

90.     def calC(self):   

91.         local = [0,0,0]   

92.         for i in range(6):   

93.             for j in range(2):   

94.                 local[j] += self.tanks_rest[i] * self.density * self.tanks_location[i][

j]   

95.         for i in range(6):   

96.             j = 2   

97.             local[j] += self.tanks_rest[i] * self.density * (self.tanks_location[i][j]

- self.tanks_size[i][j]/2 +   

98.                                                      self.tanks_rest[i] / (self.tanks_s

ize[i][0]*self.tanks_size[i][1]*2))   

99.         all_mess = self.mess