随着航空业进入工业4.0时代,提高航班运行效率与安全性成为其第一要务。据国际航空运输协会(International Air Transport Association,IATA)统计,截至2022年年底,全球拥有超过100万个航空货物集装器流通,这些集装器支撑起了高时效性货物、高价值货物、精密货物运输的供应链体系,是经济全球化和快速发展的重要引擎。然而,集装器(unit load device,ULD)装载不规范引发的安全事故频发。IATA数据表明,多数飞机地面破损事故(包括飞机结构破损、倾覆、起飞擦尾等)与ULD装载过程相关,欧洲航空安全局(European Union Aviation Safety Agency,EASA)报告指出,不合理的ULD装载严重影响飞机安全性、机动性和经济性。因此,深入研究航空集装器配载与平衡问题(weight and balance problem,WBP)具有重要的现实意义。
针对航空集装器配载与平衡问题的优化研究,主要围绕建模方法、目标函数、约束条件、算法四个方面。
建模方法方面,根据现有文献,混合整数线性规划模型(mixed integer linear programming,MILP)被多数国内外研究者采用。Sun等
[1 ] 针对内陆散货堆场的存储空间分配问题,提出了MILP模型,并基于逻辑型Benders分解的高效求解算法,将问题分解为处理作业调度的主问题和处理空间分配的子问题。Desai等
[2 ] 针对客改货飞机的航空货物装载问题,建立了0-1混合整数非线性规划模型,采用CPLEX求解器求解线性化的0-1混合整数规划模型。史永胜等
[3 ] 针对货机主货舱集装器装载优化问题,建立混合整数规划模型,采用遗传算法求解,通过实数编码、轮盘赌选择、单点交叉和逆变异等操作寻找最优装载方案。Do Nascimento等
[4 ] 采用整数规划与约束规划相结合的混合建模方法,通过三个模型迭代求解单容器装载问题,采用CPLEX求解器进行迭代优化,当发现不可行解时通过添加割平面约束排除该解集,直至找到最优可行解。Limbourg等
[5 ] 开发了一种MILP模型,以减少集装器在飞机目标重心处的堆积;并以最小化转动惯量为目标,设置零燃油重心偏差为0.01作为约束条件进行限制,减少了集装器在重心方向上的堆积,提高了飞机飞行的机动性,减少了燃油消耗,降低了成本。由此可见,研究者大都使用MILP模型,目前大多数模型并未考虑ULD与舱位之间的物理适配约束,或采用了复杂的不等式约束表达,增加了模型求解难度。
目标函数方面,随着航空货运装载的需求变大,装载要求更加严格,通常会要求同时优化多个目标。故国内外研究者近年来的研究多围绕多目标优化。赵向领等
[6 ] 为优化研究航空公司所遇到的飞机载重平衡问题,建立了以业载量最大和重心偏移量最小为目标的多目标载重平衡模型,在普通遗传算法的基础上添加了进化逆转操作算子,最后运用波音757-200F作仿真优化实例。陈丙成等
[7 ] 构建了一个多航空集装器装箱优化模型,按照遗传算法确定的航空集装器使用顺序逐个装满集装器,将多集装器装箱问题转化为单集装器装箱的序列决策问题;目标函数主要包括使用最少数量的航空集装器完成货物装载、优先选择运输成本较低的集装器类型,以及最大化舱位利用率。Zhao等
[8 ] 针对飞机质量平衡问题提出了MILP模型,该模型旨在最大化总载荷质量并最小化重心偏差,通过加权方式合并为单一目标函数。Sahun等
[9 ] 提出了一个基于启发式规则的货机装载规划模型,采用可视化计算机建模方法,将货舱划分为15个区段并按照多段航线的卸货优先级进行货物分组。目标是减少装载操作时间和优化货物分布,通过最小化装载时间来提高地面服务效率。采用基于规则的启发式决策算法,结合专家经验制定装载规则。Wong等
[10 ] 采用MILP建模方法,构建了集成物理操作系统、优化仿真模块、虚拟现实系统和物联网传感器的闭环数字孪生系统。目标函数为最大化三个加权性能指标的综合得分,分别为货运收入、重心位置、卸货优先级。可以看出,虽然目前有各种航空公司根据自身需求制定的不同目标函数,但是单一的目标函数存在局限性,多目标优化缺少对目标函数的优先级和递进关系的研究。
约束条件方面,航空集装器装载时需要满足各种约束条件,以确保飞机飞行的安全性。Tseremoglou等
[11 ] 提出了一个航空货运动态装载决策支持框架,建模方法上采用预测-优化两阶段方法,涵盖了包括概率容量约束、ULD几何边界约束、质量容量约束、货物稳定性支撑约束、旋转方向约束以及订单完整性约束等约束条件,后通过改进的极值点启发式算法求解三维装箱问题。Sahoo等
[12 ] 构建了一个航空货运整合与合并优化的混合整数线性规划模型,建模方法上采用多目标优化框架,同时考虑货物整合、合并、易腐性处理和碳排放约束;目标函数为最小化总运输成本,约束条件涵盖二元分配变量约束、整合与合并逻辑约束、目标交付时间约束、资源容量约束(体积和质量)、碳排放减少约束以及集装箱温度平衡约束。Brandt
[13 ] 针对航空货物装载规划问题,采用分解建模方法将复杂问题拆分为飞机配置、组装调度、货物装箱和重心平衡四个子问题,分别建立了相应的数学优化模型。目标函数以最小化总成本为核心,同时考虑最大化载货量和收入。约束条件涵盖了容量限制(质量/体积)、时间窗口、三维几何约束(物品不重叠、位置合规)、飞机重心平衡、物品兼容性、装载稳定性等多维度限制。赵向领等
[14 ] 将民航货机装载问题建模为二次指派问题,建立了双目标混合整数规划模型,并采用约束法将双目标转化为单目标优化问题;以基本约束(集装器与舱位的一一对应关系)、质量限制(最大无油质量、起飞质量和业载限制)、安全性约束(无油重心、起飞重心和着陆重心包线范围)以及结构性约束(各舱位和分区的载荷限制),提出了改进的遗传算法,结合最优保存策略确保算法的全局收敛性。Mazur等
[15 ] 针对强异构货物在货盘上的装载,开发一种物理包装顺序,来保证货盘的平稳性。张长勇等
[16 ] 提出了一种以货柜空间利用率最大为目标,考虑货物装载顺序、体积、质量、重心、不重叠等多种实际约束。王素欣等
[17 ] 在针对货物三维装箱问题建立三维装箱模型,加入了考虑转弯时重心约束。可以看出,现有研究并未同时考虑多重约束,且约束需要大量不等式,造成了模型的复杂性。以及多数研究忽略横向平衡约束对飞机稳定性的影响,对于飞行过程中的重心变化考虑不足。
算法方面,航空集装器装载问题属于NP-Hard(non-deterministic polynomial hard)问题,求解NP-Hard问题多使用启发式算法。Nayak等
[18 ] 通过在飞机上安排不同类型、尺寸和质量的货物来充分利用其有效载荷能力,减少燃油消耗,以实现收益最大化;同时提出了量子方法来解决最大化飞机的有效载荷能力和最小化跨多个中途站的装卸操作,主要为基于量子退火和基于门的算法。Emde等
[19 ] 针对集装器不同集合设计不同的启发式算法以减少装载集装器使用的成本。Ma等
[20 ] 设计一种混合遗传算法优化客机腹舱装载空间。Bombelli等
[21 ] 开发了一种自适应邻域搜索算法实现集装器的装卸。陈元文
[22 ] 针对以空间利用率最大为目标的装载问题,根据问题特点,设计改进遗传算法。李国举等
[23 ] 针对航空货运“干-支”联动配载优化问题,建立了以“干-支”最大装载率为优化目标,采用遗传算法开展了 “干-支”联动货物配载的优化研究。可以看出,航空货运集装器的装载问题是一个复杂的多目标组合优化问题,涉及大量的约束条件。然而,由于装载规模庞大、考虑因素复杂且计算量大等问题,仍然缺乏有效的优化方法来解决实际运输中货机装载的各种问题。
综上,现有关于标准航空集装器的装载优化研究仍存在以下问题:优化目标较为单一,与航空公司的实际运控需求存在脱节,以及集装器与舱位之间的适配约束较为宽松。针对上述不足,设计三个递进式的优化目标——业载最大化、俯仰转动惯量最大化,以及综合权衡归一化俯仰转动惯量与重心偏差,并构建相应的优化模型与策略。在综合考虑飞机承载能力、货舱内部空间等物理与几何约束的基础上,通过引入舱位-集装器适配矩阵,设计初始种群生成机制,以产生高质量初始解,并采用改进的遗传算法进行装载优化研究。最后,基于B747-400F飞机不同规模的集装器装载实例,对优化结果进行对比分析,验证本文所提出模型与算法的可行性和高效性,从而实现飞行安全性、收益性、机动性与经济性的全局平衡。
1 标准集装器数学模型建立
1.1 问题描述和假设
标准航空集装器装载问题为:在一定约束条件限制下,将一批待装载的标准集装器装入飞机货舱的舱位上,以获得最大的业载利用率、最小的绕俯仰轴的转动惯量和最小的重心偏差。
基于航空装载问题的复杂性,先做如下假设。
(1)集装器到达同一站点。
(2)集装器在舱位上沿着机身纵向放置。
(3)将集装器看作规则的长方体,将形心位置作为重心位置。
1.2 装载问题求解策略
1.2.1 航司运控目标的转化
为使航空集装器的配载优化模型和优化客观、科学地反映航空公司航司运控目标优先级,通过优化模型的目标函数和约束条件对航司运控目标进行定义与量化,具体如下。
(1)飞行安全性。飞机安全性是航司运行控制的首要目标,通过飞机的业载限制、舱位与集装器的适配约束、重心包络约束、结构约束和横向不平衡等强约束。将飞机始终限定在正常性能范围内,确保从推出、滑行到着陆全过程的安全性。这些强约束贯穿优化策略的三个阶段,保证各阶段解及全局最优解的可行性。
(2)飞机机动性。飞机机动性指飞机快速改变运动状态的能力,包括速度、高度和方向机动性。本文中重点优化货机的俯仰方向机动性。在业载和发动机性能一定时,通过最小化俯仰转动惯量,使货机能快速调整俯仰角以获得更大升力,这对实现越障以及飞机复飞至关重要。同时,通过横向不平衡约束限制左右业载差,保证横向机动性。
(3)运输经济性与收益。通过定义飞机业载最大化的目标函数来实现运输收益的最优。设定符合航司经济利益的目标起飞重心,最小化实际重心与目标重心偏差,优化运输经济性。同时,基于能量转化原理,优化最小化俯仰转动惯量不仅能提升机动性,还能减少飞机姿态改变所需能量,降低燃油消耗,进一步提升航司收益。
1.2.2 航司实际运控与运营目标的优化策略
针对现有方法将多个目标函数加权转化为单一目标,而因各目标之间量纲不一致、权重难确定导致优化过程对目标函数不敏感,优化结果不能客观反映真正优化意图的问题。
如
图1 所示为三阶段递进优化方法流程图。按飞行安全性、收益、机动性和经济性在航司实际运行控制和运营中的优先级,以及与飞机配载属性参数之间的关系,设计了一种三阶段递进优化策略。
该策略首先以飞机业载最大为目标,以飞机承载能力、重心偏移等安全相关限制为强约束,通过优化得到满足所有强约束且业载利用率最大的ULD装载方案。
根据上一阶段的最优解包含的ULD集合更新舱位-集装器适配矩阵,以上一步的最优装载方案中的ULD集合为装载对象,以俯仰转动惯量最大为目标,在满足所有安全性约束的条件下优化得到业载最大情况下的最大俯仰转动惯量,并基于此最大值对俯仰转动惯量目标函数进行归一化处理。该处理可确保归一化后的俯仰转动惯量目标函数的取值在[0,1]区间内,且可以和其他目标函数统一量纲。
最后,以归一化的俯仰转动惯量和重心偏离目标位置的平均气动弦长百分比最小为目标,根据上一阶段的最优解包含的ULD集合更新舱位-集装器适配矩阵,以上一步的最优装载方案中的ULD集合为装载对象,在满足所有安全性约束的条件下,优化得到载重利用率最大情况下的俯仰转动惯量和重心偏差加权和最小的ULD装载方案,实现航司以飞行安全性为基石,机动性、经济性和收益最大化的实际运控目标。
1.3 优化模型构建
1.3.1 决策变量
使用0-1矩阵来定义集装器集合与舱位集合的装载关系。0-1矩阵的元素取值为
(1) xij =$\left\{\begin{array}{ll}0,& \forall \mathit{i}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}\}\\ 1,& \forall \mathit{j}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}\}\end{array}\right.$
式(1)中:N ULD 为待装载的ULD数量; N POS 为需要装载的位置数量;xij 为第i 个集装器与第j 个舱位之间的装载关系。如果第i 个集装器装载到第j 个舱位上,则xij =1,否则xij =0。由此,可定义决策变量矩阵X 为
(2) X =$\left[\begin{array}{llll}{\mathit{x}}_{11}& {\mathit{x}}_{12}& \dots & {\mathit{x}}_{1{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\\ {\mathit{x}}_{21}& {\mathit{x}}_{22}& \dots & {\mathit{x}}_{2{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\\ \mathit{︙}& \mathit{︙}& & \mathit{︙}\\ {\mathit{x}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}1}& {\mathit{x}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}2}& \dots & {\mathit{x}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\end{array}\right]$
1.3.2 目标函数
以下分优化策略的三个阶段来分别定义相应的目标函数。
1)第一阶段的目标函数
业载利用率最大化的目标函数为
(3) $\left\{\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}{\mathrm{f}}_{1}\left({\mathrm{x}}_{\mathrm{i}\mathrm{j}}\right)=\stackrel{{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}{\sum _{\mathrm{i}=1}}\stackrel{{\mathrm{N}}_{\mathit{P}\mathit{O}\mathit{S}}}{\sum _{\mathrm{j}=1}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}{\mathrm{w}}_{\mathrm{i}}=\mathrm{S}\mathrm{U}\mathrm{M}\left({\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}\right)\\ {\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}=[{\mathrm{w}}_{1}{\mathrm{w}}_{2}\dots {\mathrm{w}}_{{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}]\end{array}\right.$
式(3)中:f 1 (xij )为第一阶段的目标函数,即业载的目标函数;wi 为第i 个集装器的质量;W ULD 为所有集装器的质量矩阵。
为了便于后续适应度函数的定义,同时结合改进的遗传算法的应用需求,将问题的目标函数转化为目标函数最小化来处理,上述目标函数转化为
(4) minf' 1 (xij )=1-$\frac{{\mathit{f}}_{1}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}{{\mathit{W}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}$
式(4)中:W max 为飞机在飞行计划确定情况下的最大业载量; f' 1 (xij )为经过归一化处理后的第一阶段目标函数。
2)第二阶段的目标函数
俯仰转动惯量最大化的目标函数为
(5) $\left\{\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}{\mathrm{f}}_{2}\left({\mathrm{x}}_{\mathrm{i}\mathrm{j}}\right)=\stackrel{{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}{\sum _{\mathrm{i}=1}}\stackrel{{\mathrm{N}}_{\mathit{P}\mathit{O}\mathit{S}}}{\sum _{\mathrm{j}=1}}{\mathrm{w}}_{\mathrm{i}}({\mathrm{a}}_{\mathrm{j}}-{\mathrm{I}}_{\mathit{D}}{)}^{2}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}={\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}{\mathit{A}}_{\mathit{p}}\\ {\mathit{A}}_{\mathit{p}}=[{\mathrm{a}}_{1}{\mathrm{a}}_{2}\dots {\mathrm{a}}_{{\mathrm{N}}_{\mathit{P}\mathit{O}\mathit{S}}}{]}^{\mathit{T}}\end{array}\right.$
式(5)中:f 2 (xij )为第二阶段的目标函数,即表示俯仰转动惯量的目标函数;aj 为第j 个舱位的平衡力臂; I D 为目标重心的力臂;Ap 为第p 个舱位的平衡力臂列向量。
出于同样的目的,参照第一阶段目标的转化,将第二阶段的目标函数进行转化,即
(6) minf' 2 (xij )=1-$\frac{{\mathit{f}}_{2}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}{{\mathit{f}}_{2}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}$
式(6)中:${\mathit{f}}_{2}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ xij )为各舱位装载满足适配约束且最大质量的集装器时的全机俯仰转动惯量;f' 2 (xij )为经过归一化处理后的第二阶段目标函数。
3)第三阶段的目标函数
起飞重心偏差最小化的目标函数表达式为
(7) min${\mathit{f}}_{3}^{\mathrm{c}\mathrm{g}}$ (xij )=|G TOW -G TAR |
式(7)中:${\mathit{f}}_{3}^{\mathrm{c}\mathrm{g}}$ xij )为第三阶段的目标函数,即起飞重心偏差目标函数;G TOW 和G TAR 分别为实际起飞重心和目标重心,用平均气动弦长百分比表示。
对俯仰转动惯量和重心偏差这两个目标函数分别进行归一化,即
(8) $\left\{\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{\mathit{f}}_{2}^{\mathrm{i}\mathrm{m}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)=\frac{{\mathit{f}}_{2}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}{{\mathit{f}}_{2}^{\mathrm{o}\mathrm{p}\mathrm{t}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}\\ \mathrm{m}\mathrm{i}\mathrm{n}\mathit{f}{\mathit{\text{'}}}_{3}^{\mathrm{c}\mathrm{g}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)=\frac{{\mathit{f}}_{3}^{\mathrm{c}\mathrm{g}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}{{\mathit{f}}_{{3}_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{c}\mathrm{g}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)}\end{array}\right.$
式(8)中:${\mathit{f}}_{2}^{\mathrm{i}\mathrm{m}}$ ${\mathit{f}}_{2}^{\mathrm{o}\mathrm{p}\mathrm{t}}$ f ${\mathit{\text{'}}}_{3}^{\mathrm{c}\mathrm{g}}$ xij )为起飞重心偏差最小化归一化处理目标函数;${\mathit{f}}_{{3}_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{c}\mathrm{g}}$ xij )为起飞重心包络线的各条边界与目标重心之间的最大差值。
综上,第三阶段的目标函数定义为
(9) minf 3 (xij )=ρ ${\mathit{f}}_{2}^{\mathrm{i}\mathrm{m}}$ (xij )+(1-ρ )f ${\mathit{\text{'}}}_{3}^{\mathrm{c}\mathrm{g}}$ (xij )
式(9)中:ρ 为权重值,取ρ =0.5。
1.3.3 约束条件
1)位置约束
每个舱位最多只能装载一个集装器,约束条件定义为
(10) $\left\{\begin{array}{l}{\mathrm{C}}_{1}^{}:\mathit{S}\mathit{U}\mathit{M}\left({\mathit{X}}_{\mathit{i}}\right)\le 1\\ \mathrm{j}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}\}\end{array}\right.$
式(10)中:${\mathrm{C}}_{1}^{}$ Xi 为X 的第i 列向量。
每个集装器最多只能放置一个舱位上,约束条件定义为
(11) $\left\{\begin{array}{l}{\mathrm{C}}_{2}^{}:\mathit{S}\mathit{U}\mathit{M}\left({\mathit{X}}_{\mathit{j}}\right)\le 1\\ \mathrm{i}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}\}\end{array}\right.$
式(11)中:C 2 为第2个约束条件,是关于集装器的约束条件;Xj 为X 的第j 列向量。
2)舱位与集装器的适配约束
定义集装器集合和舱位集合之间的适配约束矩阵为
(12) Y =$\left[\begin{array}{llll}{\mathit{y}}_{11}& {\mathit{y}}_{12}& \dots & {\mathit{y}}_{1{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\\ {\mathit{y}}_{21}& {\mathit{y}}_{22}& \dots & {\mathit{y}}_{2{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\\ \mathit{︙}& \mathit{︙}& & \mathit{︙}\\ {\mathit{y}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}1}& {\mathit{y}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}2}& \dots & {\mathit{y}}_{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}\end{array}\right]$
式(12)中:yij 为第i 个集装器与第j 个舱位之间的适配关系。
若两者适配,则yij 取值为1,否则取值为0。而且,若yij 为0,yij 必为0;若yij 为1,则yij 可以为0或1,两者之间的约束关系可表示为
(13) $\left\{\begin{array}{l}{\mathrm{C}}_{3}^{}:{\mathit{x}}_{\mathit{i}\mathit{j}}\le {\mathit{y}}_{\mathit{i}\mathit{j}}\\ \forall \mathrm{i}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}\}\\ \forall \mathrm{j}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}\}\end{array}\right.$
式(13)中:${\mathrm{C}}_{3}^{}$
3)最大业载限制
最大业载限制通常由最大零燃油质量、最大起飞质量和最大着陆质量确定,且这三种质量约束关系为
(14) C 4 :$\left\{\begin{array}{l}{\mathit{W}}_{\mathrm{T}\mathrm{O}}={\mathit{W}}_{\mathrm{O}\mathrm{E}}+{\mathit{W}}_{\mathrm{F}}+\stackrel{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}}{\sum _{\mathit{i}=1}}\stackrel{{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}{\sum _{\mathit{j}=1}}{\mathit{w}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}\mathit{j}}\le {\mathit{W}}_{\mathrm{M}\mathrm{T}\mathrm{O}}\\ {\mathit{W}}_{\mathrm{L}}={\mathit{W}}_{\mathrm{O}\mathrm{E}}+{\mathit{W}}_{\mathrm{R}\mathrm{F}}+\stackrel{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}}{\sum _{\mathit{i}=1}}\stackrel{{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}{\sum _{\mathit{j}=1}}{\mathit{w}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}\mathit{j}}\le {\mathit{W}}_{\mathrm{M}\mathrm{L}}\\ {\mathit{W}}_{\mathrm{Z}\mathrm{F}}={\mathit{W}}_{\mathrm{O}\mathrm{E}}+\stackrel{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}}{\sum _{\mathit{i}=1}}\stackrel{{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}{\sum _{\mathit{j}=1}}{\mathit{w}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}\mathit{j}}\le {\mathit{W}}_{\mathrm{M}\mathrm{Z}\mathrm{F}}\end{array}\right.$
式(14)中:C 4 为第4个约束条件;W TO 为飞机起飞质量;W OE 为飞机空机质量;W F 为起飞油量;W L 为飞机着陆质量;W RF 为飞机备用油量;W ZF 为飞机零燃油质量。
实际最大业载的限制为
(15) $\left\{\begin{array}{l}\stackrel{{\mathit{N}}_{\mathrm{U}\mathrm{L}\mathrm{D}}}{\sum _{\mathit{i}=1}}\stackrel{{\mathit{N}}_{\mathrm{P}\mathrm{O}\mathrm{S}}}{\sum _{\mathit{j}=1}}{\mathit{w}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}\mathit{j}}\le \stackrel{-}{\mathit{W}}\\ \stackrel{-}{\mathit{W}}=\mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{M}\mathrm{P}{\mathrm{L}}_{0},\mathit{M}\mathit{P}{\mathrm{L}}_{1},\mathit{M}\mathit{P}{\mathrm{L}}_{2},\mathit{M}\mathit{P}{\mathrm{L}}_{3})\end{array}\right.$
(16) $\left\{\begin{array}{l}\mathrm{M}\mathrm{P}{\mathrm{L}}_{0}={\mathit{W}}_{\mathrm{M}\mathrm{T}\mathrm{O}}-{\mathit{W}}_{\mathrm{O}\mathrm{E}}-{\mathit{W}}_{\mathit{L}}\\ \mathrm{M}\mathrm{P}{\mathrm{L}}_{1}={\mathit{W}}_{\mathrm{M}\mathrm{L}}-{\mathit{W}}_{\mathrm{O}\mathrm{E}}-{\mathit{W}}_{\mathrm{R}\mathrm{F}}\\ \mathrm{M}\mathrm{P}{\mathrm{L}}_{2}={\mathit{W}}_{\mathrm{M}\mathrm{Z}\mathrm{F}}-{\mathit{W}}_{\mathrm{O}\mathrm{E}}\end{array}\right.$
式(16)中:MPL0 、MPL1 、MPL2 为飞机规定载荷限制;MPL3 为航空公司视本公司的具体情况,并结合飞机载重平衡手册而设置的一个区别于MPL0 、MPL1 、MPL2 的最大业载。
4)结构约束
飞机的结构约束可分为舱位载荷约束、区域载荷约束、上下舱联合线载荷约束、前货舱累计载荷约束和后货舱累计载荷约束。
(1)舱位载荷约束:每个舱位放置集装器的质量都不能超过该舱位的最大载荷,即
(17) C 5 :W ULD Xj ≤${\mathit{W}}_{\mathit{j}}^{\mathrm{P}\mathrm{O}\mathrm{S}}$ , ∀ j ∈1,2,…,N POS
式(17)中:${\mathrm{C}}_{5}^{}$
(2)区域载荷约束:某区域中舱位上装载的集装器总质量不能超过该区域的最大载荷,即
(18) C 6 :SUM(W ULD X A k W A k k ∈{1,2,…,l }
式(18)中:C 6 为第6个约束条件,是关于区域载荷约束;W A k k 个区域限重;X A k k 个区域的所有舱位对应的决策变量子矩阵,定义G A k k 个区域中的舱位编号向量,可根据G A k X 中抽取相应的列向量来构建X A k
(3)上、下舱联合线载荷约束:根据飞机手册中提供的联合线载荷的区域划分情况,建立各线载荷区域的舱位编号向量,用G L t X 中抽取相应的列向量来构建X L t
(19) $\left\{\begin{array}{l}{\mathrm{C}}_{7}:{{\mathit{W}}^{\mathrm{L}}}_{\mathit{t}}{{\mathit{X}}^{\mathrm{L}}}_{\mathit{t}}\le \overline{{{\mathrm{W}}^{\mathit{L}}}_{\mathrm{t}}}\\ {{\mathit{W}}^{\mathrm{L}}}_{\mathit{t}}=[{{\mathrm{w}}^{\mathit{L}}}_{\mathrm{t}1}{{\mathrm{w}}^{\mathit{L}}}_{\mathrm{t}2}\dots {{\mathrm{w}}^{\mathit{L}}}_{\mathrm{t}{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}]\\ \mathit{t}\in \{\mathrm{1,2},\dots,\mathit{m}\}\end{array}\right.$
式(19)中:${\mathrm{C}}_{7}^{}$ W L t t 个线载荷区域的线载荷向量;m 为上、下舱联合线载荷约束区域的数量;$\overline{{{\mathit{W}}^{\mathrm{L}}}_{\mathit{t}}}$ t 个线载荷区域的限值。
(4)前货舱(从机头向后)累积载荷约束:同时考虑上、下货舱沿机体纵轴方向的累计载荷限制,具体定义为
(20) $\left\{\begin{array}{l}{\mathrm{C}}_{8}:{{\mathit{W}}^{\mathrm{F}}}_{\mathit{q}}{{\mathit{X}}^{\mathrm{F}}}_{\mathit{q}}\le {{\stackrel{-}{\mathrm{W}}}^{\mathit{F}}}_{\mathrm{q}}\\ {{\mathit{W}}^{\mathrm{F}}}_{\mathit{q}}=[{{\mathrm{w}}^{\mathit{F}}}_{\mathrm{q}1}{{\mathrm{w}}^{\mathit{F}}}_{\mathrm{q}2}\dots {{\mathrm{w}}^{\mathit{F}}}_{\mathrm{q}{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}]\\ \mathit{q}\in \{\mathrm{1,2},\dots,\mathit{n}\}\end{array}\right.$
式(20)中:C 8 为第8个约束条件,是关于前货舱(从机头向后)累积载荷约束;W F q q 区域的载荷向量;定义G F q q 个前货舱(从机头向后)累积载荷约束区域中的舱位编号向量,根据每个集装器实际装载舱位与前货舱累积载荷约束区域的对应关系,依次从X 中抽取相应的列向量来构建X F q n 为前货舱(从机头向后)累积载荷约束区域的数量。
(5)后货舱(从机尾向后)累积载荷约束:同时考虑上、下货舱沿机体纵轴方向的累计载荷限制,具体定义为
(21) $\left\{\begin{array}{l}{\mathrm{C}}_{9}:{{\mathit{W}}^{\mathrm{B}}}_{\mathit{r}}{{\mathit{X}}^{\mathrm{B}}}_{\mathit{r}}\le {{\stackrel{-}{\mathrm{W}}}^{\mathit{B}}}_{\mathrm{r}}\\ {{\mathit{W}}^{\mathrm{B}}}_{\mathit{r}}=[{{\mathrm{w}}^{\mathit{B}}}_{\mathrm{r}1}{{\mathrm{w}}^{\mathit{B}}}_{\mathrm{r}2}\dots {{\mathrm{w}}^{\mathit{B}}}_{\mathrm{r}{\mathrm{N}}_{\mathit{U}\mathit{L}\mathit{D}}}]\\ \mathit{r}\in \{\mathrm{1,2},\dots,\mathit{b}\}\end{array}\right.$
式(21)中:C 9 为第9个约束条件,关于后货舱(从机尾向后)累积载荷约束;W B r r 区域的载荷向量。定义G B r r 个后货舱(从机尾向后)累积载荷约束区域中的舱位编号向量,根据每个集装器实际装载舱位与后货舱累积载荷约束区域的对应关系,依次从X 中抽取相应的列向量来构建X B r
5)重心包络约束
(1)MAC重心的计算。
飞机的起飞重心平衡力臂、着陆平衡力臂和零燃油平衡力臂计算公式分别为
(22) $\left\{\begin{array}{l}{\mathrm{a}}_{{\mathrm{G}}_{\mathit{T}}}=\frac{{\mathrm{a}}_{\mathit{O}\mathit{E}}{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{a}}_{\mathit{F}}({\mathrm{W}}_{\mathit{F}}+{\mathrm{W}}_{\mathit{T}\mathit{F}})+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\left(\mathit{X}{\mathit{A}}_{\mathit{p}}\right)}{{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{W}}_{\mathit{T}\mathit{F}}+{\mathrm{W}}_{\mathit{F}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}}\\ {\mathrm{a}}_{{\mathrm{G}}_{\mathit{T}\mathit{O}}}=\frac{{\mathrm{a}}_{\mathit{O}\mathit{E}}{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{a}}_{\mathit{F}}{\mathrm{W}}_{\mathit{F}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\left(\mathit{X}{\mathit{A}}_{\mathit{p}}\right)}{{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{W}}_{\mathit{F}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}}\\ {\mathrm{a}}_{{\mathrm{G}}_{\mathit{L}}}=\frac{{\mathrm{a}}_{\mathit{O}\mathit{E}}{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{a}}_{\mathit{R}\mathit{F}}{\mathrm{W}}_{\mathit{R}\mathit{F}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\left(\mathit{X}{\mathit{A}}_{\mathit{p}}\right)}{{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathrm{W}}_{\mathit{R}\mathit{F}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}}\\ {\mathrm{a}}_{{\mathrm{G}}_{\mathit{Z}\mathit{F}}}=\frac{{\mathrm{a}}_{\mathit{O}\mathit{E}}{\mathrm{W}}_{\mathit{O}\mathit{E}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\left(\mathit{X}{\mathit{A}}_{\mathit{p}}\right)}{{\mathit{W}}_{\mathit{O}\mathit{E}}+{\mathit{W}}_{\mathit{U}\mathit{L}\mathit{D}}\mathit{X}}\\ {\mathit{A}}_{\mathit{p}}=[{\mathrm{a}}_{1}{\mathrm{a}}_{2}\dots {\mathrm{a}}_{{\mathrm{N}}_{\mathit{P}\mathit{O}\mathit{S}}}{]}^{\mathit{T}}\end{array}\right.$
式(22)中:W TF 为飞机从推出到起飞时因滑行所需消耗的燃油质量;a G T a G TO a G L a G ZF a OE 为操作空机质量对应的平衡力臂;a F 为燃油的整体平衡力臂;a RF 为备用燃油的平衡力臂;Ap 为第p 个舱位的平衡力臂列向量。
为了确保飞机的平衡与安全,重心必须保持在重心包络线的安全范围内。在实际使用过程中,飞机重心包络线横坐标通常表示为%MAC,航空中常用c 来表示平均气动弦长(meocn aerodynamic chord, MAC),即
(23) % MAC=$\frac{(\mathit{b}-\mathit{l})\times 100}{\stackrel{-}{\mathit{c}}}$
式(23)中: %MAC为飞机重心,单位是平均气动弦长的百分比;b 为集装器配载后的全机重心的平衡力臂长度;l 为平均气动弦长前缘的平衡力臂长度。将之前计算所得起飞重心平衡力臂、着陆平衡力臂、零燃油平衡力臂和飞行过程重心平衡力臂代入式(23)便可得到相应的飞机滑行重心G T 、起飞重心G TO 、着陆重心G L 、零燃油重心G ZF 和全机重心(G AIR ),其中G AIR =[${\mathit{G}}_{\mathrm{A}\mathrm{I}\mathrm{R}}^{1}$ ${\mathit{G}}_{\mathrm{A}\mathrm{I}\mathrm{R}}^{2}$ ${\mathit{G}}_{\mathrm{A}\mathrm{I}\mathrm{R}}^{\mathit{s}}$
(2)重心包络约束的定义。
为了便于阐述重心包络约束的构建方法,以B747-400F飞机的重心包络线图为例,详述该约束的定义过程如下。
如
图2 所示为B747-400F飞机重心可行性包络线图,
W T 、
W TO 、
W L 、
W ZF 四种不同质量下的重心边界围成了四个包络。
为将
图2 中的 B747-400F飞机重心可行性包络线转换为重心的偏移约束,按四种质量情况分别计算重心关于飞机总重。共有四段分段函数,以
W TO 包络线函数为例,表示为
(24) $\left\{\begin{array}{l}{\mathit{G}}_{\mathrm{T}\mathrm{O}\mathrm{W}}^{\mathrm{m}\mathrm{i}\mathrm{n}}=\left\{\begin{array}{l}13\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }160000\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<295742\mathrm{k}\mathrm{g}\\ 15.5\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }295742\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<352894\mathrm{k}\mathrm{g}\\ (0.0000637{\mathit{W}}_{\mathrm{T}\mathrm{O}}-7)\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }352894\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<396893\mathrm{k}\mathrm{g}\end{array}\right.\\ {\mathit{G}}_{\mathrm{T}\mathrm{O}\mathrm{W}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\left\{\begin{array}{l}(0.0000776{\mathit{W}}_{\mathrm{T}\mathrm{O}}+11.8847)\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }160000\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<272155\mathrm{k}\mathrm{g}\\ 33\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }272155\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<352894\mathrm{k}\mathrm{g}\\ (-0.000511{\mathit{W}}_{\mathrm{T}\mathrm{O}}+226.1)\mathrm{\%}\mathrm{M}\mathrm{A}\mathrm{C},\\ \mathit{ }\mathit{ }378296\mathrm{k}\mathrm{g}\le {\mathit{W}}_{\mathrm{T}\mathrm{O}}<396893\mathrm{k}\mathrm{g}\end{array}\right.\end{array}\right.$
式(24)中:W TO 包络线函数表达式中的质量区间按相应包络线上的拐点来确定,对应于每个质量区间,重心的函数表示式分为上限函数和下限函数。
将已知的实际W T 、W TO 、W L 、W ZF 分别代入以上不同质量的飞机重心分段函数表达式中,便可以计算出相应的重心上限和下限。
在此基础上,对起飞重心、着陆重心和零燃油重心的包络约束定义为
(25) C 10 :${\mathit{G}}_{\mathrm{T}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ ≤G T ≤${\mathit{G}}_{\mathrm{T}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$
(26) C 11 :${\mathit{G}}_{\mathrm{T}\mathrm{O}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ ≤G TO ≤${\mathit{G}}_{\mathrm{T}\mathrm{O}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$
(27) C 12 :${\mathit{G}}_{\mathrm{L}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ ≤G L ≤${\mathit{G}}_{\mathrm{L}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$
(28) C 13 :${\mathit{G}}_{\mathrm{Z}\mathrm{F}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ ≤G ZF ≤${\mathit{G}}_{\mathrm{Z}\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$
(29) $\begin{array}{l}{\mathit{C}}_{14}:{\mathit{G}}_{\mathrm{T}\mathrm{O}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathit{G}}_{\mathrm{A}\mathrm{I}\mathrm{R}}^{\mathit{s}}\le {\mathit{G}}_{\mathrm{T}\mathrm{O}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\\ \mathit{s}\in \{\mathrm{1,2},\dots,\mathit{c}\}\end{array}$
式中:C 10 为飞机滑行重心包络约束;C 11 为飞机起飞重心包络约束;C 12 为飞机着陆重心包络约束;C 13 为飞机零燃油重心包络约束;C 14 为飞机全机重心包络约束。
6)横向不平衡约束
假设装载的集装器均为标准集装器,且重心即为形心,则在水平面内,舱位的几何重心与集装器的重心重合,如
图3 所示为横向不平衡约束图。可以通过限制左右侧的质量差异来控制飞机横向的平衡。
(30) C 15 :|SUM(W ULD X R )-SUM(W ULD X L )|≤${\overline{\mathit{W}}}_{\mathrm{U}}$
式(30)中:C 15 为飞机横向不平衡约束;X L 、X R 分别为左侧舱位、右侧舱位与集装器之间的装载关系矩阵;${\overline{\mathit{W}}}_{\mathrm{U}}$
2 改进遗传算法设计
民航货机装载优化问题求解规模较大,约束条件复杂,从行业经验和航司实际运营来看,货运配载优化需要按优先级和时序构建分层优化目标。对于航空公司而言,在保证飞机安全的前提下,业载率最大化是其优先考虑的目标。而机动性和经济性两个目标的优化需要在业载率最高的装载方案基础上开展,即通过调整业载率最高的集装器集合在货舱中的位置来进行后两个目标的优化。对于以上问题,虽然可以使用非支配排序遗传算法Ⅱ(non-dominated sorting genetic algorithm Ⅱ,NSGA-Ⅱ)等多目标算法获得帕累托最优解集,但在航司实际决策时仍需明确各目标之间的主次关系,以此来确保优化策略与运营目标一致,这与NSGA-Ⅱ算法的特点是不相适应的。因此,使用改进后遗传算法来分阶段优化业载率等三个目标。遗传算法初始种群的优劣程度会影响优化的效率和质量,需要对每个步骤做相应的改进,包括初始种群的生成、适应度函数、交叉和变异等方面的改进,以提高算法的计算能力。
2.1 基于适配约束的初始种群生成方法
如
图4 所示为种群初始值生成流程图。为了提高优秀初始解的生成速度,分两种情况根据适配矩阵生成初始种群的个体,具体步骤如下。
1)N ULD ≥N POS
(1)定义初始种群为P 1 =[$\begin{array}{llll}{\mathit{X}}^{1}& {\mathit{X}}^{2}& \dots & {\mathit{X}}^{\mathit{I}}\end{array}$ I 为每代种群的个体数量,X 1 为初始种群中的第一个个体,其维度为N ULD ×N POS ,将适配矩阵赋值给X 1 ,即X 1 =Y 。
(2)从
X 1 的第
j (初始情况下
j =1)列中提取取值为1的元素,构建二维数组。该数组的第一列为所提取元素在适配矩阵中的行编号,第二列的元素全部为1,如
图4 所示。
(3)从上一步的二维数组中随机抽取一行,如[id ,1],其中id 为随机抽取元素在X 1 第j 列中的行编号。为了确保个体满足第1.3.3节约束条件式(11)C 1 和式(12)C 2 ,在X 1 的第j 列中,将第y i d j
(4)同样,为了满足第1.3.3节约束条件式(11)C 1 和式(12)C 2 ,在X 1 的第id 行中,将除y i d j X 1 的第j +1至N POS 列中,任意一列的全部元素均为零,则撤销在步骤(3)和步骤(4)对X 1 第j 列和第id 行的操作,返回步骤(3),重新进行随机抽取;否则,进入下一步。
(5)若j >N POS ,则终止流程,初始种群的第一个个体生成完毕。参照该方法依次生成初始种群的其他个体;否则,执行j =j +1,回到步骤(2)。
2)N ULD <N POS
(1)定义初始种群为${\mathrm{P}}_{}^{1}$ X 1 X 2 … XI ],其中I 为每代种群的个体数量,X 1 为初始种群中的第一个个体,其维度为N ULD ×N POS ,将适配矩阵赋值给X 1 ,即X 1 =Y 。
(2)从
X 1 的第
i (初始情况下,
i =1)行中提取取值为1的元素,步骤(2)参照第一种情况,构建二维数组如
图4 所示。该数组的第一列为所提取元素在适配矩阵中的列编号,第二列的元素全部为1。
(3)从上一步的二维数组中随机抽取一行,如[jd ,1],其中jd 为随机抽取元素X 1 在第i 列中的列编号。为了确保个体满足第1.3.3节约束条件式(11)C 1 和式(12)C 2 ,在X 1 的第i 行中,将第y ij d
(4)同样,为了满足前两项约束条件,在X 1 的第jd 行中,将除y ij d X 1 的第i +1至N ULD 列中,任意一列的全部元素均为零,则撤销在步骤(3)和步骤(4)对X 1 第i 行和第jd 列的操作,返回步骤(3),重新进行随机抽取;否则,进入下一步。
(5)若i >N ULD ,则终止流程,初始种群的第一个个体生成完毕。参照该方法依次生成初始种群的其他个体;否则,执行i =i +1,回到步骤(2)。
经过第一阶段的优化计算,可以得到满足所有约束条件且业载最大的最优解,其维度为1×N POS (N POS 为舱位数量)。根据得到的业载最大最优解,更新舱位-集装器适配矩阵,由上文初始解生成方法,得到第二阶段的初始解P 2 和第三阶段的初始解P 3 。
2.2 适应度函数和终止条件
由第1.3.3节关于约束条件定义,不同类型的约束涉及的限制条件不同。如适配约束的数量为多个,而最大业载限制仅有一个条件限制。由于各约束均为需同时满足的强约束,其重要性相同。然而,约束数量差异导致罚函数取值范围不一致,使适应度函数对不同约束的敏感程度不一,无法客观反映梯度变化,进而干扰算法寻优过程,造成收敛变慢和解的质量下降。
为此,在综合考虑目标函数取值范围和各类约束数量的基础上,对各类约束的罚函数进行归一化处理,以位置约束、适配约束和最大业载限制的适应度函数定义为例,有
(31) $\left\{\begin{array}{l}P_{1}^{j}\left(x_{i j}\right)=\left\{\begin{array}{ll}0, & \text { 不满足约束 } C_{1} \\1, & \text { 满足约束 } C_{1}\end{array}\right. \\P_{1}\left(x_{i j}\right)=\frac{\sum_{j=1}^{N_{\mathrm{POS}}} P_{1}^{j}\left(x_{i j}\right)}{N_{\mathrm{POS}}}, \quad j \in\left\{1,2, \cdots, N_{\mathrm{POS}}\right\}\end{array}\right.$
(32) $\left\{\begin{array}{l}P_{2}^{j}\left(x_{i j}\right)=\left\{\begin{array}{ll}0, & \text { 不满足约束 } C_{2} \\1, & \text { 满足约束 } C_{2}\end{array}\right. \\P_{2}\left(x_{i j}\right)=\frac{\sum_{j=1}^{N_{\mathrm{ULD}}} P_{2}^{i}\left(x_{i j}\right)}{N_{\mathrm{IIID}}}, \quad j \in\left\{1,2, \cdots, N_{\mathrm{ULD}}\right\}\end{array}\right.$
(33) $P_{3}\left(x_{i j}\right)=\left\{\begin{array}{ll}0, & \text { 满足约束 } C_{3} \\\frac{-\mathrm{SUM}\left\{\mathrm{DLAG}\left[\boldsymbol{X}^{\mathrm{T}}(\boldsymbol{Y}-\boldsymbol{X})\right]\right\}}{N_{\mathrm{POS}}}, & \text { 不满足约束 } C_{3}\end{array}\right.$
(34) $P_{4}\left(x_{i j}\right)=\left\{\begin{array}{ll}0, & \text { 满足约束 } C_{4} \\1, & \text { 不满足约束 } C_{4}\end{array}\right.$
以式(31) ~式(34)为例,设置与约束条件一一对应的15个罚函数公式。
综上所述,结合前文对各阶段目标函数的定义,分别定义三个阶段的适应度函数为
(35) $\begin{array}{c}{\mathit{F}}_{\mathit{n}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)=\mathit{f}\mathit{\text{'}}{}_{\mathit{n}}\left({\mathit{x}}_{\mathit{i}\mathit{j}}\right)+\stackrel{15}{\sum _{\mathit{i}=1}}{\mathit{P}}_{\mathit{n}},\mathit{ }\mathit{n}=\mathrm{1,2},3\end{array}$
式(35)中:Fn (xij )分别为三个阶段的适应度函数。
为了防止陷入局部最优,且提高优化的稳定性和最优解的质量。同时采用适应度无显著变化和最大迭代次数作为优化的终止条件,其中,适应度变化的表征参数为连续代最优适应度值的平均变化量,计算方法为
(36) $\left\{\begin{array}{l}{\mathit{T}}_{\mathit{e}}^{{\mathit{N}}_{\mathit{T}}-1}={\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}}\left(\mathrm{X}\right)-{\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}-1}\left(\mathrm{X}\right)\\ {\mathit{T}}_{\mathit{e}}^{{\mathit{N}}_{\mathit{T}}-2}={\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}-1}\left(\mathrm{X}\right)-{\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}-2}\left(\mathrm{X}\right)\\ \mathit{ }\mathit{ }\mathit{ }\mathit{︙}\\ {\mathit{T}}_{\mathit{e}}^{1}={\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}-{\mathit{N}}_{\mathit{T}}+1}\left(\mathrm{X}\right)-{\mathit{F}}_{\mathit{e},\mathrm{o}\mathrm{p}\mathrm{t}}^{\mathit{g}-\mathit{N}}\left(\mathrm{X}\right)\\ {\mathit{T}}_{\mathit{e}}=\frac{\stackrel{{\mathit{N}}_{\mathit{T}}}{\sum _{\mathit{i}=1}}\left|{\mathit{T}}_{\mathit{e}}^{\mathit{i}}\right|}{{\mathit{N}}_{\mathit{T}}-1}\\ \mathit{g}\in \{\mathrm{1,2},\dots,{\mathit{N}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\},\mathit{e}\in \{\mathrm{1,2},3\}\end{array}\right.$
式(36)中:g 、e 、N max 分别为当前的迭代代数、优化阶段编号和最大迭代次数;${\mathit{F}}_{\mathrm{e},\mathit{o}\mathit{p}\mathit{t}}^{\mathit{g}}$ X )为第e 优化阶段的第g 代的最优适应度值。
由此定义优化终止条件的数学描述为
式(37)中:te 为连续代最优适应度值的平均变化量阈值,结合目标函数和罚函数的取值范围,一般将该参数设置为10-4 ~ 10-5 。而将NT 设定为N max 的10%~20%。
2.3 选择算子、变异与交叉
为了提高算法的全局收敛性和计算效率,避免由于遗传操作破坏进化过程中适应度最高的染色体,本文改进的遗传算法采用最优保存策略对模型优胜劣汰操作。该策略首先对当前群体按照适应度进行排序,将适应度最好的前10%个体直接保留到下一代,然后随机生成10%的新个体与剩余90%个体组合形成临时种群,对该临时种群按适应度值从优到劣排序后选取前90%的个体进行交叉与变异运算,最终子代群体由直接保留的前10%精英个体和经过交叉变异操作后产生的个体共同组成,从而在保证种群多样性的同时确保优秀基因。
决策变量的行、列向量元素之间分别存在位置约束。因此,采用现有遗传算法中的随机变异进行个体基因的直接改变,可能会导致变异后个体质量变差,致使种群的优质个体变少,导致收敛速度慢,最优解质量下降等问题。
针对该问题,提出基于位置约束的变异方法,该方法遵循位置约束来确定变异基因的数值,可有效防止变异基因不满足位置约束的情况出现,具体操作如下。
1)N ULD ≥N POS
(1)将0-1决策变量矩阵转化为一维数组M ,该向量的元素数值为ULD编号,元素的顺序为舱位编号。
(2)进行基因变异,即一维数组中的ULD编号突变时需要遵循位置约束,当变异基因是一维数组中顺序为j 的元素(即改变舱位j 上放置的ULD)时,该元素的取值是一维数组中未出现的ULD编号,满足i ∉M ,i ∈{1,2,…,N ULD },这样可确保变异运算满足约束C 1 和C 2 。
(3)在i ∉M ,i ∈{1,2,…,N ULD }中随机选取一个数值作为变异基因的取值。
2)N ULD <N POS
除基因变异中,变异基因为M 中顺序为i 的元素(即改变第i 个ULD的放置舱位)外,其余步骤与N ULD ≥N POS 步骤保持不变。
为了避免交叉过程中出现重复基因,选择基于位置交叉的方法,这样可使新染色体适应度更高,且可以避免交叉后染色体中存在重复基因。具体操作过程如下。
(1)随机生成一个0-1的浮点数,若该数大于给定的交叉概率,则父代染色体与配对染色体进行基于位置的交叉。
(2)随机选择一个父代中的若干个基因作为非交叉基因,位置可不连续。
(3)生成一个子代,在子代中,非交叉位置的基因与其父代相同。
(4)将配对父代中与非交叉基因相同的元素保留,其他元素按原始顺序移至子代的交叉基因位置。
3 方法验证
3.1 三阶段递进优化策略模型的有效性验证
装载集装器使用的货机为B747-400F,使用的数据均来源于B747-400F质量和平衡控制以及装载手册
[24 ] 。B747-400F质量限制数据如
表1 所示。上下舱联合线载荷限值如
表2 所示。前后货舱累积载荷限值如
表3 所示。
Limbourg等
[5 ] 关于集装器装载问题的研究成果是航空货运装载优化领域的经典文献之一,其提出的优化方法在航空公司得到了成功应用。为验证本文所提模型的有效性及算法的可行性,采用文献[
5 ]中的标准测试数据集进行装载优化实验,并与文献[
5 ]的优化结果作对比,以验证模型的有效性和算法可行性。
文献[
5 ]与本文模型存在部分差异。文献[
5 ]中通过对转动惯量的优化减少货物在重心范围内的堆积,提高飞机的机动性和经济性,仅考虑将零燃油重心偏差作为一个约束来保证飞机的安全性。而本文模型综合实际装载需求,不仅对俯仰转动惯量进行优化,同时优化飞机最大业载以及飞机起飞重心与目标重心的偏差以实现对飞机安全性、机动性和经济性的优化。相对于文献[
5 ],本文中增加了适配约束、舱位和区域质量约束、最大业载限制,特别是在重心约束中,同时对起飞重心、着陆重心和零燃油重心进行约束。在模型求解方法上,文献[
5 ]使用的是经典分支切割法,本文中采用的是改进遗传算法。
为了便于进行对比验证,保证优化目标一致,在相同约束的基础上加入本文约束,使用同一组集装器数据。在验证过程中,本模型为三阶段递进优化,最终验证说明在业载最大的情况下飞机的俯仰转动惯量最小,飞机的起飞重心与目标重心的偏差最小。
由于文献[
5 ]所验证的该算例所选飞机为B747-400飞机,该货舱拥有42个装载舱位,案例1待装载集装器数量为26个,案例2待装载集装器数量为42个。B747-400飞机上下货舱的俯仰图如
图6 所示,待装载集装器信息及与集装器的适配情况如
表4 ~
表6 所示。
将本模型中第三阶段优化的权重系数
ρ 设置为1,最终目标函数与文献[
5 ]一样仅优化最小俯仰转动惯量,并将重心约束中的零燃油重心约束进一步缩紧,设置为与重心28%MAC的偏差小于等于0.01%。 设置起飞油量为101 500 kg,航程油量为82 300 kg。
下面对算例应用改进遗传算法进行求解,程序主要变量设置如下:种群规模为100,最大迭代次数3 000,交叉概率0.9,变异概率0.1。得到三阶段递进遗传算法的收敛图和装载方案如
图7 、
图8 所示。
为了验证算法的可行性,通过上文的约束条件对本次装载结果进行验证。
(1)位置约束。如
图8 所示为最优装载方案对比图。灰色部分为与文献[
5 ]装载方案不同的舱位。通过装载时产生的适配矩阵
Y ,对装载方案验证是否满足适配约束。根据
表7 、
表8 所示适配关系
yij 与
xij 对比可知,
yij =0时,
xij =0;
yij =1时,
xij =0或1,满足适配约束。
(2)载荷约束。区域载荷约束验证结果如
图9 所示。经验证文献[
5 ]案例1、案例2和本文优化结果主货舱和下货舱装载后的实际区域载荷均未超过每个区域的最大载荷限值,满足约束。
(3)上下舱联合线载荷约束。上下舱联合线载荷约束验证结果如
图10 所示,经验证文献[
5 ]案例1、案例2结果和本文优化结果实际装载的线载荷均未超过上下舱联合线载荷的最大限值,满足约束。
(4)前后货舱累积载荷约束。前后累积载荷约束验证结果如
图11 所示,在前货舱的限制区域,文献[
5 ]案例1、案例2和本文结果实际装载后的载荷未超过前货舱累积载荷的最大限值;在后货舱的限制区域,实际装载的累积载荷未超过后货舱累积载荷的最大限值,满足约束。
5)重心范围约束
文献[
5 ]案例1、案例2重心包络验证结果如
图12 、
图13 所示。文献[
5 ]对
G ZF 设置一个扩展区域,装载的
G ZF 可以保持在限制区域内。案例1优化后得到的
W ZF 的质量为224 237 kg,
G ZF 为28.000%MAC,位于拓展区域内,满足约束。文献[
5 ]优化后得到的
G TO 和
G L 重心分别为25.11%MAC、30.72%MAC,优化后得到的
W TO 和
W L 的质量为325 737、243 437 kg。均在各包络线范围内,满足约束。
案例2优化后得到
W ZF 的质量为279 539 kg,
G ZF 为27.999%MAC,位于拓展区域内,满足约束。文献[
5 ]优化后得到的
G TO 和
G L 重心分别为25.52%MAC、30.20%MAC,优化后得到的
W TO 和
W L 的质量为381 039、298 739 kg。均在各包络线范围内,满足约束。
6)横向不平衡约束
由
表9 可知,文献[
5 ]案例1、文献[
5 ]案例2实际装载后左右装载质量的差异分别为15 940、2 957 kg,限值为33 548 kg,满足约束。
如
表9 所示为文献[
5 ]与本文装载结果对比, 以俯仰转动惯量最小化为优化目标,案例1本文得出的转动惯量最优值能达到1.36×10
10 kg·in
2 ,与文献[
5 ]达到同等量级,且比文献[
5 ]多优化21.8%,零燃油质量下的重心与目标值偏差均能低于0.01%,由于针对位置、重心和位置方面添加了更严格的约束经过三阶段优化后,得到的重心偏差比文献[
5 ]得出的偏差减少0.007,偏差更小。
案例2本文得出转动惯量最优能达到2.98×10
10 kg·in
2 ,与文献[
5 ]达到同等量级,且比文献[
5 ]多优化了3.87%,零燃油质量下的重心与目标重心的偏差均能低于0.01%,由于本文针对位置、重心和质量方面添加了更严格的约束,经过三阶段优化后,得到的重心偏差比文献[
5 ]得出的偏差减少0.002,偏差更小。综合以上案例1、案例2可反映出本文的优化模型优化效果更好。
为了评价最优解适应度函数值相对于初始解适应度函数值的提升,按案例1、案例2各阶段的单次运行分别计算提升情况,统计适应度函数值的最大提升百分比、每个案例各优化阶段的平均提升百分比和变异系数,结果如
表10 、
表11 所示。
如
表10 、
表11 所示,为验证本文模型及算法效果,引入上述指标用于量化分析适应度函数值的变化幅值和离散程度。适应度函数值变化百分比、平均变化百分比以及变异系数数值越小,表示初始解质量越高,且算法稳定性越高。
由
表10 、
表11 可见,案例1、案例2在第一、二、三阶段中均保持了较高的优化效果,同时适应度函数值最高提升百分比与平均提升百分比差值较小,说明良好的初始解有效提升后续阶段优化结果。
由
表10 、
表11 可知,案例1、案例2变异系数均在10%以下,表示优化后的结果离散程度小,数据较为集中,印证优化后的适应度函数值具有良好的稳定性,反映出本文初始解质量高、模型有效。
综上所述,通过与文献[
5 ]的模型求解对比,验证了本文模型有效性和算法可行性,求解结果属于全局最优。
3.2 模型、算法的泛化性和有效性验证
某快递公司向合作的航空公司预定货机来装载货物,飞机单程运输,从A市到B市,飞机型号为B747-400F,如
图14 、
图15 所示分别为B747-400F飞机及货舱图、ULD类型图。以某货单的集装器为例,随机抽取该航空公司某一日B747-400F的运输数据,起飞油量101 500 kg,航程油量82 300 kg。取其中大、中、小三个不同规模的集装器作为案例,来验证本文模型与算法处理集装器装机问题的能力,每个集装器的类型均适用于B747-400F的装载。
本算例要求分别将三批不同规模集装器装载到货舱中,进行三次优化,分别是业载最大实验、俯仰转动惯量最大化实验、俯仰转动惯量最小和重心偏移最小实验,然后分析对飞机安全性、机动性和经济性的影响。本文所计算案例时,设置目标重心%MAC=26。程序主要相关变量设置:种群规模为100。为提高运行结果准确性,最大迭代次数提高至3 000,交叉概率为0.9,变异概率为0.1,每个算例均独立运行100次,最后取平均时间。
3.2.1 结果与分析
为了进一步说明本文所提出的模型和算法的有效性,选取不同规模的标准集装器进行装载优化。案例1~案例3三个目标优化的结果如
表12 ~
表14 所示,分别表示业载最大实验优化结果、俯仰转动惯量最大实验优化结果、俯仰转动惯量最小、重心偏移最小实验优化结果。
根据
表12 ~
表14 所示的案例1~案例3三阶段优化结果,采用递进式三阶段优化策略,引入表征飞机机动性能、经济性能和安全性能的关键指标——俯仰转动惯量、载重率及重心位置,对优化效果进行分析。
由
表12 ~
表14 可知,案例1~案例3在第一阶段优化中,载重率最大可达到100%,反映了以载质量最大为优化目标取得了理想效果。案例1~案例3第二阶段得到的最大俯仰转动惯量是在三阶段中计算俯仰转动惯量最大值时进行归一化处理。案例1、2、3的第二阶段优化结果,相较于第一阶段和第三阶段,俯仰转动惯量最大增幅达到43.21%。案例1、2、3第三阶段优化中,实际重心位置与目标重心位置(26%MAC)的偏差为0,实现了最优配载。在此基础上,俯仰转动惯量相较于第二阶段均有所降低,最大降幅达43.21%。
综上所述,经本文提出的优化模型和算法求解,获得了业载最大、俯仰转动惯量最小化和重心偏差最小化的配载方案。为验证优化后配载方案的正确性,需进一步检验其是否满足所有约束条件。
(1)位置约束:如
图16 所示为三组实验舱位装载图,通过装载时产生的适配矩阵
Y ,对装载方案验证是否满足适配约束。根据适配关系对比可知,
yij =0时,
xij =0;
yij =1时,
xij =0或1,不发生冲突,满足适配约束。
(2)区域载荷约束、上下舱联合载荷约束、前后货舱累积载荷约束、重心包络约束:由
图17(a) 、
图17(b) 可知,三组实验实际装载的区域载荷均在区域载荷最大限值,满足约束;由
图17(c) 可知,三组实验实际装载的线载荷未超过上下舱联合线载荷限值,满足约束;由
图17(d) 、
图17(e) 可知,三组实验的实际装载计算得出前后累积载荷均不超过限值,满足约束;由
图18(f) 可知,三组实验得出的
G TO 、
G L 和
G ZF 三种重心均在
W MTO 、
W ML 和
W MZF 与%MAC限值组成包络线内,满足约束。
(3)横向不平衡约束:由
图16 所示装载图可知,三组装载实验质量差异分别为6 335、242、3 574 kg,不超过12 648 kg,满足约束。
3.2.2 算例结果分析
1)初始解质量
为了证明初始解质量提升在提高遗传算法优化过程稳定性和效率方面的效果,设计了大、中、小三种规模的案例,就单个案例分三个优化阶段分别实施了100次优化,每个案例各阶段的单次优化最优解对应的适应度函数值初始值和优化值对比如
图18 所示。
为了评价最优解适应度函数值相对于初始解适应度函数值的提升,可以按每个案例各阶段的单次运行来分别计算提升情况,统计适应度函数值的最大提升百分比、每个案例各优化阶段的平均提升百分比和变异次数。如
表15 ~
表17 所示分别为案例1~案例3适应度函数值最大提升百分比以及案例1~案例3各优化阶段的平均提升百分比。
综合
图18 以及
表15 ~
表17 可知,在ULD数量与舱位数量不同的三种案例中,各案例中适应度函数值最高与平均提升百分比差异均较小(<15%),且最终变异系数均低于10%,表明优化结果离散度小,验证了算法的稳定性与收敛性。
案例1中因ULD数量少于舱位,第一阶段优化空间受限,第二、三阶段适应度函数值分别提升41.32%与23.31%。案例2中ULD略大于舱位,三阶段平均提升幅度依次为63.47%、45.69%、21.20%,
案例3中ULD远大于舱位,第一阶段平均提升达99.96%,后续阶段仍保持较高改进幅度,表明算法在适度冗余条件下表现更优。
综合以上结果可知,高质量的初始解不仅在第一阶段表现优异,还在后续阶段提供了更优的初始解,形成了“优质初始解→优化效果显著→更优中间解→持续优化改善”的正向循环。
2)算法有效性
如
图19 所示为案例1~案例3三阶段迭代曲线,可见案例1~案例3通过改进的遗传算法的收敛图收敛效果较好,可在200代内迅速收敛。
3)鲁棒性
三组案例的集装器类型与质量均不同,具有强异构性,但根据实验结果可以看出,无论集装器数量有多少,通过建立的数学模型,迭代过程中均能在200代前得到案例结果,均能较快地得到实验结果,模型具有较强的鲁棒性。
4 结论
针对航空标准集装器装载优化研究中目标单一等问题,设计了面向航空公司多重运营目标的三阶段递进配载优化方法,依次以业载最大、俯仰转动惯量最大、俯仰转动惯量最小和重心偏移最小为目标函数,使用改进的遗传算法对集装器的配载方案进行了求解,并以B747-400F飞机的集装器配载实例验证了策略、模型与算法的可行性与有效性,得出以下结论。
(1)三阶段递进优化策略可实现航司飞行安全性、收益性、机动性和经济性等目标的归一量化和优先级排序。在不同规模的案例验证中,该策略均能在保证飞行安全的前提下,将收益性、机动性与经济性先后优化至最佳,这不仅保证了优化过程与航司运控目标优先级的一致性,还能确保优化结果符合航司运营目标,具备良好的工程适用性与优化效率。
(2)利用二进制适配矩阵定义所有类型ULD与具体机型舱位之间的适配关系,在此基础上设计了基于随机抽样的初始解生成器。算例结果证明,基于适配关系和随机抽样的初始解生成器在严格保证舱位约束的基础上提高了初始解的质量,进而提升了遗传算法的全局寻优能力、收敛速度和稳定性。
(3)所设计的罚函数归一化方法可以均衡各类约束权重、客观定义各类约束对适应度函数的影响;引入最优保存选择策略、PBX交叉算子和基于位置约束的变异方法,可有效解决传统遗传算法在处理位置约束优化问题时的不足。从优化效果来看,改进遗传算法均在200代内稳定收敛,在保证解的质量的同时,显著提高了收敛速度和计算效率。
后续的研究可以进一步贴近航司实际需求,考虑实际重心位置的奇异集装器的装载,根据航班类型、航线特点和运营优先级动态调整安全性、机动性和经济性的权重,实现实时优化系统等,以及采用深度强化学习方法来开展集装器配载优化。
PDF
(21687KB)
