CN114912068B - 一种dc牵引网络的潮流分析方法 - Google Patents

一种dc牵引网络的潮流分析方法 Download PDF

Info

Publication number
CN114912068B
CN114912068B CN202210614934.9A CN202210614934A CN114912068B CN 114912068 B CN114912068 B CN 114912068B CN 202210614934 A CN202210614934 A CN 202210614934A CN 114912068 B CN114912068 B CN 114912068B
Authority
CN
China
Prior art keywords
train
equation
voltage
node
power flow
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210614934.9A
Other languages
English (en)
Other versions
CN114912068A (zh
Inventor
吕闻龙
蔡亮成
黄德青
陈鹏豪
冯灏钰
肖添
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN202210614934.9A priority Critical patent/CN114912068B/zh
Publication of CN114912068A publication Critical patent/CN114912068A/zh
Application granted granted Critical
Publication of CN114912068B publication Critical patent/CN114912068B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/30—Circuit design
    • G06F30/32—Circuit design at the digital level
    • G06F30/33—Design verification, e.g. functional simulation or model checking
    • G06F30/3323—Design verification, e.g. functional simulation or model checking using formal methods, e.g. equivalence checking or property checking
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Electricity, gas or water supply
    • H—ELECTRICITY
    • H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J1/00—Circuit arrangements for dc mains or dc distribution networks
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
    • Y02E60/60—Arrangements for transfer of electric power between AC networks or generators via a high voltage DC link [HVCD]

Abstract

本发明公开了一种DC牵引网络的潮流分析方法，具体为：通过结合并重构潮流方程得到一个只需给定
Figure DDA0003673878890000011
和
Figure DDA0003673878890000012
可递推地计算得到U的方程，再提出了一个等效的fc的计算图，以链式计算代替了矩阵运算，不断迭代完成潮流分析。本发明具有二次的收敛速度和线性的单轮迭代计算复杂度。

Description

一种DC牵引网络的潮流分析方法
技术领域
本发明属于DC牵引网络中的潮流分析技术领域，尤其涉及一种DC牵引网络的潮流分析方法，
背景技术
现代地铁系统由直流牵引系统驱动，该系统应用大量结点和电路来供给和传输电力。为了获得安全、经济的电力输送过程，潮流分析是获取和预测电压、电流和负载值的有效途径，其关键在于求解网络的潮流方程。由于牵引系统非常庞大且复杂，因此需要高效且快速收敛的潮流方程求解方法。
为此，DC潮流方程的求解方法已被广泛研究，其中Gauss-Seidel方法[1]是广为人知且广泛使用的方法，但是它需要大量的迭代。为减少迭代次数，Newton-Raphson方法[2]具有二次的收敛速度，但是它在每轮迭代中都需要额外地计算Jacobian矩阵及其逆矩阵。同时，一种基于LU分解的方法[3]也可以有效地减少潮流计算的迭代次数。
此外，在不使用Jacobian矩阵的情况下，[4]使用连续逼近方法可以减少计算负担，并且引入了充分的条件以保证解的唯一性和方法的收敛性。泰勒级数展开被用于将潮流计算问题线性化[5]，它可以有效降低计算量。还有一些方法[6-8]可高效地求解具有径向结构的网络的潮流，但是它们要求网络中只有一个独立电压源，因此无法适用于DC牵引网络。
即使上述方法减少了每轮迭代所需的计算，但由于它们简化了潮流方程且没有利用到梯度信息，它们的收敛性能会变差。据作者所知，对于一般的DC潮流计算问题，现有的方法难以同时具有高效的收敛速度和计算复杂度。但通过利用牵引网络的贯通式结构特性，可将整个网络分解为多个独立的梯形网络并分别求解各梯形网络的潮流[8]，从而减小了问题的规模。对于每个梯形网络，Jacobian矩阵可通过链式法则表示为一系列微分的乘积，但是该方法只适用于两辆列车同向行驶的网络。另外，对于两辆列车相向行驶的网络可以无需迭代地求解它的潮流[9]。但是，一般的DC牵引网络要求其上下行轨道上各有n辆和m辆列车，此时这两种方法不再适用。
参考文献：
[1]C.Li,S.K.Chaudhary,T.Dragiˇcevi′c,J.C.Vasquez,and J.M.Guerrero,“Power flow analysis for dc voltage droop controlled dc microgrids,”2014 IEEE11th International Multi-Conference on Systems,Signals&Devices(SSD14),pp.1–5,2014.
[2]C.Liu,B.Zhang,Y.Hou,F.F.Wu,and Y.Liu,“An improved approach for ac-dc power flow calculation with multi-infeed dc systems,”IEEE Transactions onPower Systems,vol.26,pp.862–869,2011.
[3]艾芊,电力系统稳态分析.清华大学出版社,2014.
[4]A.Garc′es,“Uniqueness ofthe power flow solutions in low voltagedirect current grids,”Electric Power Systems Research,vol.151,pp.149–153,2017.
[5]O.D.Montoya,L.F.Grisales-Nore～na,D.Gonzalez-Montoya,C.A.Ramos-Paja,and A.Garc′es,“Linear power flow formulation for lowvoltage dc powergrids,”Electric Power Systems Research,2018.
[6]O.D.Montoya,L.F.Grisales-Nore～na,and W.J.Gil-Gonz′alez,“Triangular matrix formulation for power flow analysis in radial dc resistivegrids with cpls,”IEEE Transactions on Circuits and Systems II:Express Briefs,vol.67,pp.1094–1098,2020.
[7]L.F.Grisales-Nore～na,O.D.Garzon-Rivera,C.A.Ram′1rez-Vanegas,O.D.Montoya,and C.A.Ramos-Paja,“Application ofthe backward/
forward sweep method for solving the power flow problem in DCnetworks with radial structure,”Journal ofPhysics:Conference Series,vol.1448,p.012012,jan 2020.
[8]O.D.Montoya,“On the existence ofthe power flow solution in dcgrids with cpls through a graph-based method,”IEEE Transactions on Circuitsand Systems II:Express Briefs,vol.67,pp.1434–1438,2020.
[9]W.J.Davis,The tractive resistance of electric locomotives andcars.General Electric,1926.
[10]B.P.Rochard and F.Schmid,“A review of methods to measure andcalculate train resistances,”Proceedings of the Institution of MechanicalEngineers,Part F:Journal of Rail and RapidTransit,vol.214,pp.185–199,2000.
[11]J.W.Simpson-Porco,F.D¨orfler,and F.Bullo,“On resistive networksofconstant-power devices,”IEEE Transactions on Circuits and Systems II:Express Briefs,vol.62,pp.811–815,2015.
[12]Z.Liu,X.Zhang,M.Su,Y.Sun,H.Han,and P.Wang,“Convergence analysisof newton-raphson method in feasible power-flow for dc network,”IEEETransactions on Power Systems,vol.35,pp.4100–4103,2020.
[13]X.Yang,X.Li,B.Ning,and T.Tang,“An optimisation method for trainscheduling with minimum energy consumption and travel time in metro railsystems,”Transportmetrica B:Transport Dynamics,vol.3,pp.79–98,2015.
[14]S.Su,T.Tang,and C.Roberts,“A cooperative train control model forenergy saving,”IEEE Transactions on Intelligent Transportation Systems,vol.16,pp.622–631,2015.
[15]Q.Wang,X.Feng,J.Zhu,and Z.Liang,“Simulation study on optimalenergy-efficient control ofhigh speed train considering regenerative brakeenergy,”China Railway Science,vol.36,pp.96–103,2015.
[16]Z.Tian,S.Hillmansen,C.Roberts,P.Weston,N.Zhao,L.Chen,and M.Chen,“Energy evaluation of the power network of a dc railway system withregenerating trains,”IET electrical systems in transportation,vol.6,pp.41–49,2016.
发明内容
为了实现高效地求解DC牵引网络的潮流方程，本发明提供一种DC牵引网络的潮流分析方法。
本发明的一种DC牵引网络的潮流分析方法，包括以下步骤：
步骤1：构建DC牵引网络模型。
采用第三轨供电的DC牵引网络由变电站、悬链线、回程轨道和列车组成，其中列车运行在上行或下行方向的轨道；两个变电站sL和sR分别位于网络的前后端，其位置分别为xL和xR，在两个变电站之间，有n辆列车运行在上行方向轨道，其位置分别为
Figure BDA0003673878870000031
有m辆列车运行在下行方向轨道，其位置分别为
Figure BDA0003673878870000032
绘制DC牵引网络的等效电路，其中每辆列车等效为一个恒功率源，每个变电站等效为一个带有串联电阻Rs的恒定直流电压源，其电阻值和电压值分别为变电站的内阻和空载电压；等效电路中，
Figure BDA0003673878870000033
为对应列车的功率，由下式计算；
P(v)＝(Ft/ηt-αηbFb)v (4)
其中，0＜ηt,ηb＜1分别为列车在牵引模式和制动模式下的机电效率的有效值，0＜α＜1为列车制动能量的利用率；Ft和Fb分别为列车所受的牵引力和制动力，v为列车的速度。
Figure BDA0003673878870000034
为对应列车所处悬链线的电压，UL和UR为两个变电站的输出电压；
Figure BDA0003673878870000035
为对应的列车与列车之间或列车与变电站之间的悬链线和回程轨道的等效电导值，由下式计算：
Figure BDA0003673878870000041
其中，ρc和ρr分别为悬链线和回程轨道的电阻率。
步骤2：推导潮流方程。
DC牵引网络的等效电路的结点电压方程如下式：
Figure BDA0003673878870000042
式中，
Figure BDA0003673878870000043
等号右侧表示每个结点的输入电流，由Tellegen第一定理推导为：
Figure BDA0003673878870000044
若令
Figure BDA0003673878870000045
以及
Figure BDA0003673878870000046
则将式(6)表达为矩阵形式GU＝I，其中G表示导纳矩阵，由于各结点的位置已知，G中的每一项均可由式(5)得到。
步骤3：求解潮流方程。
将式(6)和(7)结合并重构得到：
Figure BDA0003673878870000047
其中，
Figure BDA0003673878870000048
当
Figure BDA0003673878870000049
和
Figure BDA00036738788700000410
给定时，U可被递推地计算得到；因此，将
Figure BDA00036738788700000411
和
Figure BDA00036738788700000412
视为决策变量，即定义
Figure BDA00036738788700000413
分析可知，U满足潮流方程(6)和(7)当且仅当Uc满足
Figure BDA00036738788700000414
即Uc为函数
Figure BDA00036738788700000415
的根；fc的Jacobian矩阵Jc由下式计算：
Figure BDA0003673878870000051
定义在计算图中有从结点wi指向wj的箭头，则wj称为wi的子结点，同时wi称为wj的父结点；
Figure BDA0003673878870000052
为wi的子结点集，对应地，wj的父结点集表示为
Figure BDA0003673878870000053
在每一次迭代中具有前向传播和反向传播两个步骤。
在前向传播中，使用式(8)计算fc(Uc)，同时计算每一个中间变量结点wj对其父结点集中每个结点的偏导，即
Figure BDA0003673878870000054
由下式给出：
Figure BDA0003673878870000055
在反向传播中，根据计算图以及微分的链式法则，fc对每个结点的偏导数可由式(11)计算，其计算顺序与前向传播过程中的相反。
Figure BDA0003673878870000056
由fc的定义，可得式(11)的初始条件为：
Figure BDA0003673878870000057
接着，根据式(9)即得到Jc，根据Newton-Raphson方法的计算框架，若给定初始电压向量
Figure BDA0003673878870000058
则第k次迭代的更新公式为：
Figure BDA0003673878870000059
最终，根据是否满足终止条件来决定将结束还是进入下一次迭代。
迭代终止条件为：
k≥kmax or||ΔU||≤(ΔU)tol
其中，kmax为最大迭代次数，(ΔU)tol为最小电压变化步长容量，||·||为向量的L2范数；
Figure BDA0003673878870000061
为电压变化步长。
本发明的有益技术效果为：
本发明具有二次的收敛速度和线性的单轮迭代计算复杂度。本发明提出的CNR方法是一个高效的针对DC牵引网络的潮流计算方法，它是NR方法的改进。基于潮流方程提出了一个等效的计算图，无论网络中的列车如何分布，CNR方法都只需要两个变量，并且以链式计算代替了矩阵运算。基于北京地铁亦庄线的实际数据进行的安利分析，CNR方法不仅将单轮迭代的计算复杂度降低到了O(s)，而且还具有二次的收敛速度。
附图说明
图1为DC牵引网络的物理模型。
图2为DC牵引网络的等效电路图。
图3为函数fc的计算图。
图4为CNR方法的伪代码。
图5为北京地铁亦庄线的拓扑结构图。
图6为列车时刻表的各时间点关系示意图。
图7为CNR方法中第一辆车的电压曲线。
图8为仿真各方法的收敛性能曲线图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步详细说明。
本发明的一种DC牵引网络的潮流分析方法，具体为：
步骤1：构建DC牵引网络模型。
单车的运行模型可基于牛顿第二定律描述如下：
Ma＝Ft-Fb-R0(v)-Rr(s) (1)
其中，M,s,v,a分别为列车的质量，位置，速度和加速度，Ft和Fb分别为列车所受的牵引力和制动力，可由列车实时运行状态得到。R0表示列车所受的基本运行阻力，它可通过(2)中的Davis公式[9]计算。Rr表示列车所受的坡度附加运行阻力，它通过式(3)计算：
R0(v)＝c0+c1v+c2v2 (2)
Rr(s)＝Mgsin(θ) (3)
最终，列车的功率可由式(4)得到。当列车处于牵引模式时，功率为正值；当列车处于制动模式时，功率为负值。
P(v)＝(Ft/ηt-αηbFb)v (4)
其中，0＜ηt,ηb＜1分别为列车在牵引模式和制动模式下的机电效率的有效值，0＜a＜1为列车制动能量的利用率，它由铁路系统中能量储存装置的电容等因素决定。
采用第三轨供电的DC牵引网络由变电站、悬链线、回程轨道和列车组成，其中列车运行在上行或下行方向的轨道。一般的DC牵引网络的物理模型如图1所示。图中两个变电站sL和sR分别位于网络的前后端，其位置分别为xL和xR，在两个变电站之间，有n辆列车运行在上行方向轨道，其位置分别为
Figure BDA0003673878870000071
有m辆列车运行在下行方向轨道，其位置分别为
Figure BDA0003673878870000072
绘制DC牵引网络的等效电路如图2所示，其中每辆列车等效为一个恒功率源，每个变电站等效为一个带有串联电阻Rs的恒定直流电压源，其电阻值和电压值分别为变电站的内阻和空载电压；等效电路中，
Figure BDA0003673878870000073
为对应列车的功率，它可由式(4)计算。
Figure BDA0003673878870000074
为对应列车所处悬链线的电压，UL和UR为两个变电站的输出电压；
Figure BDA0003673878870000075
为对应的列车与列车之间或列车与变电站之间的悬链线和回程轨道的等效电导值，由下式计算：
Figure BDA0003673878870000076
其中，ρc和ρr分别为悬链线和回程轨道的电阻率。此外，记
Figure BDA0003673878870000077
以简化公式的形式。
步骤2：推导潮流方程。
DC牵引网络的潮流分析需要在已知变电站参数(Us,Gs)和列车状态(位置和功率)的条件下求解电压UL,
Figure BDA0003673878870000078
通常使用结点电压法对网络潮流进行分析[11]，DC牵引网络的等效电路的结点电压方程如下式：
Figure BDA0003673878870000081
式中，
Figure BDA0003673878870000082
等号右侧表示每个结点的输入电流，由Tellegen第一定理推导为：
Figure BDA0003673878870000083
若令
Figure BDA0003673878870000084
以及
Figure BDA0003673878870000085
则将式(6)表达为矩阵形式GU＝I，其中G表示导纳矩阵，由于各结点的位置已知，G中的每一项均可由式(5)得到。此外，本发明使用s＝n+m+2表示潮流方程中变量的个数。
步骤3：求解潮流方程。
将式(6)和(7)结合并重构得到：
Figure BDA0003673878870000086
其中，
Figure BDA0003673878870000087
当
Figure BDA0003673878870000088
和
Figure BDA0003673878870000089
给定时，U可被递推地计算得到。具体的，通过式(8a)可计算得到UL，然后按照i＝1,2,...,n和j＝1,2,...,m的顺序利用式(8b)和(8c)，可依次计算得
Figure BDA00036738788700000810
和
Figure BDA00036738788700000811
最终利用式(8d)得到UR。到此，U的每一个分量都可被计算，只需给定
Figure BDA00036738788700000812
和
Figure BDA00036738788700000813
因此，将
Figure BDA00036738788700000814
和
Figure BDA00036738788700000815
视为决策变量，即定义
Figure BDA00036738788700000816
如前所述，U可由Uc计算得到，所以它们是等效的。分析可知，U满足潮流方程(6)和(7)当且仅当Uc满足
Figure BDA00036738788700000817
即Uc为函数
Figure BDA00036738788700000818
的根。根据前文所述的由
Figure BDA00036738788700000819
和
Figure BDA00036738788700000820
计算U的过程，可将函数fc的计算图绘制于图3。fc的Jacobian矩阵Jc由下式计算：
Figure BDA0003673878870000091
注意Jc是无法直接计算得到的，本发明使用了自动微分技术以计算Jc，为了讲述其过程，需要介绍两个重要的定义：
定义1：在计算图中有从结点wi指向wj的箭头，则wj称为wi的子结点，同时wi称为wj的父结点。需注意一个结点的子结点或父结点的数量可以大于1。
定义2：
Figure BDA0003673878870000092
为wi的子结点集，对应地，wj的父结点集表示为
Figure BDA0003673878870000093
在每一次迭代中具有前向传播和反向传播两个步骤。
在前向传播中，使用式(8)计算fc(Uc)，同时计算每一个中间变量结点wj对其父结点集中每个结点的偏导，即
Figure BDA0003673878870000094
由下式给出：
Figure BDA0003673878870000095
在反向传播中，根据计算图以及微分的链式法则，fc对每个结点的偏导数可由式(11)计算，其计算顺序与前向传播过程中的相反。
Figure BDA0003673878870000096
由fc的定义，可得式(11)的初始条件为：
Figure BDA0003673878870000097
接着，根据式(9)即得到Jc，根据Newton-Raphson方法的计算框架，若给定初始电压向量
Figure BDA0003673878870000101
则第k次迭代的更新公式为：
Figure BDA0003673878870000102
最终，根据是否满足终止条件(见如图4所示的伪代码)来决定将结束还是进入下一次迭代。
迭代终止条件为：
k≥kmax or||ΔU||≤(ΔU)tol
其中，kmax为最大迭代次数，(ΔU)tol为最小电压变化步长容量，||·||为向量的L2范数；
Figure BDA0003673878870000103
为电压变化步长。
本方法的主要计算量在于fc和Jc的计算，它们都是线性复杂度的。同时，本方法也具有二次的收敛速度。需要指出的是，要使CNR方法收敛到目标解需要一定的条件。首先，保证潮流方程(6)和(7)的解存在所需的条件在文献[4]中给出。在此基础上，可应用文献[12]提出的NR方法的收敛条件推导得CNR方法收敛的条件。
实施例
北京地铁亦庄线是一个典型的第三轨方式供电线路，其供电制式为DC-750V，因此可以用第二节中的模型对其描述。亦庄线共23.3km长，具有14座车站和6个变电站(即5个供电区间)，它的拓扑结构图绘制于图5。各站点的位置信息和列车运行时刻表如表1[13]。
表1列车运行时刻表
Figure BDA0003673878870000104
Figure BDA0003673878870000111
其中，从列车在当前站点启动时开始计时，0～t1时段内为牵引工况，此时只受牵引力；t1～t2时段内为巡航工况，此时合力为0；t2～t3时段内为制动工况，此时只受制动力。在t3时刻列车到达下一站点，停留至t4时刻时再次启动并重新计时。该过程的示意图绘制于图6。根据地铁系统的运行特点，假设每辆列车从起点站出发，沿上行轨道运行至终点站，然后折返并沿着下行轨道回到起点站，并不断循环。在亦庄线的高峰时段(17:30–18:30)，共14辆列车在轨道上循环地运行，其初始位置如表2所示。
表2各列车的初始位置
Figure BDA0003673878870000112
此外，亦庄线的坡度信息由表3给出[13]：
表3亦庄线的坡度信息表
位置 0 160 470 970 1370 1880 2500 2700 3170 3570
坡度 -2 -3 10.4 3 -8 3 -2 -3 8.2 2
位置 3940 4200 4800 5200 5800 6050 6370 6770 7150 7415
坡度 -20.4 -24 0 -2 -3.2 0 3.3 2.8 -15.6 9
位置 7675 8376 8736 9036 9366 9806 10206 10606 10866 11426
坡度 0 5 -2 0 -2 5 3 0 2 -3
位置 11826 12116 12736 13116 13526 13926 14546 15176 15476 16006
坡度 0 3.5 -1.8 0 -0.5 1.5 -1 6 0 -8
位置 16326 16696 17136 17816 18136 18486 19186 19426 19776 20121
坡度 -3 5 1.4 0 15.5 24 -3 10.1 2 -3
位置 20796 21231 21481 21681 22066 22416 22728
坡度 3 2 20 3 -18.9 2 \
其中位置的单位为(m)，坡度的单位为弧度的千分数。上表数据读取方式为：0m至160m处的坡度为-2×10-3，160m至470m处的坡度为-3×10-3，以此类推。根据坡度信息和式(3)，可得到列车所受的坡度附加阻力。
亦庄线使用DKZ32型电动车组，它的最高限速为80km/h。根据文献[14]，可将列车的基本运行阻力、制动力和牵引力特性列于式(14)–(16)。式中，v的单位为km/h，Ft,Fb,R0,Rr的单位为kN。
Figure BDA0003673878870000121
Figure BDA0003673878870000122
R0(v)＝2.965-0.064v+3.858×10-4v2 (16)
最终根据文献[15,16]，可将列车和变电站的参数列于表4中。
表4测试环境的详细参数表
Figure BDA0003673878870000123
根据表1、表3以及式(14)～(16)、式(3)可得每辆列车在任意时刻的受力情况，再根据表2给出的列车初始位置和列车运动方程(1)可得列车任意时刻的位置x和速度v.进而根据x和变电站位置(在图5中标出)可判断各列车所在供电区间的编号，将v和列车受力情况用于式(4)可得列车功率P。在给定时刻，牵引网络的潮流可在5个供电区间中独立地分析，对于其中一个供电区间的分析过程如下：
选出该供电区间内的列车子集，记x和P为其位置和功率。将x代入到式(5)中可得该供电区间的导纳矩阵G。最终设置终止条件kmax和(ΔU)tol，即可应用潮流分析方法来计算网络的潮流。
仿真比较：
比较了几种高效的数值方法：经典Newton-Raphson方法(NR)[2]，线性逼近方法(LA)[5]，连续逼近方法(SA)[4]以及本文提出的压缩Newton-Raphson方法(CNR)。仿真测试在64位Windows10操作系统的台式机上运行，采用频率为4.10GHz的Inter(R)core(TM)i5-9300H的处理器和16GB的内存。为了公平地比较不同方法的性能，本发明为各方法设置相同的终止条件：kmax＝100和(ΔU)tol＝10-10。为验证方法的有效性，使用CNR方法求解得到的第一辆车的电压曲线绘制于图7，可见该方法收敛到的解在一个适当的区间内(700-850V)。
本发明使用Q-convergence定义分析各方法的收敛速度的阶数，设序列(Uk)收敛至U*，则若存在阶数q≥1使得式：
Figure BDA0003673878870000131
成立，其中M为正数，则称序列以q阶速度收敛。在实验中可通过序列的||Uk-U*||曲线比较其收敛速度，如图8。可见NR和CNR方法为二次收敛，然而LA和SA方法为线性收敛，且LA方法收敛得比SA方法更快。
各方法的平均迭代次数列于表5。当s从5增加至6时，LA和SA方法的平均迭代次数增长幅度较大，而NR和CNR方法则只有轻微的增长。
表5平均迭代次数
Figure BDA0003673878870000132
各方法的平均每轮迭代的运行时间列于表6。
表6单轮迭代运行时间
Figure BDA0003673878870000133
结果表明SA方法是单轮迭代计算效率最高的方法，因为它不需要计算Jacobian矩阵，也不含矩阵求逆运算。然而如前所述，SA方法仍含有矩阵运算，因此其计算复杂度为O(s2)。当s从5增加到6时，SA方法的计算时间增加了近一半。相比之下，不含任何矩阵运算的CNR方法的计算复杂度仅为O(s)，s的增加只导致其计算时间有轻微的增长。然而，由于自动微分操作导致CNR方法的计算复杂度的常系数很大，所以当s较小时，其线性复杂度的优势不能充分体现。

Claims (2)

1.一种DC牵引网络的潮流分析方法，其特征在于，包括以下步骤：
步骤1：构建DC牵引网络模型：
采用第三轨供电的DC牵引网络由变电站、悬链线、回程轨道和列车组成，其中列车运行在上行或下行方向的轨道；两个变电站sL和sR分别位于网络的前后端，其位置分别为xL和xR，在两个变电站之间，有n辆列车运行在上行方向轨道，其位置分别为
Figure FDA0003673878860000011
有m辆列车运行在下行方向轨道，其位置分别为
Figure FDA0003673878860000012
绘制DC牵引网络的等效电路，其中每辆列车等效为一个恒功率源，每个变电站等效为一个带有串联电阻Rs的恒定直流电压源，其电阻值和电压值分别为变电站的内阻和空载电压；等效电路中，
Figure FDA0003673878860000013
为对应列车的功率，由下式计算；
P(v)＝(Ft/ηt-αηbFb)v (4)
其中，0＜ηt,ηb＜1分别为列车在牵引模式和制动模式下的机电效率的有效值，0＜α＜1为列车制动能量的利用率；Ft和Fb分别为列车所受的牵引力和制动力，v为列车的速度；
Figure FDA0003673878860000014
为对应列车所处悬链线的电压，UL和UR为两个变电站的输出电压；
Figure FDA0003673878860000015
为对应的列车与列车之间或列车与变电站之间的悬链线和回程轨道的等效电导值，由下式计算：
Figure FDA0003673878860000016
其中，ρc和ρr分别为悬链线和回程轨道的电阻率；
步骤2：推导潮流方程：
DC牵引网络的等效电路的结点电压方程如下式：
Figure FDA0003673878860000017
式中，
Figure FDA0003673878860000018
等号右侧表示每个结点的输入电流，由Tellegen第一定理推导为：
Figure FDA0003673878860000021
若令
Figure FDA0003673878860000022
以及
Figure FDA0003673878860000023
则将式(6)表达为矩阵形式GU＝I，其中G表示导纳矩阵，由于各结点的位置已知，G中的每一项均可由式(5)得到；
步骤3：求解潮流方程：
将式(6)和(7)结合并重构得到：
Figure FDA0003673878860000024
其中，
Figure FDA0003673878860000025
当
Figure FDA0003673878860000026
和
Figure FDA0003673878860000027
给定时，U可被递推地计算得到；因此，将
Figure FDA0003673878860000028
和
Figure FDA0003673878860000029
视为决策变量，即定义
Figure FDA00036738788600000210
分析可知，U满足潮流方程(6)和(7)当且仅当Uc满足
Figure FDA00036738788600000211
即Uc为函数
Figure FDA00036738788600000212
的根；fc的Jacobian矩阵Jc由下式计算：
Figure FDA00036738788600000213
定义在计算图中有从结点wi指向wj的箭头，则wj称为wi的子结点，同时wi称为wj的父结点；
Figure FDA00036738788600000214
为wi的子结点集，对应地，wj的父结点集表示为
Figure FDA00036738788600000215
在每一次迭代中具有前向传播和反向传播两个步骤；
在前向传播中，使用式(8)计算fc(Uc)，同时计算每一个中间变量结点wj对其父结点集中每个结点的偏导，即
Figure FDA00036738788600000216
由下式给出：
Figure FDA0003673878860000031
在反向传播中，根据计算图以及微分的链式法则，fc对每个结点的偏导数可由式(11)计算，其计算顺序与前向传播过程中的相反；
Figure FDA0003673878860000032
由fc的定义，可得式(11)的初始条件为：
Figure FDA0003673878860000033
接着，根据式(9)即得到Jc，根据Newton-Raphson方法的计算框架，若给定初始电压向量
Figure FDA0003673878860000034
则第k次迭代的更新公式为：
Figure FDA0003673878860000035
最终，根据是否满足终止条件来决定将结束还是进入下一次迭代。
2.根据权利要求1所述的一种DC牵引网络的潮流分析方法，其特征在于，所述迭代终止条件为：
k≥kmax or||ΔU||≤(ΔU)tol
其中，kmax为最大迭代次数，(ΔU)tol为最小电压变化步长容量，||·||为向量的L2范数；
Figure FDA0003673878860000036
为电压变化步长。
CN202210614934.9A 2022-06-01 2022-06-01 一种dc牵引网络的潮流分析方法 Active CN114912068B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210614934.9A CN114912068B (zh) 2022-06-01 2022-06-01 一种dc牵引网络的潮流分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210614934.9A CN114912068B (zh) 2022-06-01 2022-06-01 一种dc牵引网络的潮流分析方法

Publications (2)

Publication Number Publication Date
CN114912068A CN114912068A (zh) 2022-08-16
CN114912068B true CN114912068B (zh) 2023-03-10

Family

ID=82771643

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210614934.9A Active CN114912068B (zh) 2022-06-01 2022-06-01 一种dc牵引网络的潮流分析方法

Country Status (1)

Country Link
CN (1) CN114912068B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115498681B (zh) * 2022-11-08 2023-03-24 清华大学 柔性直流牵引供电系统的近似最优潮流计算方法及装置
CN116388186B (zh) * 2023-06-06 2023-08-25 清华大学 交流牵引供电系统的潮流计算方法、装置、设备及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106253289A (zh) * 2016-08-12 2016-12-21 成都运达科技股份有限公司 一种车网耦合的地铁供电系统电力潮流计算方法
CN106655196A (zh) * 2017-02-15 2017-05-10 河海大学 一种轨道交通的供电网潮流计算方法
CN108183487A (zh) * 2018-01-23 2018-06-19 中国石油大学(华东) 一种基于线性数学模型的配电网潮流快速分析方法
CN109995037A (zh) * 2017-12-29 2019-07-09 湖南工业大学 计及交直流耦合的牵引供电系统潮流分析方法、系统及存储介质
CN110053521A (zh) * 2019-03-08 2019-07-26 北京交通大学 城市轨道交通牵引供电系统及车-网配合参数优化方法
CN110932282A (zh) * 2019-12-25 2020-03-27 福州大学 一种增广直角坐标下的基于vsc内部修正方程矩阵与交替迭代法的潮流计算方法
CN112350326A (zh) * 2020-10-23 2021-02-09 株洲中车时代电气股份有限公司 轨道交通牵引供电系统及其控制方法、系统及相关组件

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106253289A (zh) * 2016-08-12 2016-12-21 成都运达科技股份有限公司 一种车网耦合的地铁供电系统电力潮流计算方法
CN106655196A (zh) * 2017-02-15 2017-05-10 河海大学 一种轨道交通的供电网潮流计算方法
CN109995037A (zh) * 2017-12-29 2019-07-09 湖南工业大学 计及交直流耦合的牵引供电系统潮流分析方法、系统及存储介质
CN108183487A (zh) * 2018-01-23 2018-06-19 中国石油大学(华东) 一种基于线性数学模型的配电网潮流快速分析方法
CN110053521A (zh) * 2019-03-08 2019-07-26 北京交通大学 城市轨道交通牵引供电系统及车-网配合参数优化方法
CN110932282A (zh) * 2019-12-25 2020-03-27 福州大学 一种增广直角坐标下的基于vsc内部修正方程矩阵与交替迭代法的潮流计算方法
CN112350326A (zh) * 2020-10-23 2021-02-09 株洲中车时代电气股份有限公司 轨道交通牵引供电系统及其控制方法、系统及相关组件

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Multi-train modeling and simulation integrated with traction power supply solver using simplified Newton-Raphson method;Thanatchai Kulworawanichpong;《Journal of Modern Transportation》(第24期);241-251 *

Also Published As

Publication number Publication date
CN114912068A (zh) 2022-08-16

Similar Documents

Publication Publication Date Title
CN114912068B (zh) 一种dc牵引网络的潮流分析方法
Domínguez et al. Energy savings in metropolitan railway substations through regenerative energy recovery and optimal design of ATO speed profiles
Calderaro et al. Optimal siting and sizing of stationary supercapacitors in a metro network using PSO
CN104332996A (zh) 一种评估电力系统可靠性的方法
CN101533061A (zh) 基于稀疏pmu配置的大型输电网络故障定位方法
CN107689459A (zh) 一种有轨电车用燃料电池阵列系统的效率优化控制方法
He et al. Optimal control of metro energy conservation based on regenerative braking: A complex model study of trajectory and overlap time
CN105303021A (zh) 城市轨道交通直流牵引供电系统自适应实时动态数学建模方法
Mahyiddin et al. Fuzzy logic energy management system of series hybrid electric vehicle
Calderaro et al. An algorithm to optimize speed profiles of the metro vehicles for minimizing energy consumption
Chen et al. Cooperative eco-driving of multi-train under dc traction network
Diab et al. A complete DC trolleybus grid model with bilateral connections, feeder cables, and bus auxiliaries
Ahmadi et al. Improving energy‐efficient train operation in urban railways: employing the variation of regenerative energy recovery rate
Zhang et al. Research on multi‐train energy saving optimization based on cooperative multi‐objective particle swarm optimization algorithm
Chuang et al. Design of optimal coasting speed for saving social cost in mass rapid transit systems
CN110309550B (zh) 一种基于势能场与网络效率的高速列车系统可靠性分析方法
CN103544373A (zh) 一种智能配电网清洁性评估方法
Guo et al. Adaptive fuzzy sliding mode control for high‐speed train using multi‐body dynamics model
Alnuman et al. Electrical modelling of a metro system
Doan et al. The design of an optimal running curve for train operation based on a novel parameterization method aiming to minimize the total energy consumption
Battistelli et al. Generalized approach to design supercapacitor-based storage devices integrated into urban mass transit systems
Prema et al. Application of hybrid neuro-wavelet models for effective prediction of wind speed
Sato et al. A method of generating energy-efficient train timetable including charging strategy for catenary-free railways with battery trains
CISMARU Energy efficient train operation using simulated annealing algorithm and simulink model
Nomura et al. Scheduling method for minimum energy consumption considering constraints of time intervals between local and express trains

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant