基于子集模拟的建设项目离散型工期 - 成本优化算法研究

王 家， 韩淙吉， 陈雅含， 李炎彪

(湖南大学 土木工程学院， 湖南 长沙 410082)

工期和成本是建设项目管理的重要目标，且两者相互联系、相互制约。项目工期的压缩，往往伴随着劳动力的增加、更高效率施工设备和施工方法的采用等，进而导致项目成本的增加[1]。因此，在进度计划管理决策时，需综合考虑工期与成本间的合理权衡，考察项目的工期 - 成本优化问题。同时，建设项目的资源(如施工设备、劳动力)分配通常具有离散型的特征，项目中各工序的实施方式只能从有限个执行模式中选取[2～5]。因此，离散型工期 - 成本优化问题更为贴合工程实际，受到项目管理者和研究人员的广泛关注。

离散型工期 - 成本优化问题属于组合优化问题，其求解方法主要分为精确算法和启发式算法。其中，精确算法以整数规划为代表，旨在寻求问题的准确最优解[6]。但离散型工期 - 成本优化问题为NP困难问题，其复杂程度随项目中工序数量的增加而急剧上升[7]。当优化问题中项目的工序数量较多时，精确算法并不适用。此时，适宜采用启发式算法，如遗传算法[8～11]、粒子群算法[12]、蚁群算法[13,14]等，来寻找问题的最优解或近似最优解。但是，启发式算法的求解具有随机性，每次运行获得的最优解不一定相同(不稳定)，现有的启发式算法在最优解获取稳定性上仍有较大的改进空间。

为此，本文针对建设项目的离散型工期 - 成本优化问题，提出一种基于子集模拟法的新型启发式算法。同时，为解决算法中可行域及收缩区域中随机样本点的高效产生难题，本文建议采用“带有反射壁的随机游动”的马尔科夫链蒙特卡罗方法。通过算例验证，与应用较广的遗传算法相比，本文建议的优化算法性能更好，在最优解的获取稳定性上有较大改进。

1 建设项目离散型工期 - 成本优化问题模型

建设项目离散型工期 - 成本优化问题中，涉及的项目工序执行模式为有限个。因此，考虑一包含N项工序的建设项目，并设任一工序i的执行模式总数为Ki，i=1,…,N。项目工序的不同执行模式，可能对应于不同的施工工艺，也可能对应于不同的的人力、机械投入数量，进而导致不同的工序成本和工序工期。例如，表1为某桥梁工程项目中大孔板梁安装工序的3种备选执行模式，对应于3种不同的施工工艺，其对应的成本和工期亦各不相同。

表1 大孔板梁安装工序的备选执行模式

假设工序i在第j种执行模式下的持续时间和直接成本分别为ti,j和ci,j。优化模型中，一般采用xi,j作为执行模式的指示变量，即当工序i选择执行模式j时，xi,j=1，否则xi,j=0。因任一工序同时只能执行一种模式，所以指示变量xi,j需满足：

(1)

项目的实施方案，由项目中所有工序的执行模式所确定，可表示为包含所有指示变量的矩阵X=[xij,i=1,…,N,j=1,…,Ki]。例如，表2为某建设项目的实施方案，该项目包含6项工序，每项工序的执行模式总数分别为5,4,3,4,4,5。由指示变量xi,j的含义可知，该项目实施方案中，工序1～6分别选择第3、第2、第2、第4、第1及第4种执行模式。

表2 某建设项目实施方案矩阵X

对应项目的某一实施方案X，工序i的持续时间ti和直接成本ci可表示为：

(2)

(3)

项目的总工期T(X)，可采用关键路径法，由关键路径上所有工序的持续时间求和得到[15]：

(4)

式中:A为关键路径上的工序集合。考虑项目的要求完工工期TD时，项目的实施方案X需满足约束条件：

T(X)≤TD

(5)

项目的总成本C(X)包含直接总成本Cd(X)和间接总成本Ci(X)。其中，直接总成本Cd(X)对应于所有工序的直接成本之和。间接总成本Ci(X)一般假定为项目总工期的线性函数，通过引入单位工期的间接成本l来计算[16]。因此，对应于某一项目实施方案，项目的总成本C(X)可表示为：

C(X)=Cd(X)+Ci(X)

=Cd(X)+l×T(X)

(6)

在此基础上，建设项目的离散型工期 - 成本优化模型可描述为如下的0-1型整数规划问题[17]：

(7)

建设项目的实施方案，也可采用向量的表示形式，写为Z=[Z1,Z2,…,ZN]。此时，向量Z中的元素Zi，i=1,…,N代表工序i的执行模式编号，为[1,Ki]内的自然数。由项目实施方案的矩阵X与向量Z的含义可知，X与Z间存在一一对应关系。例如，表2的矩阵X对应的实施方案中，工序1～6分别选择第3、第2、第2、第4、第1及第4种执行模式，因此该方案的向量表示为Z=[3,2,2,4,1,4]。由向量Z的含义可知，其对应的矩阵X自动满足公式(1)对应约束条件。因此，在项目实施方案的向量表示下，建设项目的离散型工期 - 成本优化模型可写为：

(8)

其中，为得到实施方案Z下的项目总工期T(Z)和总成本C(Z)，可首先求得Z对应的实施方案矩阵X(Z)，进而利用公式(4)和公式(6)计算。

2 基于子集模拟的离散型工期 - 成本优化算法

子集模拟法(Subset Simulation)是由Au和Beck提出的、针对高维空间中可靠度估算分析的高效数值模拟算法[18]。该方法通过引入一组中间失效事件，将偶发事件的小概率表示为一组较大的条件概率的乘积，进而将小概率估算问题转化为较大条件概率的估计问题，以提高数值模拟算法的计算效率。通过建立可靠度问题与优化问题间的内在联系，子集模拟法在结构优化问题中取得了较好应用[19]。但是，在建设项目工期 - 成本优化问题中，尚未见到基于子集模拟法的相关应用。

2.1 基于子集模拟的优化框架

针对公式(8)的工期 - 成本优化问题，其可行域为Ω0={Z:T(Z)≤TD,1≤Zi≤Ki,i=1,…,N}。基于子集模拟的优化算法中，通过引入一系列与目标函数相关的、递减的边界值c1>c2>…>cJ，在可行域Ω0内定义了一组逐渐缩小的区域Ωj={Z:Z∈Ω0,C(Z)≤cj}，j=1,…,J。其中，定义的区域间存在子集关系，即Ω0⊇Ω1⊇Ω2⊇…⊇ΩJ。基于子集模拟的优化算法中，通过逐次在可行域Ω0及逐渐缩窄的区域Ωj(j=1,…,J)中按均匀分布随机取样，将搜索范围逐渐缩窄到最优解附近的较小区域，以得到问题的最优解。

基于子集模拟的优化算法中有两个重要参数：指定区域中的抽样样本数M和条件概率参数p0。参数M和p0的取值，需保证Mp0和1/p0均为正整数。条件概率参数p0决定了迭代中考察区域的收缩快慢，如p0过大，迭代中考察区域的收缩变慢，优化算法需要较多的迭代次数。同时，Mp0对应于每一个迭代中收缩区域的初始样本点数量，如p0过小，则需要较大的样本数量M才能保证Mp0的数值不至过小，以达到对收缩区域的有效搜索。综合目前研究成果，p0一般在0.05～0.3区间取值，进而可根据计算资源的承受能力确定抽样样本数量M。

图1～4描述了参数M=25和p0=0.2时的上述流程。针对图1中随机产生(满足可行域Ω0内的均匀分布)的25个抽样样本点，图2将它们的目标函数值按升序排列。此时，定义下一收缩区域Ω1的临界值c1，可由这25个目标函数值中第Mp0=5小的数值确定，见图2。对应上述确定的c1，可行域Ω0内随机产生的M=25个抽样样本点，有Mp0=5个落在收缩区域Ω1={Z:Z∈Ω0,C(Z)≤c1}内，见图3。为产生收缩区域Ω1内的M=25个抽样样本点，以目前的Mp0=5个样本点为种子，以Ω1内的均匀分布为目标概率分布，采用MCMCS方法产生5条马尔科夫链链条，见图4。每条链条以其中一个种子为起点，并额外产生(1/p0-1)=4个状态点，以确保总样本数为Mp0(1/p0-1+1)=M=25。

图1 可行域Ω0内随机产生的抽样样本点

图2 定义收缩区域Ω1的边界值c1

图3 收缩区域Ω1及落入其中的初始样本点

图4 收缩区域Ω1的抽样样本点

2.2 收缩区域Ωj={Z:Z∈Ω0,C(Z)≤cj}内随机样本点的产生

如2.1节所述，基于子集模拟的优化算法中，需根据收缩区域Ωj内的少量样本点，利用MCMCS方法产生要求数量的更多样本点。本文采用“带有反射壁的随机游动”方法[21]来实现马尔科夫链的产生。具体而言，为产生收缩区域Ωj内Z1为起点的马尔科夫链{Z1,Z2,…}，采用如下Zk=[Zk,1,Zk,2,…,Zk,N]到Zk+1=[Zk+1,1,Zk+1,2,…,Zk+1,N]的迭代操作：

(1)从样本点Zk=[Zk,1,Zk,2,…,Zk,N]中，随机选择一可变元素Zk,r(对应工序r的执行模式总数Kr>1)。

2.3 可行域Ω0内随机样本点的产生

2.4 基于子集模拟的离散型工期 - 成本优化算法总结

针对建设项目的离散型工期 - 成本优化问题，本文基于子集模拟法，建议的优化算法的步骤总结如下：

(1)确定优化算法中采用的参数，即在指定区域中随机抽样的样本数量M、条件概率参数p0。

(4)算法采用的迭代终止原则为：消耗的计算资源达到设定的上限值。如不满足迭代终止条件，令j=j+1，并重复步骤3。

3 案例分析

为检验基于子集模拟的优化算法的性能，本文考虑文献[22]中的两个离散型工期 - 成本优化问题。两个问题分别包含18项工序和63项工序，为项目离散型工期-成本优化问题的经典算例。同时，为增加算例的复杂程度，本文将18项工序的算例网络首尾相接、重复3次，以构成包含54项工序的新算例。考虑到版面限制，项目中各工序在不同执行模式下的工期和成本信息在此不再列出，请参考文献[22]中的表1及表3。

针对算例问题的算法实现，笔者采用MatLab2020a自主编程完成，使用的计算机硬件配置为处理器Intel Core i7-7700、内存16 G。

3.1 算例一(包含54项工序，真实最优解已知)

图5 算例一项目的单代号网络图

图6 建议优化算法中最优目标函数值随迭代阶段的变化

为检验基于子集模拟的建议优化算法的求解稳定性，表3给出了建议优化算法100次独立运行求解后的统计结果。由第2节优化算法的阐述可知，算法中的计算资源主要用于给定实施方案Z下的项目总工期T(Z)和总成本C(Z)的分析。因此，迭代终止条件中的计算资源限制取60000次的项目工期和成本分析。目前，针对工期 - 成本优化问题的启发式算法间的对比研究较少，学界对各种启发式算法的优劣未达成共识。同时，遗传算法因其自行概率搜索、运算并行性、应用不依赖问题种类的强鲁棒性等特点，在工期 - 成本优化问题中的应用更为广泛[8～11]。因此， 本文选择遗传算法与建议算法进行对比分析。考虑到种群规模m、交叉概率参数pc和变异概率参数pm对遗传算法性能的影响，本文依据这些参数的一般取值范围，进行了大量不同参数组合下的遗传算法性能检验。检验发现，针对本算例的离散型工期 - 成本优化问题，遗传算法在参数取值m=100，pc=0.1和pm=0.01时的求解性能最优。表3提供了上述参数取值下遗传算法100次独立运行求解后的统计结果。其中，遗传算法求解采用相同的计算资源限制，即60000次的项目工期和成本分析。

表3 基于子集模拟的建议优化算法与遗传算法独立运行100次后的统计结果

由前人的研究成果可知，该问题的真实最优解对应目标函数值为813810[22]。由表3可知，基于子集模拟的建议优化算法100次运行求解中，均找到了问题的真实最优解(对应目标函数值813810)。相比之下，遗传算法的100次运行求解中，仅有78次获得了问题的真实最优解。由此可见，基于子集模拟的建议优化算法获取最优解的稳定性更高。

3.2 算例二(包含63项工序，真实最优解未知)

图7 算例二项目的单代号网络图

表4提供了上述参数取值下建议优化算法与遗传算法100次独立运行求解后的统计结果。对比可知，建议优化算法获得最优目标函数值的平均值为5468054，小于遗传算法的平均值5487892。同时，建议优化算法获得最优目标函数值的最大值(5493430)、最小值(5453770)均优于遗传算法获得最优目标函数值的最大值(5537150)、最小值(5457580)，且最优目标函数值的标准差更小。图8给出了两种算法100次独立运行获得的最优目标函数值的分布情况。由图8可见，建议优化算法获得的最优目标函数值更小，且分布更为集中，有91%(高于遗传算法的47%)的最优目标函数值分布在区间[5450000，5483000)中，表明基于子集模拟的建议优化算法获取最优解的稳定性更高。

表4 基于子集模拟的建议优化算法与遗传算法独立运行100次获得的最优目标函数值的统计结果

图8 两种优化算法独立运行100次获得的最优目标函数值的分布

4 结 论

针对进度计划管理决策中常见的离散型工期 - 成本优化问题，本文提出了基于子集模拟法的优化算法。同时，本文采用“带有反射壁的随机游动”的马尔科夫链蒙特卡罗方法，重点解决了优化算法中可行域及收缩区域中随机样本点的产生问题。此外，本文利用包含54项工序和63项工序的两个算例，检验了建议优化算法的性能，并与应用较广的遗传算法进行了对比。在包含54项工序，真实最优解已知的算例中，建议优化算法获取真实最优解的频率为100%，高于遗传算法的78%。而在包含63项工序，真实最优解未知的算例中，建议优化算法获取的最优目标函数值更优，且更为集中。研究表明，相较于遗传算法，基于子集模拟的建议优化算法在最优解的获取稳定性上有较大改进。