基金项目:湖北能源集团科技攻关项目“湖北能源气化长江工程战略布局规划研究”,EN0T-ZX-F2018-100。
● Received: 2024-05-09● Revised: 2024-07-16● Online: 2024-08-15
1.武汉软件工程职业学院(武汉开放大学)商学院;2.湖北能源集团股份有限公司
1.Business School, Wuhan Vocational College of Software and Engineering (Wuhan Open University); 2.Hubei Energy Group Co., Ltd.
natural gas demand, random forest (RF), Sobol, BorutaShap, honey badger algorithm (HBA), Levy flight, XGBoost
DOI: 10.6047/j.issn.1000-8241.2025.01.011
英国能源研究所发布的《世界能源年统计年鉴(2023)》显示,中国天然气消费量仍呈现出高需求态势,这对天然气需求量合理预测与分配提出了更高要求[1]。天然气需求受诸多因素影响,其时序数据主要呈现出非线性、波动性等特征。目前,天然气需求预测主要采用时间序列法、灰色模型、多元线性回归等传统预测方法与机器学习算法。
在传统天然气需求预测领域,李洪兵等[2]从天然气供给、需求、市场、社会4个维度筛选出10个预警指标并构建了天然气供需预测预警指标体系,采用GM(1,1)模型对中国天然气供需发展趋势进行了预测。贺志明等[3]采用多因素回归模型筛选天然气需求关键影响因素,提出了多重影响因素天然气需求预测理论模型。林文辉等[4]通过研究构建 GM(1, 1)-ARIMA(Autoregressive Integrated Moving Average)-LR(Logistic Regression)天然气集成预测模型,基于国家统计局1973—2022年的历史数据集对天然气产量进行预测研究。Li等[5]基于灰色关联度法筛选关键影响因素,构建了有效驱动因素的逐步回归双对数需求函数模型,并对中国2020—2035年天然气需求量进行了预测。Chen等[6]构建具有外生变量的函数自回归模型,对德国天然气需求进行短期预测。Shahbaz等[7−9]采用多元框架模型,分别研究了经济增长、劳动力、资本等要素与天然气需求的相互关系。近年来,机器学习模型在天然气需求研究领域得到广泛关注。Aras[10]采用遗传算法对土耳其住宅天然气消费量月度数据进行预测研究。佟敏等[11]利用灰色关联度法、平均影响值法及主成分分析法进行特征筛选与降维,采用BP(Back Propagation)模型分别结合3种特征降维方式对2019—2025年天然气需求量进行预测。温泉等[12]采用灰色关联度与皮尔逊相关系数法进行特征筛选,构建了改进 Tent-APSO(Adaptive Particle Swarm Optimization)-LSTM(Long Short Term Memory)预测模型,并对2021—2030年中国天然气需求量进行预测。田文才等[13]利用小波变换对天然气负荷数据进行分解,引入Sobol序列改进麻雀搜索算法,对LSTM预测模型进行超参寻优,有效提升了模型的预测精度。王可睿等[14]采用相关系数法筛选影响天然气负荷的气象特征因素,构建了 LSTM-TCN(Temporal Convolutional Network)-Attention天然气负荷组合预测模型,对陕西省西安市某天然气站的居民用气量进行了预测研究。颜珂等[15]构建了 CEEMD(Complementary Ensemble Empirical Mode Decomposition)-Bayes-LSTM天然气需求预测模型,采用2019—2021年共1 066天的西部某城市天然气负荷数据对模型进行实证研究。
目前,天然气需求预测相关研究主要集中于天然气日需求、年需求方向,鲜有基于月度需求预测的研究成果,高维特征序列数据缺失、时间粒度不一致等现实问题也鲜有考量。此外,传统预测方法往往关注单一标度区间的变化趋势,忽略了受多种因素叠加影响的数据序列在不同标度变化趋势的研究。而多重分形消除趋势波动分析法(Multi-Fractal Detrended Fluctuation Analysis, MF-DFA)具有多标度分析能力,不仅可精确揭示数据中的长程相关性,还能有效消除数据中的非平稳趋势,目前已广泛应用于GDP、交通流、径流量等多领域时序数据的长程相关性与特征分形研究,但尚未应用于天然气需求量预测领域。在此,首先引入MF-DFA 法对天然气需求时序数据进行分析,采用二次插值法及随机森林(Random Forest, RF)插值法处理特征序列数据缺失及时间粒度不一致问题。对比选择BorutaShap算法进行特征序列筛选,采用 Sobol低差异序列、改进密度因子、莱维飞行策略提升蜜獾优化算法(Honey Badger Algorithm, HBA)的寻优效能,对极限梯度提升(eXtreme Gradient Boosting, XGBoost)模型中决策树数量、决策树深度、学习速率等决定模型拟合能力的超参组合进行寻优,并将迭代完成后的最优超参组合赋予XGBoost模型进行预测研究。
去趋势波动分析(Detrended Fluctuations Analysis,DFA)是由Peng等[16]在研究核糖核酸行为机理过程中提出的一种分析幂律函数关系的方法,广泛应用于多领域非平稳时间序列长程相关性的研究工作之中。Kantelhardt等[17]在此基础上提出了MF-DFA,该方法不仅可检测长程相关性、确定其标度不变性,还能判断序列是否具有多重分形属性并确定多重分形特征。
对于研究对象影响因素众多的情况,特征序列筛选是降低模型输入维度、最大限度保留特征数据关键信息以及剔除非关键影响因素的重要环节。Boruta算法[18]是基于RF分类算法的包装器特征选择方法, SHAP(Shapley Addictive exPlanation)则是基于Shapely值理论,对比不同条件下边际贡献均值的特征筛选方法。BorutaShap算法[19]融合了Boruta算法与Shapley值的优势,在处理高维与非线性数据领域表现优异,在计算速度与所产生特征子集的质量方面均优于原单一排列重要性筛选方法。因此,拟采用BorutaShap算法进行特征序列筛选。
Chen等[20]以梯度提升决策树(Gradient Boosting Decision Tree, GBDT)为基础提出了 XGBoost模型,其属于集成学习算法,在并行计算效率、缺失值处理、预测性能方面优势明显,已广泛应用于各类机器学习竞赛与工程中。XGBoost的核心是通过迭代建立一系列弱学习器(决策树)并组合为一个强学习器,实现对数据的拟合与预测,其预测值函数与目标函数为:
式中:O为目标函数; ˆyη 为第η个样本的预测值;k为第k个回归树;K为回归树数量;fk为回归树函数;F为回归树的函数空间;zη为第η个样本点的特征向量;l(yη,ˆyη) 为损失函数; γ 为叶子节点数的惩罚因子;V为叶子节点数目;λ为正则化项的惩罚因子;ωj为第j个叶子节点的权重;N为样本数量。
HBA由Hashim等[21]提出,主要模拟蜜獾挖掘与寻找蜂蜜的动态搜索行为。其结构简单,收敛速度较快,在 CEC17(Congress on Evolutionary Computation)测试函数集选择多维数测试条件下,测试结果在绝大多数测试项中均优于常用优化算法,因此选用HBA优化XGBoost模型。
1)初始化阶段。种群粒子中第i个蜜獾的位置xi表达式如下:
式中:ui、li分别为种群粒子i的搜索空间上限与下限;r1为取值区间[0,1]的随机数。
2)定义密度因子。气味密度与猎物的集中强度有关,也受猎物与第i个蜜獾间的距离影响,即:
式中:Ii为猎物对于第i个蜜獾的气味密度;S为猎物集中强度;di为猎物与第i个蜜獾间的距离;xp为猎物迄今为止的最佳位置;r2为取值区间[0,1]的随机数。
3)更新密度因子。密度因子 ∂ 的表达式如下:
式中:C为常量,C≥1(默认设置C=2);T为当前迭代次数;Tmax为最大迭代次数。
4)跳出局部最优。为避免HBA过早陷入局部最优的搜索状态,本阶段引入一个基于[0,1]分布随机数的随机变量G,从而改变搜索方向并尽量扩大搜索范围。
5)更新代理位置。在挖掘阶段蜜獾新位置xn的表达式为:
式中:β为蜜獾获得食物的能力(β≥1,默认设置β=6);r3、r4、r5、r6分别为取值区间[0,1]的不同随机数。
若r6>0.5,则将 G为−1代入蜜蜂吸引阶段表达式:
式中:r7为取值区间[0,1]的随机数。
低差异序列在给定空间内较伪随机序列分布更加均匀,在任意长度的子序列都能均匀填充整个函数空间。Sobol序列作为一种低差异序列,在种群初始化应用中分布均匀、覆盖率佳、计算周期短、采样速度快,在处理高维数据中具有更好的效果[13]。因此,拟采用Sobol序列对种群粒子进行初始化。设种群规模为100,维度取2,种群粒子取值范围均为[0,1],Sobol序列产生的随机数Φ∈[0,1],种群粒子初始化公式为:
Sobol低差异序列较随机序列产生的初始种群粒子分布更均匀,差异性较小且分布质量更高(图1、图2)。因此,拟采用Sobol序列对种群粒子进行初始化。
图1 Sobol序列初始化种群图Fig. 1 Population initialization of Sobol sequence
图2 随机序列初始化种群图Fig. 2 Population initialization of random sequence
由式(8)~式(10)可知,蜜獾在更新代理位置的挖掘阶段与蜜蜂吸引阶段均受密度因子 ∂ 的影响,随着式(7)中迭代次数的增加,其易陷入局部最优的搜索状态。为增大搜索范围并提高搜索空间的历遍性,研究者尝试引入自适应权重、非线性控制参数调整策略等均衡全局搜索与局部搜索能力[22−25]。
改进密度因子 ∂1、∂2、∂3推导公式如下:
式中:r8为取值区间[0,1]的随机数。
HBA中密度因子 ∂ 迭代变化趋势较平缓,∂1、∂3前期迭代递减速度较快,而后期递减较慢,易产生局部搜索不完整的情况;∂、∂1迭代结束时未能收敛至坐标轴0点,可能导致局部搜索不完整;∂2前期迭代递减较缓,加强了全局搜索能力,而迭代中期递减速度加快,提高了收敛速度与跳出局部最优的能力,迭代末期则递减速率放缓,加强了局部搜索能力(图3)。为了平衡HBA的全局搜索与局部搜索能力,拟采用∂2作为改进后的密度因子。
图3 不同密度因子迭代变化曲线Fig. 3 Iterative variation curves of different density factors
莱维飞行是一种特殊的随机游走策略,其游走步长分布为服从重尾特征的概率分布[26]。在搜索空间中莱维飞行呈现频繁短距离小范围搜索与少量长距离大范围开发状态,该特性使其更易于帮助常规优化算法跳出局部最优。
为了进一步提升HBA跳出局部最优的能力,拟针对式(8)、式(10)迭代后产生的解进行莱维飞行策略扰动,扰动后新求解公式如下:
式中:L(D)为莱维飞行分布策略。
多策略优化HBA-XGBoost预测模型各阶段流程(图4)共包括以下6个步骤:①设定HBA种群数量、最大迭代次数、对应XGBoost模型超参寻优的维数及搜索空间上下限。②引入Sobol序列初始化HBA种群粒子,计算初始适应度函数与种群位置。③定义密度因子并引入∂2作为改进密度因子函数,定义搜索方向F取值区间。④基于随机数r6在[0, 1]区间内的取值确定F取值,按照HBA中逻辑判断将G取值代入式(8)或式(10)进行计算,获得当前迭代个体位置。⑤引入莱维飞行策略扰动当前解并得到新解,对比扰动前后适应度函数变化,择优对比记录。⑥满足设定最大迭代次数后,将最优超参组合赋予XGBoost预测模型,并对预测结果进行误差分析。
图4 多策略融合HBA-XGBoost预测模型流程图Fig. 4 Flow chart of HBA-XGBoost forecast model with multi-strategy fusion
研究对象为中国天然气月度需求量数据,历史数据主要来源于国家发改委、国家统计局、国家能源局、交通运输部、中国物流与采购联合会、中国汽车工业协会、中国碳核算数据库等多家单位官网发布的统计数据,时间跨度为2016年1月—2023年12月,共计96个月。针对天然气月度需求量预测研究所涉影响因素特点,重点分析相关领域研究成果中影响因素指标构建与筛选方式,综合考虑拟将GDP、进出口总额当期值、工业生产者出厂价格指数、工业增加值同比增长、居民消费价格指数、制造业采购经理指数、一次性能源消费总量(以吨标准煤当量tce计算)、电力消费量、煤炭消费量、二氧化碳排放量、环境污染治理投资总额、天然气产量、天然气出口量、天然气进口量、天然气重卡终端累计保有量、公路货运量16项指标作为影响因素进行分析。
由天然气消费总量年度分类统计数据可知,工业、制造业及电力热力企业天然气消费占比位居前3,交通运输业消费占比逐年提升。受国家大力发展电动乘用汽车政策与内河船舶“油改气”补贴政策结束的影响,天然气动力乘用汽车及内河“油改气”双燃料动力船舶领域的天然气需求量暂不纳入考虑范畴。受国际油价近年大幅波动影响,在中国商用车领域,天然气重卡已成为老旧重卡汰旧换新的主要选择,其保有量逐年递增,因此拟将天然气重卡终端累计保有量、公路货运量两类数据作为候选影响因素。此外,人口数量属年度统计数据且人口变动数据统计周期较长,暂不纳入特征序列筛选范畴。
采用MF-DFA法分析天然气月度需求量数据序列,其总体呈逐步递增态势。为了验证序列是否具有多重分形特征与长程相关性,采用最小二乘法m(m≤3)阶多项式进行拟合。随分形阶数q值(q∈[-10,10])变化的广义Hurst指数h(q)的变化曲线(图6)表明,当q取2时,h(q)分别为0.899 4、0.850 4、1.334 8;当q≠2时,h(q)随q的增大而逐渐减小,但均大于0.5。因此研究对象序列具备长程正向相关性且可以用于预测下一阶段变化趋势[27]。天然气月度需求量数据序列的多重分形谱(图7)为单峰函数,当q取2时,不同阶数的多重分形谱宽度 ∆ɑ 分别为 1.599 9、 1.347 3、1.213 7,而多重分形谱最大概率子集的分形维与最小概率子集的分形维间的差异∆f(ɑ)分别为 1.368 3、1.040 7、0.893 8;由∆f(ɑ)>0可知,研究对象序列波动率处于最大值的机会较处于最小值的机会占优,进一步说明研究对象序列具备多重分形特征。
图6 广义Hurst指数随q的变化曲线Fig. 6 Variation curve of generalized Hurst index with respect to q
图7 天然气月度需求数据序列的多重分形谱函数曲线Fig. 7 Multifractal spectrum function curve for monthly natural gas demand data sequences
上述影响因素特征序列中部分数据统计时间粒度不同,且个别序列存在数据缺失的情况(主要集中于1月、2月),考虑到模型输入数据质量与完整度对模型预测精度的影响,对特征序列数据集进行预处理。
GDP、二氧化碳排放量、环境污染治理投资总额等特征序列数据时间粒度通常为年度、季度。为统一数据时间粒度,需要对年度、季度数据进行低频向高频转换,即将其全部转换为月度数据。参考同类研究文献转换方式,拟采用EViews10软件,同时选择二次插值选项进行高频数据填充。
实际应用中,大多数模型不能处理包含缺失值的数据集,因此需要对缺失值进行插补。常用的插补方法包括均值插补、中位数插补(异常值处理)、众数插补、分段3次样条插值、拉格朗日插值、模型预测缺失值插补等。
电力消费量与煤炭消费量特征序列数据缺失率最大,其余存在数据缺失的特征序列数据缺失率为2%~8%(图8)。RF缺失值插补属于模型预测缺失值插补法,虽然其在处理较大数据集时耗时较长,但能较好地捕捉数据间的非线性关系,且对异常值具有一定的鲁棒性。以天然气产量缺失数据插补(图9)为例,将所有包含缺失值的特征序列都采用 RF法进行插补。
图8 天然气需求预测各特征序列数据缺失比例图Fig. 8 Proportions of missing data in various feature sequences of natural gas demand forecasts
图9 天然气产量缺失值随机森林插补图Fig. 9 Random forest interpolation for missing values of natural gas production
数据预处理环节对初始特征序列数据进行低频转换与插补处理,分别采用Boruta、SHAP、BorutaShap 3种方法进行特征序列筛选。经Boruta筛选后选择等级排序位列1级的二氧化碳排放量、天然气进口量、天然气产量、煤炭消费量、工业生产者出厂价格指数、进出口总值当期值、一次性能源消费总量、天然气重卡终端累计保有量8项重要特征与等级排序位列2级的1项待定特征(GDP)构建新特征序列集合;经SHAP筛选后考虑特征序列影响因素对测试集的影响程度,根据SHAP值大小拟选取天然气重卡终端累计保有量、天然气产量、天然气进口量、煤炭消费量、进出口总值当期值、GDP、居民消费价格指数、工业增加值同比增长、电力消费量、环境污染治理投资总额10项作为筛选后的新特征序列集合;经BorutaShap筛选,选择11项关键特征与1项暂定特征构建新的特征序列集合(图 10)。
理论上,模型输入特征序列维度越高越有助于提升预测精度,但各特征序列的影响程度与影响方向却有所不同,通过特征筛选可剔除具有干扰的特征序列,最大程度保留初始数据集的关键特征信息。对比不同特征序列组合(表1),方案2较方案1的预测精度有所提升;对比Boruta、SHAP及BorutaShap 3种方式特征筛选后的新特征序列数据集,方案5构建的新特征序列数据集预测精度最佳。因此,模型拟选择BorutaShap法对特征序列进行降维。
经数据预处理与BorutaShap特征筛选后的新特征数据集包含12个特征序列,每个序列包含96个月度数据,以2016—2022年(84个月)数据为训练集,2023年(12个月)数据作为测试集。为评估新模型的预测效果,拟 选取 XGBoost、 LSSVM(Least Squares Support Vector Machine)、 LSTM、 GRU(Gate Recurrent Unit)、RF、BP等非包含超参寻优算法的基础预测模型及近年新提出的包含 HBA、算数优化算法(Arithmetic Optimization Algorithm, AOA)、非洲秃鹰优化算法(African Vultures Optimization Algorithm, AVOA)、天鹰优化算法(Aquila Optimizer, AO)、金枪鱼优化算法(Tuna Swarm Optimization, TSO)等未采用改进策略优化的超参寻优算法XGBoost模型进行对比,选择MAE、MAPE及RMSE等对预测结果进行误差分析,采用决定系数R2评价模型的拟合效果。上述对比模型中,非包含超参寻优算法的基础预测模型采用默认参数设定,其中LSTM、GRU、BP中隐含层节点数采用经验公式√a+b+c (其中,a为模型输入数据维度,b为输出维度,c为[2,11]之间的整数)设定。此外,上述包含未采用改进策略优化的超参寻优算法XGBoost模型,除优化算法按默认参数设定外,其寻优维度均为3、种群规模为20、迭代次数为100次,决策树数量、决策树深度及学习速率的搜索范围分别为[1,1 000]、[1,100]、[0.000 1,1]。
图 10 BorutaShap筛选后新的特征序列图Fig. 10 New feature sequences after screening by BorutaShap
表1 不同方案下特征序列比选表Table 1 Comparison of feature sequences under different schemes
针对上述12种预测模型对天然气需求量进行预测误差分析(表2),该研究提出的多策略优化HBA-XGBoost模型预测误差最小, MAE、 MAPE、 RMSE、R2分别为为9.350 9、2.87%、11.335 3、0.890 9。此外,对所有未采用改进策略优化的超参寻优算法预测模型分别独立运行20次,发现各模型MAPE基本稳定在3%~4%范围内,进一步说明超参寻优能较好地提升模型预测精度,同时也克服了非超参寻优预测模型凭经验大量组合试算的不确定性。
表2 各模型预测结果评价指标列表Table 2 Evaluation indexes for forecast results from various models
对比超参寻优模型适应度变化(图 11)可知,新提出的多策略优化HBA在搜索区间内初始解覆盖范围、适应度函数下降速率等方面显著优于未采用改进策略的超参寻优算法,同时也能获取搜索迭代次数范围内的最佳超参组合,有效提升预测精度。
采用多策略优化HBA-XGBoost模型及11种对比模型分别对2023年1—12月天然气需求量进行预测(图 12),可见未采用改进策略优化的超参寻优算法预测模型预测值与实际值误差较小,而其他对比模型预测值则在3月、4月、10月上存在一定偏差。多策略优化HBA-XGBoost模型的 MAPE相较 HBA-XGBoost模型降低了0.19%,进一步提升了模型的预测精度。
图 11 超参寻优模型适应度变化曲线Fig. 11 Fitness variation curve of hyper-parameter optimization model
图 12 2023年天然气需求量实际值与各模型预测值误差对比图Fig. 12 Comparison of actual values and forecast results from various models of natural gas demand in 2023
采用MF-DFA与BorutaShap算法分别对研究对象可预测性及关键影响因素特征序列进行分析,构建了多策略优化HBA-XGBoost天然气需求预测模型,仿真计算结果证明其准确、有效。基于研究成果,得出下列结论:
1)多策略优化HBA-XGBoost天然气需求预测模型的预测精度优于对比预测模型。经Sobol低差异序列、改进密度因子及莱维飞行策略优化的HBA有效提升了XGBoost模型的超参寻优效能,进一步提升了模型预测精度,证明了新模型具有良好的可行性与有效性。
2)采用XGBoost模型分别对插值前后特征序列及经Boruta、SHAP、BorutaShap筛选后的新特征序列进行对比计算,结果显示采用BorutaShap算法可更好地剔除非关键特征因素并保留关键特征信息,不仅降低了模型特征输入维度,还提升了预测精度。
3)随着“双循环”与“双碳”战略的迅速发展,中国产业结构、能源结构及外部政策因素变化等都会对天然气需求产生深远影响。随着关键特征因素的不断变化,建立一套动态多维度特征筛选机制尤为重要。后续研究应着重于完善这一动态筛选机制,确保预测模型能够准确反映变化趋势,为天然气行业的可持续发展与能源政策的制定提供坚实的数据支撑。
主管单位:国家石油天然气管网集团
有限公司
主办单位:国家管网集团北方管道有
限责任公司
编辑出版:《油气储运》编辑部
通信地址:河北省廊坊市金光道51号
邮政编码:065000
主 编:张对红
副 主 编:关中原
电 话:0316-2176173
0316-2072055
传 真:0316-2072240
发行范围:公开发行
国内发行:河北省廊坊市邮政局
国外发行:中国出版对外贸易总公司
(北京782信箱)100011
国际刊号:ISSN 1000-8241
国内刊号:CN 13-1093/TE
邮发代号:18-89
广告许可:1310004000002