刘 炜 ，刘童童 ，王 辉 ，李鲲鹏 ，张 戬 ，桑国阳 ，吴拓剑

（1. 西南交通大学电气工程学院，四川 成都 611756；2. 广州地铁设计研究院股份有限公司，广东 广州 510010）

直流牵引供电计算的准确性对城市轨道供电系统设计、系统能耗评估、供电系统故障还原有着重要作用. 国内外研究人员对直流牵引供电系统的仿真计算做了大量研究，研究内容包含列车运行图编制优化、列车运行过程模拟和负荷过程仿真计算等.

在运行图编制方面，彭其渊等[1-2]考虑到发线运行方案，建立了在单、复线下均适用的列车运行图优化模型. 在列车运行过程模拟方面，主要是针对单列车准点节能运行操纵工况求解，文献[3-5]结合列车时刻表，提出两阶段寻优方法，分别优化列车区间运行操纵策略和站间运行时分分配方案，但计算量大，求解效率不高；文献[6-8]研究了在连续坡道和限速条件下的列车准点节能运行；文献[9-10]在研究列车控制工况最优切换时机时引入伴随变量，针对不同运行情况给出了最优切换规则. 而在城轨交通供电系统负荷过程仿真计算方面，文献[11]中建立了地面储能装置的通用模型，给出了不同工作状态下牵引供电系统等效电路的求解方法；文献[12]中考虑牵引所多运行状态，设计了城轨供电系统交-直-交交替迭代潮流求解算法. 以上研究均是以平铺运行图为核心对直流牵引供电系统的负荷过程进行仿真分析.

在供电系统设计阶段多以单列车运行曲线为基础，按照发车间隔生成平铺运行图，在供电计算运行图截面法的基础上模拟多列车运行，计算牵引网网压及钢轨电位等. 而在运营阶段，各列车按照时刻表运行[13]，考虑动态客流量的运营需求，须采用不同的行车运行图，如大小交路、快慢车等. 不同时段的列车运行区间、区间运行时间、停站时间均不尽相同[14-15]，如仍采用平铺运行图进行运行仿真则不能完全还原供电系统负荷的实际运行过程，也不能基于实际运行情况评估和优化再生制动能量的利用效果及分析现场频繁发生的钢轨电位异常问题，使得仿真结果与实际负荷过程差别大.

多列车运行情况下，供电系统仿真更多的是为了再现系统最高、最低网压，钢轨电位峰值，变电所负荷功率峰值和高峰小时有效值. 这些都与利用列车的真实站间运行时分还原多车运行过程直接相关，供电系统电气状态的极限值与车辆的启动、制动过程，以及各过程之间的相互配合亦均有直接关系.

本文要解决的核心问题是通过列车站间运行时分，还原列车运行轨迹，提高仿真模型准确性，更加准确地模拟供电系统负荷的实际运行过程. 基于此，本文以运行图中各列车站间运行时分为约束条件，建立列车定时节能运行的指标函数，基于固定阶梯级目标速度搜索算法优化列车操纵序列，得到多列车具有电气信息的运行轨迹. 以实迹运行图为驱动，实现还原多车运行轨迹的直流牵引供电系统负荷过程动态模拟. 实际案例分析表明，本文所提算法仿真牵引所负荷过程曲线与实测数据趋势相近，能针对供电系统特定问题还原系统运行过程，验证了模型及算法的有效性与可靠性.

1 以实迹运行图为驱动的供电系统仿真模型

1.1 运行图模型

运行图是列车在区间运行、在车站到发通停时刻的图解. 根据运营部门行车组织提供的列车时刻表，可绘制基于行车计划的运行图，再根据列车实际到发时刻，得到运营阶段的实迹运行图. 运行图简化模型如图1 所示. 以列车i在区间k运行为例，图1中：k1、k2分别为区间k的起始、终止车站节点；Di_k1、Di_k2分别为列车i到达区间k车站节点k1、k2的时刻；Fi_k1为列车i离开区间k车站节点k1的时刻；Ti_k为列车i在区间k的计划运行时间；j为站台编号，N为全线总站台数，满足j= 1，2， · ··，N.

图1 运行图简化示意Fig. 1 Simplified train operation diagram

运行图简化模型可用Ti_k和列车i在车站节点k1的停站时间Si_k1分别描述为

式中：i∈M，M为全线路运行的所有列车集合；k∈Li，Li为列车i在各区间运行时对应的弧线集合；k1，k2∈Ni，Ni为列车i经过的车站节点集合.

1.2 考虑惰行控制系数的列车定时节能控制模型

1.2.1 列车运行过程划分

研究表明，在列车制动降速前采用惰行工况是有效的节能操作，而最优惰行点的求解是节能控制的关键[16-17]. 以列车i在区间k的运行过程为例，为求解最优惰行点，引入目标速度vcm_i_k（m为目标速度所处的层级数）和惰行控制系数τi_k，重新将列车运行过程依次划分为3 个阶段，如图2 所示，阶段Ⅰ为列车速度由0 加速运行至vcm_i_k的启动过程，阶段Ⅱ为中间调速过程，阶段Ⅲ为制动降速过程，该阶段由速度曲线和考虑τi_k之后的牵引反算曲线交点决定. 图中：v为速度；x为位置；b为列车牵引反算速度曲线与线路限速的交点；b1为列车实际运行过程中速度曲线与牵引反算曲线的交点，交点位置为xi_k_b1；b2为在阶段划分过程中引入τi_k后列车速度曲线与牵引反算曲线的交点，交点位置为xi_k(τi_k，vcm_i_k)；P 表示牵引；C 表示惰行，B 表示制动.

图2 列车运行过程划分Fig. 2 Division of train operation process

根据以上定义，则τi_k可表示为

为保证安全运行，根据工程经验τi_k的取值范围为0≤τi_k≤0.25.

考虑vcm_i_k和τi_k划分列车运行区间，以列车在区间准点节能运行为目标，以运行图中各列车的运行时分为主要约束条件，建立各列车定时节能运行控制模型.

1.2.2 目标函数

定时层的目标是列车i在第k个区间以τi_k、vcm_i_k运行时的实际运行时间Tsi_k(τi_k，vcm_i_k)与Ti_k的差值δi_k(τi_k，vcm_i_k)最小，即

列车i在区间k以τi_k、vcm_i_k运行时时刻t的牵引功率可表示为

其中：μi_k(·)为牵引力使用系数，0 ≤μi_k(·)≤ 1；ηi为机电效率；vi_k(·)为速度；fmaxi_k(·)为最大牵引力.

列车i在区间k以τi_k、vcm_i_k运行时牵引能耗为

选择在满足式（4）定时目标的多组组合变量(τi_k，vcm_i_k)下使Ei_k最小的一组变量作为最终的优化结果，则列车区间运行的最终目标为

1.2.3 约束条件

为求解上述目标函数，需满足运行图中各列车在各车站到发时间、起停车速度、位置等边界条件约束；同时为保证安全运行，还需满足线路限速、乘客舒适度指标约束.

1） 边界约束

① 到发时间约束

式中：δmax为Tsi_k(τi_k，vcm_i_k)与Ti_k之间允许的最大误差.

② 起停车速度、位置约束

列车i在区间k的起点k1和终点k2分别通过各自的标志bi_k1和bi_k2的值确定是否停站（bi_k1，bi_k2= 1 表示停站，bi_k1，bi_k2= 0 表示不停站），记节点k1、k2所处的位置分别为Xi_k1、Xi_k2；vlim(·)为节点位置处的限速；vi和xi分别为列车i的速度与位置.

若bi_k1= 1 且bi_k2= 1，则边界约束可表示为

若bi_k1= 0 且bi_k2= 1，则边界约束可表示为

若bi_k1= 1 且bi_k2= 0，则边界约束可表示为

若bi_k1= 0 且bi_k2= 0，则边界约束可表示为

2） 线路限速约束

3） 乘客舒适度指标约束

列车运行过程中，加速度变化过大或工况切换过于频繁都会影响乘客舒适度[18]. 本文以加速度变化率作为评价乘客舒适度指标.(t)

式中： 、Δai_k(t)分别为列车i在区间k时刻t的加速度变化率、加速度变化量；ai_k(t)为列车i在区间k时刻t的加速度；Δt为仿真时间间隔.

1.3 直流牵引供电系统模型

直流牵引供电系统主要由牵引网、储能装置、列车、整流机组等组成[11]. 考虑到潮流算法的求解效率和计算精度，牵引网采用三层地网模型，整流机组采用恒压源-内阻模型，列车采用功率源模型，再生制动能量利用装置基于外特性建模.

还原多列车运行轨迹得到时刻t各列车的位置xi_k(τi_k，vcm_i_k，t)、功率Pi_k(τi_k，vcm_i_k，t)等，根据直流牵引供电系统等效模型，若时刻t全线节点数目为n，则在直流侧构建时刻t节点电压方程为

式中：G为直流侧供电系统节点导纳矩阵，

其中，Gfg为节点f与节点g的互导纳（f ≠ g）或自导纳（f=g）；

U为节点电压矩阵，

其中，Uf(t)为时刻t节点f的电压；

I为节点电压矩阵，

其中，If(t)为时刻t节点f的注入电流.

根据节点类型的不同，If(t)可表示为

式中：Ui_k(t)为时刻t列车i的电压；Us(t)、Rs(t)为时刻t牵引所s的输出电压、等效电阻.

2 模型求解及潮流计算

考虑惰行控制系数，以τi_k、vcm_i_k划分列车运行过程，基于固定阶梯级目标速度搜索算法完成各列车在各区间的准点节能运行，得到含有速度vi_k(t)、位置xi_k(t)、机车出力Fi_k(t)、机车取流Ii_k(t)等电气信息的实际运行轨迹，并用于城轨供电系统负荷过程仿真计算. 设：λm和σm分别为层级m对应阶梯值和最大搜索次数，m、λm和σm可按照线路条件及车辆条件具体按需设置；Δτ为惰行控制系数仿真步长.具体步骤如下：

步骤1加载列车、线路及运行图数据.

步骤2初始化仿真条件：令i= 1，k= 1，τi_k=0，m= 0，vcm_i_k=vmaxi_k，其中，vmaxi_k为列车i在区间k的最高运行速度.

步骤3利用vcm_i_k和线路中限速信息，根据文献[10]中的划分原则重新进行坡道划分，并根据τi_k、vcm_i_k划分列车区间运行过程，完成列车i在区间k的速度衔接.

步骤4计算Tsi_k(τi_k，vcm_i_k)，若Tsi_k(τi_k，vcm_i_k)＞Ti_k，则输出列车i在区间k的运行记录，k=k+ 1，返回步骤2；反之，令vcm_i_k=vaveri_k（其中，vaveri_k为列车i在区间k的平均速度），返回步骤3 完成该区间的速度衔接后，跳转至步骤5.

步骤5若Tsi_k(τi_k，vcm_i_k) ≠Ti_k，则采用固定阶梯级搜索算法，修改vcm_i_k，直至Tsi_k(τi_k，vcm_i_k) =Ti_k，具体搜索步骤为

① 初 始 化vcm_i_k=vaveri_k，Tsi_k(τi_k，vcm_i_k) =Tsi_k(τi_k，vaver_i_k)，层级m搜索次数dm= 1；

② 计算当目标速度为vcm_i_k+dmλm的时间Tsi_k(τi_k，vcm_i_k+dmλm)，若Tsi_k(τi_k，vcm_i_k+dmλm) =Ti_k，输出vcm_i_k=vcm_i_k+dmλm，转入④；

③ 若Tsi_k(τi_k，vcm_i_k+dmλm) ＞Ti_k＞Tsi_k(τi_k，vcm_i_k)，则Tsi_k(τi_k，vcm_i_k) =Tsi_k(τi_k，vcm_i_k+dmλm)，m=m+1，dm= 1，转入②，否则dm=dm+ 1，Tsi_k(τi_k，vcm_i_k) =Tsi_k(τi_k，vcm_i_k+dmλm)，转 入②；若dm＞σm，则m=m-1，dm= 1，转入②；

④ 保存列车区间运行记录Ri_k(τi_k，vcm_i_k).

步骤6令τi_k=τi_k+ Δτ，其中， Δτ为惰性控制系数步长. 若τi_k≤0.25，则令vcm_i_k=vmaxi_k返回至步骤3；反之，则选择牵引能耗值最小的Ri_k(τi_k、vcm_i_k)为列车i在区间k的最终优化结果.

步骤7令k=k+ 1，若k＜Li，则令τi_k= 0，m=0，vcm_i_k=vmaxi_k，转入步骤3；反之，转入步骤8.

步骤8令i=i+ 1，若i ＞ M，则转入步骤9；反之，令k= 1，若Ti_k=Tz_k（z为列车编号，z= 1,2,…,i-1），则Ri_k=Rz_k(τz_k，vcm_z_k)，并转入步骤7；反之，令τi_k= 0，m= 0，vcm_i_k=vmaxi_k，转入步骤3.

步骤9根据各列车优化后的运行记录，还原多列车具有电气信息的运行轨迹.

步骤10初始化迭代次数r、电压收敛精度ε、变电所状态W，仿真起止时间及仿真时长T，直流侧供电系统节点导纳矩阵G.

步骤11读取当前时刻t下各列车的位置及功率信息.

步骤12计算时刻t直流牵引供电系统潮流，根据式（15）更新各节点电压矩阵U、节点电流矩阵I.

步骤13设Ur为第r次迭代时的电压，判断是否满足 |Ur-Ur-1| ＜ε，若不满足则令r=r+ 1，返回步骤12；否则，转入步骤14.

步骤14判断牵引所状态W是否需要调整，若W不合理，则调整后返回步骤13；否则，令t=t+Δt，转入步骤15.

步骤15若t＜T，则返回步骤11；否则，输出仿真结果.

3 实例分析

3.1 算例1

以某运营线路作为仿真实例，该线路全长25.442 km，牵引所1、9、10 安装有逆变回馈装置，采用自动驾驶系统（automatic train operation，ATO）运行，仿真车辆为4 动2 拖6 编组车辆，定员载荷下车重为291.8 t，冲击限制为0.75 m/s3，结构速度为115.00 km/h，最大加速度为1.2 m/s2. 牵引所和车站位置分布如图3 所示，具体的位置信息见表1.

表1 车站位置信息Tab. 1 Information of station locations

图3 牵引所与车站位置分布Fig. 3 Locations of traction stations and stations

3.1.1 单列车运行过程分析

以第5 列车在区间4 运行过程为例验证基于惰行控制系数的固定阶梯级目标速度搜索算法. 该区间Ti_k为152 s，仿真过程中τi_k、vcm_i_k、xi_k(τi_k，vcm_i_k)详细变化关系如图4 所示. 由图4 可知，vcm_i_k随着τi_k的增加而增加，当τi_k大于0.12 之后，不能满足定时目标，故该区间最终的优化结果：τi_k= 0.12，vcm_i_k= 96.95 km/h，xi_k(τi_k，vcm_i_k) = 10.420 km.

图4 仿真结果Fig. 4 Simulation results

3 种算法区间运行过程结果对比如图5 所示，其中：CaseA1 为实测数据，CaseA2 为文献[16]中所提算法，CaseA3 为本文算法，分别用实线、短划线、点划线表示；图中：蓝色表示牵引工况，绿色表示惰行工况，红色表示惰行工况；线路坡道为实际坡度数据，仿真坡道为根据目标速度vcm_i_k及线路限速重新进行坡道划分的结果；o1、o2、o3分别为CaseA1、CaseA2、CaseA3 加速至目标速度的工况装换点；v1、v2、v3为CaseA1、CaseA2、CaseA3 的转换速度；p1、p2分别为CaseA1、CaseA2 在中间运行过程中的工况转换点.

由图5 可知：CaseA3 在列车启动阶段加速运行至vcm_i_k，工况由牵引转为惰行，转换点为o3，目标速度vcm_i_k=v3，惰行与制动的转换点为q3，xi_k（τi_k，vcm_i_k）为10.420 km；CaseA2 列车制动位置xi_k_2为9.590 km，转换点为q2；CaseA1 列车制动位置xi_k_1为10.330 km，转换点为q1.

图5 区间运行过程对比Fig. 5 Comparison of running sections

全线有8 个区间，CaseA1 总能耗为332.590 kW·h，CaseA2 总能耗为350.870 kW·h，CaseA3 总能耗为331.470 kW·h，3 种算法车公里能耗和吨公里能耗比较如表2 所示.

表2 车公里能耗与吨公里能耗实测与仿真对比Tab. 2 Comparison of measured and simulated energy consumption of trains per km and per ton-km kW·h

由表2 可知，CaseA3 仿真能耗与CaseA1 实测能耗基本吻合，验证了所提算法模拟列车运行过程的有效性与合理性，与CaseA2 相比全线节能效果提升约5.53%.

3.1.2 负荷过程仿真分析

选择早发车时段 0 5：30：00—06：30：00 、高峰小时 0 8：00：00—09：00：00 分别进行负荷过程的仿真分析. 以早发车时段为例，运行图对比如图6 所示，图中横轴为站台编号，纵轴为时间.

图6 运行图对比Fig. 6 Operation diagram comparison

课题组测量数据为牵引所1 和牵引所9 的负荷过程，牵引所监测点位置安装见文献[12]，以牵引所9 高峰小时仿真结果为例，整流机组电流仿真与实测负荷数据比较如图7 所示，逆变回馈装置电流仿真与实测数据比较如图8 所示.

图7 整流机组电流仿真与实测曲线对比Fig. 7 Comparison of current simulation and measured curve of rectifier unit

图8 逆变回馈装置电流仿真与实测曲线对比Fig. 8 Comparison of current simulation and measured curve of inverter feedback device

由图7、图8 可知：基于实迹运行图的牵引变电所负荷过程仿真结果与实测数据相比，负荷过程趋势相似，经计算 0 5：30：00—06：30：00、08：00：00 —09：00：00负荷过程仿真与实测数据Person 系数均在0.89 以上，仿真与实测数据相关性较强，本文模型及算法的有效性均得到验证.

同时对高峰小时 0 8：00：00—09：00：00 分别采用平铺运行图（CaseB2）、实迹运行图（CaseB3）仿真分析，并与实测数据（CaseB1）进行比较，1 h 数据统计见表3.

表3 高峰小时负荷过程统计表Tab. 3 Statistics of load process at peak hour

表3 中，装置节能率用每小时逆变回馈装置反馈能量与整流机组牵引能耗的比值表示，CaseB2 仿真结果与CaseB1 最大误差达17.05%，CaseB3 仿真结果与CaseB1 误差最大不超过6.85%，较CaseB2仿真结果准确度最高可提升12.91%，故采用实迹运行图更能还原供电系统负荷实际运行情况，较平铺运行图更具有参考价值.

3.2 算例2

3.2.1 工程概况

某地铁线路全长42.600 km，采用6 动2 拖8 编组车辆，最高时速100.00 km/h，最小发车间隔2.5 min，1 500.00 V 接触网供电. 牵引所及车站位置分布信息如图9 所示，全线共计车站27 座，牵引所20 座，其中车站4、8、11、13、15、17、26 位置处无牵引所. 经统计，该线路某月钢轨电位限制装置闭锁总次数达257 次，部分车站闭锁次数高达31 次.

图9 牵引所与车站位置分布Fig. 9 Locations of traction stations and stations

课题组为该线路的钢轨电位异常问题进行大量实测，选取该线路轨电位异常情况较为突出的车站8 及邻近区间进行监测. 车站8 位于9.43 km 处，钢轨电位信号取自钢轨电位限制装置（over voltage protection device，OVPD）的负母排和地母排之间，采用16 通道同步采集装置记录数据，监测时间不小于24 h.

3.2.2 实测结果及仿真分析

经测量及统计，全日钢轨电位最大值达118.16 V（ 0 8：09：19 ），最小值为-68.72 V （ 0 6：40：33 ），早高峰 0 8：09：19 该车站OVPD 动作，动作电压120.00 V，之后OVPD 持续闭合.

根据该线路供电系统图及车站位置信息，搭建供电系统模型，对高峰小时采用本文算法还原实迹运行图并进行负荷过程仿真分析，钢轨电位仿真与实测数据比较如图10 所示. 图中，钢轨电位的仿真曲线与实测过程趋势接近，在峰值阶段比中间过程吻合的情况好. 主要原因是：还原列车的实际运行过程中，不同类型列车在区间运行过程有不同的驾驶策略，本文以运行图为驱动设计的定时节能算法还原列车运行轨迹与实际列车驾驶过程仍然存在一定差异，若需得到不同驾驶策略下供电系统仿真结果，可通过采用不同的目标速度或设置不同的区间运行时间进行模拟.

由图10 可知：车站8 OVPD 在 0 8：09：19 因钢轨电位瞬时值超过120.00 V 而动作. 基于实迹运行图供电仿真截面法追踪至该时间切面，钢轨电位仿真峰值为116.60 V，该车站临近区间（5.780 ～ 10.780 km）上行有2 辆车牵引，1 辆车惰行，下行有3 辆车牵引，1 辆车制动，瞬时需求功率达11 903 kW，在08：09：19的轨电位动作过程与该时刻该牵引所临近区间多列车集中牵引取流有着直接关联.

图10 钢轨电位仿真与实测数据对比Fig. 10 Comparison of simulated and measured rail potential

4 结 论

为提高仿真模型的准确性，更准确地模拟供电系统负荷的实际运行过程，本文以实迹运行图为驱动，建立了多列车定时节能控制模型，通过目标速度搜索算法还原各列车具有电气信息的运行轨迹，实现了城轨供电系统负荷过程动态仿真分析，得到以下结论：

1） 考虑惰行控制系数的固定阶梯级目标速度搜索算法，能有效完成列车在区间准点节能运行，与已有算法相比节能效果可提升约5.53%.

2） 基于运行图的城市轨道交通供电系统负荷过程仿真，在早晚收发车时段与高峰小时，仿真结果与实测数据的Person 系数在0.89 以上，负荷过程特征值仿真与实测数据误差最大不超过6.85%，较平铺运行图仿真结果准确度最高可提升12.91%.

3） 在运营阶段，利用实迹运行图还原多列车运行过程和供电系统负荷过程对评估和优化再生制动能量利用效果、还原供电系统异常情形及分析异常产生原因、制定解决方案等有着十分重要的作用. 若需完成全天或较长时段负荷过程的仿真计算，算法效率对硬件性能有更高的要求，为提高算法效率，下一步考虑供电仿真过程的并行加速.