基于轨道几何状态的大机作业决策方法

摘要

大型养路机械作业主要包括捣固作业和打磨作业，是维持线路良好几何状态的重要作业方式。为了制定合理的捣固与打磨作业计划，保障线路能够安全运营，以年度几何不平顺平均值最优为目标，现场作业能力为约束，建立了年度捣固与打磨作业计划编制模型，并且运用改进的枚举法对模型进行求解。该模型被应用于京沪铁路进行验证，结果表明：该模型建立的大机作业计划能够提高轨道几何状态的改善率，降低年度几何不平顺平均值，为现场捣固和打磨作业的计划编制提供参考。

关键词

有砟轨道; 几何不平顺; 计划编制模型; 捣固作业; 打磨作业

近年来，我国铁路运输需求不断增大，轨道恶化速率加快，这给线路的养护维修工作带来了很大的困难^［1］。为了使得养修工作更加高效，往往需要根据线路的几何状态进行养修策略的制定，实现周期修到状态修的转变^［2］。捣固作业和打磨作业是改善有砟轨道线路状态的大机作业方式，捣固作业可以大幅度降低线路的几何不平顺值，打磨作业可以缓解线路几何不平顺值的恶化速率，从而间接改善线路的几何状态^［3］。合理地安排捣固与打磨作业可以有效地控制有砟轨道几何状态的恶化，保证线路的运营状态。因此，根据线路轨道状态建立一个年度有砟轨道捣固与打磨作业养护维修计划十分重要。

为此，国内外学者纷纷展开了针对有砟轨道养护维修计划编制模型的研究：Consilvio等^［4］分析了现场轨道线路的故障概率，在保证列车安全运行的前提下，实现最低的养护维修成本；Macedo等^［5］在考虑养修资源的同时以养护维修成本最低为目标，提出混合整数线性规划模型来进行捣固计划编制，并且运用变邻域搜索算法来求解模型；周宇^［6］建立了年末轨道状态最优综合养护计划模型，并运用了遗传算法进行求解，提高了计算效率，实现了线路养护维修的辅助决策；郭然^［7］以最低维护费用和最低大型养路机械走行费用为目标函数，针对高速铁路建立了养护维修策略最优化模型，建立了高速铁路月度养修计划模型；陶竑宇等^［8］以线路单元区段作业质量为研究对象，建立了预防性养护维修轨道管理方法模型，实现了预防性轨道养护维修计划的编制；邱俊兴等^［9］提出了综合评价指数的捣固决策方法，考虑了列车运行品质和捣固作业成本两方面进行捣固决策；徐伟昌^［10］分析了大型养路机械捣固作业轨道质量评价指数均值及方差随时间变化规律，给出了“线路维修服务水平”的概念，并基于此确定了捣固作业维修周期。然而上述学者的制定的养护维修计算模型仍存在以下不足：不能根据轨道几何状态编制合理的养修计划；仅考虑年末轨道几何状态最优而不能考虑全年轨道几何状态均值最优；仅考虑了捣固作业对轨道几何状态的影响，几乎没有文献同时进行捣固与打磨作业计划编制，作业计划编制不完善。

在上述研究的基础上，本文拟进行考虑年度轨道几何状态均值最优的有砟轨道捣固与打磨作业决策研究，主要根据相关文献^［11］中已建立的有砟轨道几何不平顺预测模型来建立预期目标函数，做到同时考虑捣固与打磨作业；进而确定符合现场作业的相关约束条件；并用改进的枚举法进行求解。同时利用建立好的模型进行了京沪线线路的实例分析：证实了模型的可靠性以及求解方法的有效性，可以为现场捣固与打磨作业安排提供参考。

1 有砟轨道养护维修计划编制模型

目前，我国现有的养护维修计划编制方法大都以周期修为主，忽略了线路不平顺的劣化特点^［12］。同时，由于全年的天窗时间是有限的，难以在有限的天窗时间合理地安排养修项目和作业区段。因此，应当结合线路状态和实际作业能力来安排养修计划，实现以状态修为主的养护维修模式。

作业计划编制模型设计的思路为：首先对线路进行作业区段的划分以及根据天窗时间确定大机作业的具体时间；其次根据轨检车数据以及各作业区段轨道几何状态恶化特征与恢复特征，形成恶化函数、捣固与打磨作业恢复函数，并组装形成恶化矩阵；然后根据恶化函数、恢复函数以及恶化矩阵推导出以年度轨道几何状态均值最低的目标函数，同时考虑现场适用性给出约束函数；最后运用改进枚举法进行求解，得到各作业区段最适宜的捣固与打磨作业天窗时间。该模型的优点是考虑到了轨道几何状态的发展，依据轨道几何状态的预测编制年度均值最低的目标函数；既考虑了捣固作业又考虑了打磨作业对轨道状态的影响，对现有研究进行了补充，使得作业计划编制更加完善。具体作业计划编制模型的构建流程如图1所示。图中，TQI表示轨道质量指数（track quality index）。

图1  有砟轨道养护维修计划编制模型流程图

Fig. 1  Flow chart of preparation model of ballast track maintenance plan

1.1　模型基本假设

为了使养护维修计划编制模型符合实际情况，并考虑到运算量和求解难度，对此模型做出如下假设：

（1）模型的基本作业区间长度为4 km。经过对作业现场一个天窗时间捣固与打磨作业可完成的区段长度进行调研，单次捣固作业平均里程为4 km，单次打磨作业平均里程为8 km；

（2）大型养路机械单次捣固作业需占用一个天窗时间，可完成一个作业区段；单次打磨作业同样需占用一个天窗时间，可完成两个连续作业区段；

（3）模型中捣固和打磨作业仅考虑线路捣固和线路打磨，不考虑道岔捣固和道岔打磨；

（4）大型养路机械捣固作业在一个天窗时间内，只能在一个作业区域工作，打磨作业只能在两个连续作业区域工作；

（5）每个基本作业区段一年最多只进行一次捣固作业和打磨作业；

（6）只考虑钢轨打磨、捣固两种作业类型对轨道不平顺值的影响，不考虑人工作业；

（7）在同一作业区段内，恶化函数与恢复函数的相关参数由历年来轨检车数据决定，不考虑其他因素；

（8）线路捣固作业仅改变无轨距的TQI值大小，而线路打磨作业仅改变无轨距的TQI值恶化速率^［13］；

（9）模型计算步长为一个天窗时间间隔；

（10）不考虑轨道检测、诊断的时间影响。

1.2　模型基本参数定义

本文基于文献［14］制定的轨道几何状态预测模型中相关理论阐述：不同区段线路的恶化特点不同，在恶化模型中应针对不同的作业区段设置不同的恶化函数，并形成恶化矩阵。对模型的基本参数定义如下：

（1）大型养路机械作业时间矩阵（决策变量）$\mathit{X}$：$\mathit{X}$为$m\times \mathrm{2}$的矩阵，$m$为基本作业区段个数，$\mathit{X}$的取值范围为$\left[\mathrm{0},m_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]$，$m_{\mathrm{m}\mathrm{a}\mathrm{x}}$为全年天窗个数。$\mathit{X}$的第一列表示捣固作业的天窗编号，第二列表示打磨作业的天窗编号。例如，若$\mathit{X}(\mathrm{7,2})=\mathrm{45}$，表示第7个基本作业区段在45号天窗进行打磨作业。

（2）修正后的恶化矩阵${\mathit{P}}_{\mathrm{m}}\text{'}$：${\mathit{P}}_{\mathrm{m}}\text{'}$为$n\times \mathrm{3}$的矩阵。${\mathit{P}}_{\mathrm{m}}\text{'}$的第一列为不平顺值，第二列为打磨前一个天窗时间间隔的恶化率，第三列为打磨后一个天窗时间间隔的恶化率，下标表示基本作业区段编号，即对于每一个基本作业区段来说，都有其独特的恶化矩阵。

（3）无轨距几何不平顺标准差${\sigma }_{i,j}^{\mathrm{N}}$，mm：因捣固与打磨作业对轨距不平顺影响不大^［14］，故考虑无轨距几何不平顺值，包括：左高低、右高低、左轨向、右轨向、水平、三角坑6个单项之和，用${\sigma }^{\mathrm{N}}$来指代。下标$i$和$j$分别表示地点和时间。例如${\sigma }_{\mathrm{6,27}}^{\mathrm{N}}$表示第6个作业区段在第72个天窗时的无轨距不平顺值；

（4）基本作业区段编号$m$和天窗时间编号$l$。

（5）恶化率$k$：用来描述线路几何不平顺的恶化情况，取值范围为［0，+∞），与时间、线路、作业情况有关。

（6）恢复率$J$：用来描述大型养路机械线路捣固作业后线路的恢复情况，取值范围为［0，1）。

1.3　目标函数推导

根据文献［14］中的恶化率和恢复率函数推导，恶化率$k$可表示为

$$k=\frac{\left({\sigma }_m^{\mathrm{N}}-C^{\mathrm{\text{'}}}\right)·{\mathrm{e}}^{b^{\mathrm{\text{'}}}·T\left(l\right)}+C^{\mathrm{\text{'}}}}{{\sigma }_m^N}\left(\mathrm{1}-\frac{p}{T}\right)+\frac{\left({\sigma }_m^{\mathrm{N}}-C\right)·{\mathrm{e}}^{b·T\left(l\right)}+C}{{\sigma }_m^{\mathrm{N}}}·\frac{p}{T}-\mathrm{1}$$

（1）

式中：${\sigma }_m^{\mathrm{N}}$表示区段$m$的无轨距不平顺值，mm；$T\left(l\right)$表示天窗l时的天窗间隔，d；$b$和$C$为待定系数；$p$为上次线路打磨作业距此时间区段的时间，d；$T$为线路打磨周期，d。

捣固作业恢复率函数可表示为

$$J\left({\sigma }^{\mathrm{N}},m\right)=D_m·\mathrm{l}\mathrm{n}\left({\sigma }^{\mathrm{N}}\right)+F_m$$

（2）

其中：参数$D_m$和$F_m$为待定系数，且与区段编号$m$相关。

无轨距TQI初始值可表示为${\sigma }_{m,\mathrm{1}}^N$。若在第一个天窗时间对该区段进行线路捣固作业，则第二个天窗时间的无轨距TQI值可表示为

$${\sigma }_{m,\mathrm{2}}^{\mathrm{N}}={\sigma }_{m,\mathrm{1}}^{\mathrm{N}}\mathrm{\text{'}}·\left(\mathrm{1}+k\right)$$

（3）

式中：${\sigma }_{m,\mathrm{1}}^{\mathrm{N}}\mathrm{\text{'}}={\sigma }_{m,\mathrm{1}}^{\mathrm{N}}·(\mathrm{1}-J)$，表示完成线路捣固作业后的无轨距TQI值（下同）。若不在第一个天窗时间进行线路捣固作业，则二个天窗时间的无轨距TQI值可表示为

$${\sigma }_{m,\mathrm{2}}^{\mathrm{N}}={\sigma }_{m,\mathrm{1}}^{\mathrm{N}}·\left(\mathrm{1}+k\right)$$

（4）

联立式（1）—（4）可得：

$${\sigma }_{m,\mathrm{2}}^{\mathrm{N}}={\sigma }_{m,\mathrm{1}}^{\mathrm{N}}·{\left(\mathrm{1}+k\right)}^{\mathrm{1}-Y\left(m,\mathrm{1}\right)}·{\left[\frac{{\sigma }_{m,\mathrm{1}}^{\mathrm{N}}\mathrm{\text{'}}}{{\sigma }_{m,\mathrm{1}}^{\mathrm{N}}}\left(\mathrm{1}+k\right)\right]}^{Y\left(m,\mathrm{1}\right)}$$

（5）

式中：$\mathit{Y}$为$m_{\mathrm{m}\mathrm{a}\mathrm{x}}\times l_{\mathrm{m}\mathrm{a}\mathrm{x}}$的二维辅助决策矩阵，$Y(i,j)$表示矩阵$\mathit{Y}$第$i$行第$j$列的取值。矩阵的取值规则如式（6）所示：

$$\begin{array}{l}Y\left(i,j\right)=\\ \left\{\begin{array}{c}\mathrm{1}\mathrm{区域}i\mathrm{在}j\mathrm{天窗}\mathrm{时间}\mathrm{进行}\mathrm{捣固}\mathrm{作业}\\ \mathrm{0}\mathrm{区域}i\mathrm{不在}j\mathrm{天窗}\mathrm{时间}\mathrm{进行}\mathrm{捣固}\mathrm{作业}\end{array}\right.\end{array}$$

（6）

决策变量$X$与$Y$的关系为：$Y\left(m,X\left(m,\mathrm{1}\right)\right)=\mathrm{1}$。

推广到一般情况，若基本作业区段$m$在天窗$l$时的无轨距不平顺值为${\sigma }_{m,l}^{\mathrm{N}}$，则下一天窗时间$l+\mathrm{1}$的无轨距不平顺值${\sigma }_{m,l+\mathrm{1}}^{\mathrm{N}}$为

$${\sigma }_{m,l+\mathrm{1}}^{\mathrm{N}}={\sigma }_{m,l}^{\mathrm{N}}·{\left(\mathrm{1}+k\right)}^{\mathrm{1}-Y\left(m,l\right)}·{\left[\frac{{\sigma }_{m,l}^{\mathrm{N}}\mathrm{\text{'}}}{{\sigma }_{m,l}^{\mathrm{N}}}\left(\mathrm{1}+k\right)\right]}^{Y\left(m,l\right)}$$

（7）

由式（7）易得，基本作业区段$m$年末无轨距不平顺值${\sigma }_{m,l_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{N}}$计算公式如下：

$${\sigma }_{m,l_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{N}}={\sigma }_{m,l}^{\mathrm{N}}{\prod }_{l=\mathrm{1}}^{l_{max}}{\left(\mathrm{1}+k\right)}^{\mathrm{1}-Y\left(m,l\right)}·{\left[\frac{{\sigma }_{m,l}^{\mathrm{N}}\mathrm{\text{'}}}{{\sigma }_{m,l}^{\mathrm{N}}}\left(\mathrm{1}+k\right)\right]}^{Y\left(m,l\right)}$$

（8）

式中：$l_{\mathrm{m}\mathrm{a}\mathrm{x}}$表示基本作业区段$m$全年天窗数量，为常数。

目标函数应该能反映区段$m$全年的几何不平顺状况，因此应计算全年$l_{\mathrm{m}\mathrm{a}\mathrm{x}}$个天窗时间下$l_{\mathrm{m}\mathrm{a}\mathrm{x}}$个${\sigma }_{m,l}^{\mathrm{N}}$的平均值，即：

$$\overline{{\sigma }_m^{\mathrm{N}}}=\frac{\mathrm{1}}{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\displaystyle \sum _{l=\mathrm{1}}^{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}}{\sigma }_{m,l}^{\mathrm{N}}$$

（9）

除了全年均值外，线路区段不平顺值的标准差也是描述区段全年整体状态的重要参数。在一般情况下，一个区段的不平顺值在全年的分布情况可大致视为正态分布。根据正态分布的$\mathrm{3}\sigma$准则，即正态分布的数据落在$\mu \pm \mathrm{3}\sigma$以外的概率仅0.27 %，此处$\mu$表示为全年无轨距TQI的均值$\overline{{\sigma }_m^{\mathrm{N}}}$。最终可以将目标函数定义为$\mu +\mathrm{3}\sigma$，该函数有以下优点：①有效地反映出作业区段全年无轨距TQI值的波动情况；②有效地规避了由于轨检车的测量误差而带来的列车运行安全问题。为避免字母造成的误解，标准差采用其计算公式来表示。则对于区段$m$来说，其目标函数$F$可表示为

$$F_{\mathrm{m}\mathrm{i}\mathrm{n}}={\displaystyle \sum _{m=\mathrm{1}}^{m_{\mathrm{m}\mathrm{a}\mathrm{x}}}}\left[\overline{{\sigma }_m^{\mathrm{N}}}+\mathrm{3}\times \sqrt[]{\frac{\mathrm{1}}{l_{max}}{\displaystyle \sum _{l=\mathrm{1}}^{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}}({\sigma }_{m,l}^{\mathrm{N}}{-\overline{{\sigma }_m^{\mathrm{N}}})}^{\mathrm{2}}}\right]$$

（10）

1.4　目标函数单调性质

根据式（10）的推导结果，对于基本作业区段$m$，目标函数表达式为

$$F_m=\overline{{\sigma }_m^{\mathrm{N}}}+\mathrm{3}\times \sqrt[]{\frac{\mathrm{1}}{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\displaystyle \sum _{l=\mathrm{1}}^{\mathrm{m}\mathrm{a}\mathrm{x}}}({\sigma }_{m,l}^{\mathrm{N}}{-\overline{{\sigma }_m^{\mathrm{N}}})}^{\mathrm{2}}}$$

（11）

式中：$F_m$表示基本作业区段$m$的目标函数表达式，$\overline{{\sigma }_m^{\mathrm{N}}}$表示基本作业区段$m$的无轨距不平顺值平均值，${\sigma }_{m,l}^{\mathrm{N}}$的递推关系式在式（7）中给出。

可以发现，目标函数$F_m$仅与捣固作业天窗编号$l_{\mathrm{1}}$、打磨作业天窗编号$l_{\mathrm{2}}$、区段不平顺初始值${\sigma }_{m,\mathrm{0}}^{\mathrm{N}}$有关。为方便后续求解最优值，需探求目标函数$F_m$与捣固作业天窗编号$l_{\mathrm{1}}$的关系，故做如下推导：

设基本作业区段$m$恶化函数为

$${\sigma }_m^{\mathrm{N}}\left(t\right)=A·{\mathrm{e}}^{bt}+C$$

（12）

该恶化函数与Quiroga和Schnieder^［15］所建立的几何不平顺预测函数类型一致。式中，自变量$t$单位为d，$A、b、C$为待定系数，${\sigma }_m^{\mathrm{N}}\left(t\right)$单位为mm。为便于之后的求导运算，假设自变量$t$连续，即存在小数天。若该区段在天窗编号为$l_{\mathrm{1}}$时进行线路捣固作业，此时自变量$t$为

$$t={\displaystyle \sum _{i=\mathrm{1}}^{l_{\mathrm{1}}}}T\left(i\right)$$

（13）

式中：函数$T(l)$为天窗时间安排函数，输入天窗编号$l$，输出与下一天窗的间隔时间，d。

设基本作业区段$m$捣固作业恢复函数为

$$J_m\left({\sigma }_m^{\mathrm{N}}\right)=D·\mathrm{l}\mathrm{n}\left({\sigma }_m^{\mathrm{N}}\right)+M$$

（14）

则在捣固作业前，任一天窗时间无轨距不平顺值为

$${\sigma }_{m,l}^{\mathrm{N}}=\left({\sigma }_{m,\mathrm{0}}^{\mathrm{N}}-C\right){\mathrm{e}}^{b\cdot {\displaystyle \sum _{i=\mathrm{1}}^l}T\left(i\right)}+C$$

（15）

同理，在捣固作业后，任一天窗时间无轨距不平顺值为

$${{\sigma }^{\mathrm{N}}}_{m,l}=\{[({{\sigma }^{\mathrm{N}}}_{m,\mathrm{0}}-C){\mathrm{e}}^{b\cdot {\displaystyle \sum _{i=\mathrm{1}}^{l_{\mathrm{1}}}}T\left(i\right)}+C]\cdot [\mathrm{1}-D\mathrm{l}\mathrm{n}(({{\sigma }^{\mathrm{N}}}_{m,\mathrm{0}}-C)e^{b\cdot {\displaystyle \sum _{i=\mathrm{1}}^{l_{\mathrm{1}}}}T\left(i\right)}+C)-F]-C\}\cdot$$

$${\mathrm{e}}^{b\cdot \left[{\displaystyle \sum _{i=\mathrm{1}}^l}T\left(i\right)-{\displaystyle \sum _{i=\mathrm{1}}^{l_{\mathrm{1}}}}T\left(i\right)\right]}+C$$

（16）

通过联立式（14）—（16）可以得到目标函数$F_m$关于捣固作业天窗编号$l_{\mathrm{1}}$的表达式。对函数$F_m(l_{\mathrm{1}})$求二阶导数，计算得到其二阶导数恒大于0，即：

$$\frac{{\mathrm{d}}^{\mathrm{2}}{F}_m}{\mathrm{d}{l_{\mathrm{1}}}^{\mathrm{2}}}>\mathrm{0}$$

（17）

在$l_{\mathrm{1}}\in [\mathrm{1},l_{\mathrm{m}\mathrm{a}\mathrm{x}}]$恒成立。由导数性质可知，一阶导数$F_m\text{'}(l_{\mathrm{1}})$单调递增，故其最多在定义域范围内有唯一零点。因此，函数$F_m(l_{\mathrm{1}})$至多有一个极小值点。这说明，若仅针对基本作业区段$m$，运用枚举法得到最佳作业天窗时间$l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$，则$F_m(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})$必定是定义域内的极小值点或区间端点。当$F_m(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})$为区间端点值时，说明函数在定义域内单调；当$F_m(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}})$为区间极小值点时，则函数在区间$[\mathrm{1},l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}]$及区间$[l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}},l_{\mathrm{m}\mathrm{a}\mathrm{x}}]$两个范围内均单调。综合以上两种情况，可以得出结论：若$l_{\mathrm{1}}$取值远离$l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$，即$|l_{\mathrm{1}}-l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}|$增大时，则目标函数$F_m(l_{\mathrm{1}})$必定会增大。

综上所述，在线路打磨作业天窗时间编号$l_{\mathrm{2}}$固定的情况下，对于一个基本作业区段$m$来说，若根据当前打磨天窗时间编号$l_{\mathrm{2}}$得到了最优捣固作业天窗时间$l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$，则函数随着$|l_{\mathrm{1}}-l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}|$的增大单调递增。至此，得到了目标函数一个重要的单调性质。

1.5　约束函数的建立

在本文制定的计算模型中，决策矩阵$\mathit{X}$取值并没有进行约束。但是现场作业具有复杂性，决策矩阵$\mathit{X}$并不能进行随意取值，故需要对模型的决策变量进行约束。根据现场作业的不同情况，约束函数分为以下几种。

（1）作业区段逻辑约束

在本文制定的计算模型中，时间区段为一个天窗时间。显然，在一个天窗时间内，对于同一个基本作业区段来说，仅能进行一项作业，完成线路捣固作业或线路打磨作业。

定义辅助决策矩阵$\mathit{Z}$：$\mathit{Z}$为$m_{\mathrm{m}\mathrm{a}\mathrm{x}}\times l_{\mathrm{m}\mathrm{a}\mathrm{x}}$的二维矩阵，$Z(i,j)$表示矩阵$\mathit{Z}$第$i$行第$j$列的取值。矩阵的取值规则如下：

$$Z\left(i,j\right)=\left\{\begin{array}{c}\mathrm{1}\mathrm{区域}i\mathrm{在}j\mathrm{天窗}\mathrm{时间}\mathrm{进行}\mathrm{打磨}\mathrm{作业}\\ \mathrm{0}\mathrm{区域}i\mathrm{不在}j\mathrm{天窗}\mathrm{时间}\mathrm{进行}\mathrm{打磨}\mathrm{作业}\end{array}\right.$$

（18）

决策变量$\mathit{X}$与$\mathit{Z}$的关系为：$\mathit{Z}(m,\mathit{X}(m,\mathrm{2}))=\mathrm{1}$。

由式（6）及式（18）得，辅助决策变量$\mathit{Y}$与$\mathit{Z}$都仅由决策变量$\mathit{X}$决定，且为一一对应关系。因此，对于任意时间区段$l$，都应满足：

$${\displaystyle \sum _{m=\mathrm{1}}^{m_{\mathrm{m}\mathrm{a}\mathrm{x}}}}\mathit{Y}\left(m,l\right)+{\displaystyle \sum _{m=\mathrm{1}}^{m_{\mathrm{m}\mathrm{a}\mathrm{x}}}}\mathit{Z}\left(m,l\right)\le \mathrm{1}$$

（19）

这表示对于同一个基本作业区段来说，一个天窗时间内仅能完成一项作业。而且由于大型养路机械数量原因，所有$m_{\mathrm{m}\mathrm{a}\mathrm{x}}$个区段至多仅能有一个区段进行作业。

（2）作业次数上限约束

现场实际作业效果表明：大型养路机械捣固作业对线路的维护效果都较为明显。根据数据的统计分析显示，在周期修为主的体制下大型养路机械线路捣固作业的周期为1~1.5年。同理，对于线路打磨作业来说，钢轨的打磨周期一般为2年左右。因此，在本文的计算模型中，线路捣固作业和打磨作业在一年的时间中，仅能安排一次。

对于任意的基本作业区段$m$，应满足：

$${\displaystyle \sum _{l=\mathrm{1}}^{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}}Y\left(m,l\right)\le \mathrm{1}$$

（20）

且

$${\displaystyle \sum _{l=\mathrm{1}}^{l_{\mathrm{m}\mathrm{a}\mathrm{x}}}}Z\left(m,l\right)\le \mathrm{1}$$

（21）

通过设置该约束可以极大简化模型的计算量，同时也符合实际情况。

（3）恶化状态上限约束

在制定年度养护维修计划时，如果在某个时刻下，线路的TQI值超过了安全限值，就会存在较大安全隐患。因此，对于所有基本作业区段$m$，在任意天窗时间下的无轨距不平顺值不应该超过安全限值${\sigma }_{\mathrm{u}}$，即对于$\forall mϵ[\mathrm{1},m_{\mathrm{m}\mathrm{a}\mathrm{x}}]$，$lϵ[\mathrm{1},l_{\mathrm{m}\mathrm{a}\mathrm{x}}]$，应有：

$${\sigma }_{m,l}<{\sigma }_{\mathrm{u}}$$

（22）

（4）线路捣固作业与打磨作业间隔时间约束

一般来说，若一个作业区段需要进行捣固作业和打磨作业时，应先进行捣固作业，再进行打磨作业。原因如下：①钢轨打磨车很容易将几何不平顺值较大的线路的钢轨打坏；②线路捣固作业会影响打磨作业后的效果。因此，如果先进行线路打磨作业，捣固作业应间隔一定时间。

即对于$\forall mϵ[\mathrm{1},m_{\mathrm{m}\mathrm{a}\mathrm{x}}]$，$l\mathrm{ϵ}[\mathrm{1},l_{max}]$，应有：

$$\mathit{Y}\left(m,l\right)-\mathit{Z}\left(m,l\right)<\mathrm{0}$$

（23）

或当不满足上式时，应满足式（24）：

$$\mathit{Y}\left(m,l\right)-\mathit{Z}\left(m,l\right)>u_{\mathrm{x}}$$

（24）

式中：$u_{\mathrm{x}}$表示线路打磨作业后进行捣固作业应间隔的天窗个数。

（5）线路打磨作业连续性约束

根据假设，该模型的基本作业区段长度为4 km。由于线路打磨作业单次作业长度为8 km，一次线路打磨作业需要对两个连续作业区段进行打磨。

即对于$\forall mϵ[\mathrm{1},m_{\mathrm{m}\mathrm{a}\mathrm{x}}]$，$lϵ[\mathrm{1},l_{\mathrm{m}\mathrm{a}\mathrm{x}}]$，若$Z(m,l)=\mathrm{1}$，则有：

$$\mathit{Z}\left(m-\mathrm{1},l\right)+\mathit{Z}\left(m+\mathrm{1},l\right)=\mathrm{1}$$

（25）

这意味着，单次线路打磨作业必须对两个连续作业区段进行。

2 捣固与打磨作业计划编制模型的求解

有砟轨道养护维修计划编制模型的求解是一个典型的最优化问题。最优化问题就是在一定条件下寻找最佳方案的问题。本节主要针对第1节中的目标函数、约束函数进行求解，主要采用的方法为改进枚举法，在确保运算结果可靠性的前提下提高运算速度。

2.1　问题描述

对于有砟轨道养护维修计划编制模型的求解来说，是一个有约束条件的、非线性的、单目标的、整数规划问题，其数学定义可以用式（26）描述：

$$\mathrm{M}\mathrm{i}\mathrm{n}f\left(x,y\right)$$

（26）

$$\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{g}_i\left(\mathit{x},\mathit{y}\right)\le \mathrm{0},i=\mathrm{1,2},\dots ,p\\ h_r\left(\mathit{x},\mathit{y}\right)=\mathrm{0},r=\mathrm{1,2},\dots ,q\end{array}\right.$$

（27）

式中：$\mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_n)$，$\mathit{y}=(y_{\mathrm{1}},y_{\mathrm{2}},\dots ,y_n)$，二者均为${\mathrm{R}}^n$空间上的$n$维决策向量；$f(x,y)$为目标函数；$g_i(x,y)$为不等式约束条件，$h_r(x,y)$为等式约束条件。

2.2　求解算法

枚举法也称穷举法，在计算时可以考察某最优化问题的所有可能情况，最终得到最优解。枚举法利用了计算机运算速度快、精确度高的特点，对要解决问题的所有可能情况，一个不漏地进行检验，从中找出符合要求的答案，因此枚举法是解决单目标最优化问题的一个很好方法。

传统枚举法有以下几个特点：①通过枚举法计算得到的结果一定是正确的；②计算效率可能做了很多无用功，效率较为低下；③如果计算量大，可能会造成时间崩溃。为了能将枚举法有效地运用在本文所建立的模型求解中，同时提高运算速度，需要对该方法进行适当改进。改进思路如图2所示。

图2  改进枚举法求解路线图

Fig. 2  Solution roadmap of improved enumeration method

整体计算思路如下：根据目标函数进行求解，首先仅考虑1.5节中约束条件（2）、（3）、（4），利用枚举法进行全局搜索求解，得到初始解。再进行search函数与permu函数反复递归，使其满足1.5节中约束条件（1）和（5）得到最终解。下面对这两个简化函数进行详细介绍。

（1）search函数

本模型需对多个基本作业区段的目标函数最优值进行求解。设基本作业区段共有$M$个，各个基本作业区段均独立，则对于每个基本作业区段来说，都可以得到一组最优解$(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^m,l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^m)$（上标表示基本作业区段编号）。此时，可以得到一组初始解${\mathit{X}}_{\mathrm{0}}$，初始解${\mathit{X}}_{\mathrm{0}}$取值为

$${\mathit{X}}_{\mathrm{0}}=\left[\begin{array}{cc}{l}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{1}}\left(\mathrm{0}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{1}}\left(\mathrm{0}\right)\\ l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{2}}\left(\mathrm{0}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{2}}\left(\mathrm{0}\right)\\ ⋮& ⋮\\ l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{M}}\left(\mathrm{0}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{M}}\left(\mathrm{0}\right)\end{array}\right]$$

（28）

式中：括号中的数字表示枚举法计算次数，0表示未计算。若初始解${\mathit{X}}_{\mathrm{0}}$符合作业区段逻辑约束（约束条件1）和线路打磨作业连续性约束（约束条件（5）），即可直接输出结果。这里考虑一般情况，在绝大部分情况下，初始解$X_{\mathrm{0}}$不符合约束条件1和约束条件（5）。

考虑到线路打磨作业需要对两个连续基本作业区段进行作业，则以相邻两个作业区段为一组，打磨作业天窗区段$l_{\mathrm{2}}$进行最优值搜寻。不妨设任意两个相邻作业区段的初始决策变量$X_{\mathrm{0}}^m=(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^m,l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^m)$，$X_{\mathrm{0}}^{m+\mathrm{1}}=(l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{m+\mathrm{1}},l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{m+\mathrm{1}})$，则对于区段$m$与区段$m+\mathrm{1}$来说，其打磨作业时间区段$l_{\mathrm{2}}$的枚举法搜寻上下限为$[l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^m,l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{m+\mathrm{1}}]$。在确定了打磨作业天窗时间后，可以利用枚举法针对捣固作业进行比选，根据1.4对目标函数性质的分析，在打磨作业时间区段$l_{\mathrm{2}}$确定的情况下，目标函数是捣固作业时间区段$l_{\mathrm{1}}$的一元函数，有唯一极值点或最值点。故可以得到第一次枚举法计算结果临时解$X_{\mathrm{1}}$，取值为

$${\mathit{X}}_{\mathrm{1}}=\left[\begin{array}{cc}{l}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{1}}\left(\mathrm{1}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{1}}\left(\mathrm{1}\right)\\ l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{2}}\left(\mathrm{1}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{2}}\left(\mathrm{1}\right)\\ ⋮& ⋮\\ l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{1}}^{\mathrm{M}}\left(\mathrm{1}\right)& l_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{2}}^{\mathrm{M}}\left(\mathrm{1}\right)\end{array}\right]$$

（29）

在令枚举计算操作为search，则有${\mathit{X}}_{\mathrm{1}}=\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h}({\mathit{X}}_{\mathrm{0}})$。

（2） permu函数

${\mathit{X}}_{\mathrm{1}}$与${\mathit{X}}_{\mathrm{0}}$的区别在于，临时解${\mathit{X}}_{\mathrm{1}}$符合约束条件（5），但不一定符合约束条件（1），而初始解${\mathit{X}}_{\mathrm{0}}$约束条件（1）、（5）均不一定符合。对于临时解${\mathit{X}}_{\mathrm{1}}$来说，捣固作业和打磨作业共占据$M+[\frac{M}{\mathrm{2}}]$个天窗时间，这些天窗时间不应重合，否则将违反约束条件（1）。若临时解${\mathit{X}}_{\mathrm{1}}$不满足约束条件1，不妨设重合的天窗时间区段集合为$L=\{l_{\mathrm{1}}^{\mathrm{c}},l_{\mathrm{2}}^{\mathrm{c}},\dots ,l_n^{\mathrm{c}}\}$，对$L$内元素进行升序重排列，天窗时间变化幅度为重复元素个数。例如，某三个作业区段均在天窗编号56进行作业，则重排列后的解集集合为$\{\mathrm{54,55,56}\}$或$\{\mathrm{55,56,57}\}$或$\{\mathrm{56,57,58}\}$。令$K$表示重排列后解集集合的集合，即$K=\{k_{\mathrm{1}},k_{\mathrm{2}},\dots ,k_n\}$。将临时解${\mathit{X}}_{\mathrm{1}}$进行重排列微扰动操作后，取其中一个集合元素排列$k_i$，得到临时解${\mathit{X}}_{\mathrm{1}}\text{'}$。

若将这次重排列微扰动操作视为permu函数（permutation，重排列）操作，则有${\mathit{X}}_{\mathrm{1}}\text{'}=\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{u}({\mathit{X}}_{\mathrm{1}})$。由目标函数单调性和零点定理可知，最优解必定在集合$K$中。同样的，对${\mathit{X}}_{\mathrm{1}}\text{'}$进行search函数操作，得到临时解${\mathit{X}}_{\mathrm{2}}$，即$X_{\mathrm{2}}=\mathrm{s}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h}({\mathit{X}}_{\mathrm{1}}\text{'})$。对得到的临时解${\mathit{X}}_{\mathrm{2}}$，若符合约束条件1，则输出${\mathit{X}}_{\mathrm{2}}$为最终解$\mathit{X}$；若不符合约束条件（1），则再次进行上述操作，直至符合约束条件（1）。

从流程图可以看出，最终解符合所有约束函数，并根据目标函数的性质可以推知，最终解为最优解。search函数操作目的是让解符合约束条件（2）、（3）、（4），而permu函数操作目的是让解符合约束条件（1）、（5），在缩小可行解范围的前提下，进行枚举法搜索，并结合目标函数性质得到符合所有约束函数的最优解$\mathit{X}$，而最优解$\mathit{X}$即为所求。至此完成了利用改进枚举法求解目标函数。

3 实例分析

为了验证改进的枚举法的有效性以及建立最优化模型的实用性，本节选取了部分线路进行有砟轨道年度养护维修计划的编制。

3.1　模型参数

本节选取共6个基本作业区段作为研究对象，共24.33 km。总里程不是25.4 km，这是因为该区段有道岔区段，应进行去除。6个区段的无轨距的不平顺初始值及打磨作业情况见表1。

表1  2019年京沪线上行K1364+000~K1389+400里程区段信息

Tab. 1  Information of K1364 + 000 to K1389 + 400 mileage section of Beijing-Shanghai Railway in 2019

作业区段序号	起点里程	终点里程	初始无轨距不平顺值 /mm	距上次打磨作业过去的时间/d
1	K1364+000	K1368+000	2.65	455
2	K1368+000	K1372+000	2.08	455
3	K1372+320	K1376+320	1.92	456
4	K1376+810	K1380+810	2.91	456
5	K1380+810	K1385+140	3.41	557
6	K1385+400	K1389+400	3.28	557

3.2　模型求解

将现场数据及计算参数输入至模型中，计算得到最终$\mathit{X}$：

$$\mathit{X}=\left[\begin{array}{cc}\mathrm{6}& \mathrm{3}\\ \mathrm{7}& \mathrm{3}\\ \mathrm{40}& \mathrm{11}\\ \mathrm{5}& \mathrm{11}\\ \mathrm{8}& \mathrm{4}\\ \mathrm{2}& \mathrm{4}\end{array}\right]$$

（30）

模型总求解时间为186 s，计算效率得到了提升。根据最终解X得到这6个区段的年度捣固与打磨作业计划安排，计算出根据该安排下的目标函数，具体值见表2。

表2  2019年京沪线上行K1364+000~K1389+400区段模型计算结果

Tab. 2  Model calculation results of K1364 + 000 to K1389 + 400 mileage section of Beijing- Shanghai Railway in 2019

区段编号$m$	捣固作业天窗编号$l_1$	打磨作业天窗编号$l_2$	目标函数/mm
1	6	3	2.24
2	7	3	2.17
3	40	11	2.79
4	5	11	2.68
5	8	4	2.67
6	2	4	2.58
平均值			2.52

3.3　结果分析

根据2019年现场大机作业安排表，整理得到现场作业天窗时间数据，根据现场作业数据得到的目标值计算结果见表3。

表3  2019年京沪铁路上行K1364+000~K1389+400区段年度作业计划安排表

Tab. 3  Annual operation schedule of K1364 + 000 to K1389 + 400 mileage section of Beijing-Shanghai Railway in 2019

区段编号$m$	捣固作业天窗编号$l_1$	打磨作业天窗编号$l_2$	目标函数（mm）
1	79	37	3.71
2	78	37	3.41
3	77	36	2.96
4	0	36	4.23
5	0	35	4.70
6	74	35	3.54
平均值			3.76

表中天窗编号为0表示不进行作业。根据现场作业计划安排，可以得到各区段全年无轨距不平顺值变化如图3所示。

图3  各区段全年无轨距不平顺值变化折线图（现场作业条件下）

Fig. 3  Broken line diagram of annual non gauge irregularity value change of each section(under site operation conditions)

可以看到，区段4和区段5在2019年并未安排捣固作业，而本文建立的模型求解结果是有捣固作业计划的。为方便对比，这里将模型计算结果中区段4和区段5决策变量捣固作业部分强制设0，得到新的模型计算结果如表4。

表4  2019年京沪铁路上行K1364+000~K1389+400区段模型计算结果（新）

Tab. 4  Model calculation results of K1364 + 000 to K1389 + 400 mileage section of Beijing- Shanghai Railway in 2019 (new)

区段编号$m$	捣固作业天窗编号$l_1$	打磨作业天窗编号$l_2$	目标函数/mm
1	6	3	2.24
2	7	3	2.17
3	40	11	2.79
4	0	11	4.23
5	0	4	4.70
6	2	4	2.58
平均值			3.12

同样的，按照新模型求解得到的作业计划安排，可以得到各区段全年无轨距不平顺值的变化图如图4所示。

图4  各区段全年无轨距不平顺值变化折线图（模型条件下）

Fig. 4  Broken line diagram of annual non gauge irregularity value change of each section(under model conditions)

从图3、图4可以看出，根据所建立的模型求解得到的作业计划安排下，各区段全年无轨距轨道不平顺值均可保持到一个较低水平。从表中可以看出，按照2019年上海工务段的大机作业安排来确定模型中的决策变量取值，得到目标函数平均值为3.76 mm；而根据最优化模型和改进的枚举法计算得到的结果，在进行同样捣固次数前提下，目标函数平均值为3.12 mm。以目标函数为指标来看，可以认为，模型计算结果比现场作业结果优化了约11.8 %；若不考虑捣固次数限制，则模型计算结果比现场作业结果优化了33 %。以无轨距轨道几何状态不平顺值的最大值为指标，在实际现场作业计划下，区段6无轨距不平顺值的最大值达到了5.49 mm，而在最优化模型求解结果的作业计划下，无轨距不平顺值的最大值仅为4.39 mm，改善了约19.9 %。

从上述结果可以看出，现场作业并未结合各作业区段轨道几何不平顺的发展状态，因此作业效果并不好。而本文所建立的大机作业决策模型求得的结果综合考虑了各个区段的起始不平顺、恶化速率、作业效果等因素，使得捣固作业和打磨作业效果更好，最终得到年度轨道几何状态均值更低，证明了模型的实用性和求解算法的有效性。

4 结论

本文根据已有文献中建立的轨道几何不平顺预测模型，完成了年度轨道几何状态最优目标函数、现场作业约束函数的推导，分析了目标函数的单调性质，建立了有砟轨道捣固与打磨养护维修计划编制数学模型。然后对单目标最优化问题的求解方法进行了讨论，利用改进的枚举法完成了对模型的求解。最后，选取京沪线上行K1364+000~K1389+400里程连续区段进行实例分析，通过对比现场作业计划和模型求解结果，对模型和算法的可靠性与适用性进行了验证，主要内容如下：

（1）建立了有砟轨道捣固与打磨作业养护维修计划编制数学模型。该模型可以根据轨道几何状态的发展进行合理的大机作业安排。同时，可以确定捣固作业时机，填补了打磨作业时机预测的空白；

（2）根据有砟轨道养护维修计划编制模型的特点，阐述了改进枚举法的思路，完成了枚举法的改进。该改进枚举法不仅能够准确计算结果，也提高了运算的效率；

（3）选取京沪线上行K1364+000~K1389+400里程连续区段进行模型验证。通过对比现场作业计划安排和模型求解结果发现：以目标函数为指标，优化了约11.8 %；以无轨距轨道几何状态不平顺值的最大值为指标，优化了约19.9 %。

通过对计算结果的分析，可以认为本文建立的有砟轨道养护维修计划编制模型及改进的枚举法计算结果可靠、真实。根据得到的模型计算结果可以为现场捣固作业和打磨作业计划编制提供参考，对现有养修策略模型进行完善。

作者贡献声明

许玉德：方案制定；过程监督；论文校修。

张 喜：模型确定；数据处理；论文写作。

刘思磊：模型确定；论文校修。

吴宣庆：模型确定；数据分析。