CN110336287B - 一种基于雅可比元素提取的配电系统三相潮流计算方法 - Google Patents
一种基于雅可比元素提取的配电系统三相潮流计算方法 Download PDFInfo
- Publication number
- CN110336287B CN110336287B CN201910695403.5A CN201910695403A CN110336287B CN 110336287 B CN110336287 B CN 110336287B CN 201910695403 A CN201910695403 A CN 201910695403A CN 110336287 B CN110336287 B CN 110336287B
- Authority
- CN
- China
- Prior art keywords
- node
- phase
- matrix
- elements
- formula
- 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
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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/04—Circuit arrangements for ac mains or ac distribution networks for connecting networks of the same frequency but supplied from different sources
- H02J3/06—Controlling transfer of power between connected networks; Controlling sharing of load between connected 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Power Engineering (AREA)
- Supply And Distribution Of Alternating Current (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于雅可比元素提取的配电系统三相潮流计算方法,利用Matlab具有丰富的函数并且擅长矩阵运算的特点来形成雅可比矩阵有效行标号数组和有效列标号数组。采用Matlab的sum函数和abs函数计算雅可比矩阵各行元素或各列元素绝对值之和,如果雅可比矩阵的某行元素或某列元素绝对值之和为0,则该行内或该列内元素全部为0;利用Matlab的find函数查找所计算各行元素或各列元素绝对值之和不等于0的行或列,进而形成雅可比矩阵有效行标号数组和有效列标号数组。由于Matlab的内置函数和矩阵运算速度很快,本发明大大缩短了形成雅可比矩阵有效行数组和有效列数组的时间,提高了计算速度,也简化了编程。
Description
技术领域
本发明涉及一种配电系统三相潮流计算方法,特别是一种基于雅可比矩阵元素提取的配电系统三相潮流计算方法。
背景技术
电力系统是由生产、变换、输送、分配和使用电能的设备(元件)组成的十分庞大而复杂的统一体,分为输电系统和配电系统两部分。输电系统主要起变换和输送电能的作用,电压等级较高;配电系统主要起变换和分配电能的作用,电压等级较低。
电力系统潮流计算是研究电力系统稳态运行的一项基本计算,它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础,如安全分析、暂态稳定分析等都要用到潮流计算。
配电系统是由配电线路、配电变压器、配电开关、配电电容器和配电负荷等组成的直接向终端用户分配电能的一个网络系统。对配电系统进行电力系统分析首先需要进行潮流计算,得到它的运行状态。与高压输电系统的三相对称运行方式不同,配电系统的负荷和网络都可能不对称,配电系统进行潮流计算时,应考虑三相不对称的特点,建立各元件的三相模型,进行三相潮流计算。
三相潮流计算的数学模型和算法设计都非常复杂,编程和调试工作量都大,三相潮流计算研究迫切需要易于编程和调试的软件开发环境。Matlab平台可以方便地处理各种矩阵和复数运算,有大量常见且实用的函数,给编程带来很大便利,同时还具有非常强大的可视化图形处理和交互式功能,为科学研究以及工程应用提供了一种高效的编程工具,目前已经成为许多科学领域的基本工具和首选平台,在各种科学和工程计算领域得到了广泛的应用。因此Matlab平台很适合三相潮流计算算法的开发和研究。但Matlab为解释型编程语言,通常对C语言等编译型编程语言非常有效的实现方法,用Matlab实现时计算速度则很慢,因此设计潮流计算方法时要充分考虑Matlab的特点。
根据电力系统节点的特点,潮流计算把电力系统节点分成3类:节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点;节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点;节点电压幅值和电压相角已知、节点有功功率和无功功率未知的节点称为平衡节点。
潮流计算的功率方程或电流方程是非线性方程,必须采用迭代方式求解。迭代方法求解潮流计算问题时,需要设定各节点电压的初值。电压初始化一般采用平启动,单相潮流计算采用平启动设定电压初值时,PV节点和平衡节点的节点电压幅值取给定值,PQ节点的节点电压幅值取1.0;所有节点电压的相角都取0.0。三相潮流计算采用平启动设定电压初值时,PV节点和平衡节点的三相节点电压幅值都取给定值,PQ节点的三相节点电压幅值都取1.0;所有节点A相节点电压的相角都取0°,B相节点电压的相角都取-120°,C相节点电压的相角都取120°。这里单位除相角外都采用标幺值。
设配电系统的节点数为n,其中PV节点数为m。为叙述方便,设PV节点的编号为1~m。
三相潮流计算的节点电压方程为:
I=YU (1)
式中,I为三相注入电流相量列向量,Y为3n×3n的导纳矩阵,U为三相节点电压相量列向量。
导纳矩阵由输电线路、变压器支路和无功补偿元件形成,为3n×3n阶矩阵。配电系统的输电线路一般为三相,也有的为两相或单相,因此,导纳矩阵中有些行或列的元素全部为0。
各相注入复功率为:
式中,为节点i的A相复功率,为节点i的B相复功率,为节点i的C相复功率,为节点i的A相电压相量,为节点i的B相电压相量,为节点i的C相电压相量,为节点i的A相注入电流相量,为节点i的B相注入电流相量,为节点i的C相注入电流相量,上标*表示共轭,i=1、2、…、n。
由式(1),得:
式中,导纳矩阵元素为节点i的A相和节点k的A相之间的互导纳,当k=i时,为节点i的A相的自导纳;为节点i的A相和节点k的B相之间的互导纳;为节点i的A相和节点k的C相之间的互导纳;为节点i的B相和节点k的B相之间的互导纳,当k=i时,为节点i的B相的自导纳;为节点i的B相和节点k的A相之间的互导纳;为节点i的B相和节点k的C相之间的互导纳;为节点i的C相和节点k的C相之间的互导纳,当k=i时,为节点i的C相的自导纳;为节点i的C相和节点k的A相之间的互导纳;为节点i的C相和节点k的B相之间的互导纳;k=1、2、…、n。
由式(2),得:
将式(4)减去式(3)得到:
为了提高配电系统运行的可靠性,中性点一般不接地。当中性点不接地时,中性点注入电流为0,即:
将式(4)代入式(6),得:
对于PQ节点,已知其每一相的有功功率和无功功率或者三相总有功功率和三相总无功功率,求其电压相量。当中性点不接地时,由中性点注入电流等于0可得,3个单相有功功率和3个单相无功功率并不是完全可控的,只能控制其中4个量,另外两个是状态量。根据中性点注入电流为0的条件,多出了中性点注入电流实部和虚部为0的两个方程,共有8个方程和6个状态量。8个方程分别是A、B、C相和中性点N的电流实部偏差和虚部偏差方程,6个状态量为A、B、C相电压实部和虚部,所以需要增加2个状态量,使方程个数和状态量个数一致。控制变量选取总有功功率Pi∑和3个单相无功功率,其值不变,为给定值;A相有功功率Pi A和B相有功功率Pi B作为状态变量,其值待求,为计算值,C相有功功率由Pi∑、Pi A和Pi B求得,也是计算值。当中性点接地时,去掉中性点N的电流实部偏差和虚部偏差方程,以及将A相有功功率Pi A和B相有功功率Pi B作为控制变量,为已知量。为了列写方程方便,中性点接地的节点也列写中性点N的电流实部偏差和虚部偏差方程,且中性点N的电流实部偏差和虚部偏差都置为0,对应方程各变量的系数也都置为0,相应节点的Pi A和Pi B的修正量也为0。
PV节点有足够的可调节的无功功率容量,使得它可以保证电压幅值的恒定,一般将有无功功率储备的发电厂和具有调节无功功率能力的变电所选作PV节点。PV节点电压幅值是给定的,分别为:这样与PQ节点相比,PV节点需要增加3个关于电压幅值偏差的方程;同时需要将控制变量 变为状态变量,其值为计算值。
电压幅值偏差公式如下:
可见,PV节点比PQ节点多了3个状态变量修正量和3个偏差量。为了设计方便,可以把PV节点的状态变量修正量和偏差量分成两部分,一部分是与PQ节点相同的部分,另一部分则是多出的3个状态变量修正量和3个偏差量。
对PQ节点和PV节点都适用的状态变量修正量为:
式中,和Δfi A分别为的实部修正量和虚部修正量,和Δfi B分别为的实部修正量和虚部修正量,和Δfi C分别为的实部修正量和虚部修正量,ΔPi A和ΔPi B分别为Pi A和Pi B的修正量,上标T表示转置。
对PQ节点和PV节点都适用的偏差量为:
PV节点附加的状态变量修正量为:
PV节点附加的偏差量为:
由式(9)和式(11),得到状态变量修正量ΔX为:
由式(10)和式(12),得到偏差量列向量ΔW为:
平衡节点不参与迭代计算,不需要计算电流偏差或电压偏差,也不需要计算电压修正量和功率修正量,迭代完成后再计算各相功率。
潮流计算的修正方程为:
ΔW=-JΔX (15)
式中,J为雅可比矩阵,ΔW为偏差量列向量,ΔX为状态变量修正量列向量。
将雅可比矩阵写成如下的分块矩阵:
式中,J0为8n×8n阶基本雅可比子矩阵,D为8m×3m阶子矩阵,F为3m×8m阶子矩阵,O1为8(n-m)×3m阶零矩阵,O2为3m×8(n-m)阶零矩阵,O3为3m×3m阶零矩阵。
基本雅可比子矩阵J0为:
式中,Jik、JDii为分块子矩阵,diag表示对角矩阵。
子矩阵Jik为:
式中,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部;分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部;分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部。
子矩阵JDii为:
式中元素分别表示如下:
式(19)中元素分别表示如下:
如果节点i的中性点接地,则式(19)的第4行和第8行、第4列和第8列元素清零。
D为分块对角阵:
D=diag(D1,D2,…,Dm) (24)
如果节点i的中性点接地,则式(25)的第4行和第8行元素清零。
F为分块对角阵:
F=diag(F1,F2,…,Fm) (26)
如图1-3所示,现有基于电流偏差型直角坐标牛顿法三相潮流计算方法,主要包括以下步骤:
A、输入原始数据和初始化电压。
B、形成三相导纳矩阵。
根据输入的输电线路、变压器和无功补偿元件参数形成三相导纳矩阵。
C、令t=0。
D、形成偏差量列向量ΔW并计算最大偏差量ΔWmax。
根据式(5)和式(7)计算各节点的注入电流相量偏差,根据式(8)计算PV节点的电压相量偏差,形成偏差量列向量ΔW,并计算最大偏差量ΔWmax,即绝对值最大的偏差量。
E、判断|ΔWmax|是否满足收敛精度ε,如果满足,转至步骤L。
F、形成基本雅可比子矩阵J0。
计算式(17)所示的基本雅可比子矩阵J0各分块子矩阵,其中J0各分块子矩阵Jik按式(18)计算;J0各分块子矩阵JDii按式(19)-(23)计算。形成基本雅可比子矩阵J0后,平衡节点对应的基本雅可比子矩阵J0行和列元素清零。
G、形成雅可比矩阵J。
在基本雅可比子矩阵J0的基础上追加与PV节点有关的附加雅可比矩阵元素形成完整的雅可比矩阵J,计算过程如下:
按式(25)计算分块对角矩阵D的各分块矩阵,并形成如式(24)所示的对角矩阵D;按式(27)计算分块对角矩阵F的各分块矩阵,并形成如式(26)所示的对角矩阵F。由基本雅可比子矩阵J0、对角矩阵D、对角矩阵F以及零矩阵O1、O2、O3形成完整的雅可比矩阵J。
H、判断t=0是否成立,如果不成立,转至步骤J。
I、形成雅可比矩阵有效行标号数组JR和有效列标号数组JC。
雅可比矩阵J中有许多整行元素全部为0的行或整列元素全部为0的列,为无效的行或列,进行计算时,需要去掉这些行和列,只取有效的行和列。具体有两种方法,第一种方法是根据Matlab按列存储矩阵数据的特点,按列取数据较快,按行取数据较慢,因此可以通过矩阵转置将判断矩阵一行中元素是否全部为0的运算转化为判断矩阵一列中元素是否全部为0的运算;第二种方法是用Matlab的find函数查找并记录雅可比矩阵一行(列)内不为0元素的列(行)号,再用length函数计算其个数,进而判断该行(列)元素是否全部为0。
选择形成雅可比矩阵有效行标号数组和有效列标号数组的方法,如果选择第一种方法,则转步骤I1;如果选择第二种方法,则转步骤I21;
I1、令V=JT;
I2、令l=1,s=0,u=0;
I3、令q=1;
I4、判断Vql=0是否成立,如果成立,转步骤I6;
I5、令s=s+1,JRs=l,转步骤I8;
I6、令q=q+1;
I7、判断q大于雅可比矩阵的列数N是否成立,如果不成立,转步骤I4;
I8、令q=1;
I9、判断Jql=0是否成立,如果成立,转步骤I11;
I10、令u=u+1,JCu=l,转步骤I13;
I11、令q=q+1;
I12、判断q大于雅可比矩阵的行数N是否成立,如果不成立,转步骤I9;
I13、令l=l+1;
I14、判断l大于雅可比矩阵的行数N是否成立,如果成立,转步骤J;否则转步骤I3;
I21、令V=JT;
I22、令l=1,s=0,u=0;
I23、令A=find(V(:,l));
I24、令q=length(A);
I25、判断q=0是否成立,如果成立,转步骤I27;
I26、令s=s+1,JRs=l;
I27、令A=find(J(:,l));
I28、令q=length(A);
I29、判断q=0是否成立,如果成立,转步骤I31;
I30、令u=u+1,JCu=l;
I31、令l=l+1;
I32、判断l大于雅可比矩阵的行数N是否成立,如果成立,转步骤J;否则转步骤I23;
J、解修正方程并修正状态变量;
考虑去掉雅可比矩阵中元素全部为0的行和列后的修正方程为
ΔWJR=-JJR,JCΔXJC (28)
式中,ΔWJR为按数组JR记录的标号提取偏差量列向量ΔW的元素形成的新列向量;JJR,JC为按数组JR记录的行标号和JC记录的列标号提取雅可比矩阵J的元素形成的新矩阵;ΔXJC为按数组JC记录的标号提取状态变量修正量ΔX的元素形成的新列向量。
利用Matlab的除法运算求解式(28)所示的修正方程,得到状态变量修正量ΔXJC如下:
ΔXJC=-JJR,JC\ΔWJR (29)
式中,\为Matlab的除法运算符号,其左侧是分母,右侧是分子。
按式(30)修正状态变量如下:
K、令t=t+1,转至步骤D。
L、输出节点及支路数据,结束。
由于配电系统的规模较大、雅可比矩阵维数很高,在Matlab编程时采用传统第一种方法或第二种方法形成雅可比矩阵有效行标号数组和有效列标号数组时,计算时间很长,难以满足三相潮流计算对速度的要求。
发明内容
为解决现有技术存在的上述问题,本发明要提出一种能大大缩短计算时间,提高计算速度的基于雅可比元素提取的配电系统三相潮流计算方法。
本发明的基本思路是充分利用Matlab具有丰富的函数并且擅长矩阵运算的特点来形成雅可比矩阵有效行标号数组和有效列标号数组。采用Matlab的sum函数和abs函数计算雅可比矩阵各行元素或各列元素绝对值之和,如果雅可比矩阵的某行元素或某列元素绝对值之和为0,则该行内或该列内元素全部为0;利用Matlab的find函数查找所计算各行元素或各列元素绝对值之和不等于0的行或列,进而形成雅可比矩阵有效行标号数组和有效列标号数组。sum函数和find函数都是针对整个矩阵的运算,不是具体某个矩阵元素的运算,运算速度非常快。
为了实现上述目的,本发明的技术方案如下:一种基于雅可比元素提取的配电系统三相潮流计算方法,包括以下步骤:
A、输入原始数据和初始化电压。
B、形成三相导纳矩阵。
根据输入的输电线路、变压器和无功补偿元件参数形成三相导纳矩阵。
C、令t=0。
D、形成偏差量列向量ΔW并计算最大偏差量ΔWmax。
根据式(5)和式(7)计算各节点的注入电流相量偏差,根据式(8)计算PV节点的电压相量偏差,形成偏差量列向量ΔW,并计算最大偏差量ΔWmax,即绝对值最大的偏差量。
E、判断|ΔWmax|是否满足收敛精度ε,如果满足,转至步骤L。
F、形成基本雅可比子矩阵J0。
计算式(17)所示的基本雅可比子矩阵J0各分块子矩阵,其中J0各分块子矩阵Jik按式(18)计算;J0各分块子矩阵JDii按式(19)-(23)计算。形成基本雅可比子矩阵J0后,平衡节点对应的基本雅可比子矩阵J0行和列元素清零。
G、形成雅可比矩阵J。
在基本雅可比子矩阵J0的基础上追加与PV节点有关的附加雅可比矩阵元素形成完整的雅可比矩阵J,计算过程如下:
按式(25)计算分块对角矩阵D的各分块矩阵,并形成如式(24)所示的对角矩阵D;按式(27)计算分块对角矩阵F的各分块矩阵,并形成如式(26)所示的对角矩阵F。由基本雅可比子矩阵J0、对角矩阵D、对角矩阵F以及零矩阵O1、O2、O3形成完整的雅可比矩阵J。
H、判断t=0是否成立,如果不成立,转至步骤J。
I、形成雅可比矩阵有效行标号数组JR和有效列标号数组JC。
雅可比矩阵J中有许多整行元素全部为0的行或整列元素全部为0的列,为无效的行或列,进行计算时,需要去掉这些行和列,只取有效的行和列。用Matlab的sum函数和abs函数计算雅可比矩阵各行元素或各列元素绝对值之和;如果各行元素或各列元素绝对值之和为0,则雅可比矩阵对应行或对应列的元素全部为0,从而找到了矩阵J中元素都为0的行或元素都为0的列。形成雅可比矩阵有效行标号数组和有效列标号数组的具体步骤如下:
I1、令矩阵A为由雅可比矩阵J元素的绝对值组成的矩阵,即:
A=abs(J) (31)
式中,abs为Matlab中对矩阵元素取绝对值的函数。
I2、利用Matlab的函数sum求矩阵A各行或各列的和,得到数组R和C:
R=sum(A,2) (32)
C=sum(A,1) (33)
式中,sum为Matlab中矩阵元素求和函数,sum的第2个参数为1时表示对各列元素求和,第2个参数为2时表示对各行元素求和。
I3、分别用式(34)和式(35)提取数组R和C中不为0元素的标号:
JR=find(R) (34)
JC=find(C) (35)
式中,find为Matlab中查找不为0的数组元素的标号。
J、解修正方程并修正状态变量。
考虑去掉雅可比矩阵中元素全部为0的行和列后的修正方程为:
ΔWJR=-JJR,JCΔXJC (36)
式中,ΔWJR为按数组JR记录的标号提取偏差量列向量ΔW的元素形成的新列向量;JJR,JC为按数组JR记录的行标号和JC记录的列标号提取雅可比矩阵J的元素形成的新矩阵;ΔXJC为按数组JC记录的标号提取状态变量修正量ΔX的元素形成的新列向量。
按式(37)解式(36)所示的修正方程,得到状态变量修正量ΔXJC如下:
ΔXJC=-JJR,JC\ΔWJR (37)
按式(38)修正状态变量如下:
K、令t=t+1,转至步骤D。
L、输出节点及支路数据,结束。
与现有技术相比,本发明具有以下有益效果:
本发明采用Matlab的sum函数和abs函数来计算矩阵各行或各列的元素绝对值的和,用find函数查找所求元素绝对值和不为0的行号或列号,进而形成雅可比矩阵J有效行数组JR和有效列数组JC。由于Matlab的内置函数和矩阵运算速度很快,与传统采用循环方法或逻辑函数方法相比,本发明方法大大缩短了形成雅可比矩阵有效行数组和有效列数组的时间,提高了计算速度,同时也简化了编程。
附图说明
本发明共有附图4张。其中:
图1是现有配电系统三相潮流计算的流程图。
图2是现有形成雅可比矩阵有效行和有效列数组的方法1的流程图。
图3是现有形成雅可比矩阵有效行和有效列数组的方法2的流程图。
图4是本发明方法配电系统三相潮流计算的流程图。
具体实施方式
下面结合附图对本发明进行进一步地说明,按照图4所示配电系统三相潮流计算流程对4个大型实际配电系统算例进行了计算。
算例1:有269个节点,264条线路,7个变压器,220个负荷节点,2个电源节点,0个补偿电容器节点。
算例2:有665个节点,653条线路,13个变压器,605个负荷节点,2个电源节点,1个补偿电容器节点。
算例3:有1447个节点,1437条线路,12个变压器,548个负荷节点,2个电源节点,5个补偿电容器节点。
算例4:有2661个节点,2647条线路,16个变压器,1498个负荷节点,2个电源节点,5个补偿电容器节点。
采用本发明和几种对比方法对这4个实际系统算例进行了计算,各种潮流计算均采用了Matlab的稀疏技术,计算时收敛精度为0.00001。几种对比三相潮流计算算法分别为:
方法1:采用传统方法1形成雅可比矩阵有效行数组和有效列数组的三相潮流计算算法。
方法2:采用传统方法2形成雅可比矩阵有效行数组和有效列数组的三相潮流计算算法。
方法3:本发明方法。
3种对比方法计算各算例时有效数组形成时间和潮流计算时间见表1~表4,潮流计算的计算时间不包括数据读入和输出及支路功率计算的时间。
表1几种三相潮流算法计算算例1时的时间比较
潮流计算方法 | 有效数组形成时间(s) | 潮流计算时间(s) | 时间占比 |
方法1 | 13.75322 | 14.76090 | 93.17% |
方法2 | 0.10973 | 1.11541 | 9.84% |
方法3 | 0.00275 | 1.00768 | 0.27% |
从表1可见,对于具有269个节点的系统,传统方法1形成雅可比矩阵有效行数组和有效列数组的时间较长,占整个潮流计算的计算时间的93.17%;传统方法2形成雅可比矩阵有效行数组和有效列数组的时间虽然大大减少,但占整个潮流计算的计算时间也高达9.84%,对形成雅可比矩阵有效行数组和有效列数组这一简单功能的计算时间也显得太长了。本发明形成雅可比矩阵有效行数组和有效列数组的时间非常短,仅占整个潮流计算的计算时间的0.27%。
表2几种三相潮流算法计算算例2时的时间比较
从表2可见,对于具有665个节点的系统,传统方法1形成雅可比矩阵有效行数组和有效列数组的时间很长,竟高达127s,占整个潮流计算的计算时间的97.55%,无法接受;传统方法2形成雅可比矩阵有效行数组和有效列数组的时间虽然大大减少,但占整个潮流计算的计算时间也高达12.82%,对形成雅可比矩阵有效行数组和有效列数组这一简单功能的计算时间也显得太长了。本发明形成雅可比矩阵有效行数组和有效列数组的时间非常短,仅占整个潮流计算的计算时间的0.09%。
表3几种三相潮流算法计算算例3时的时间比较
潮流计算方法 | 有效数组形成时间(s) | 潮流计算时间(s) | 时间占比 |
方法2 | 1.823120 | 14.86029 | 12.27% |
方法3 | 0.00433 | 13.04117 | 0.03% |
对于具有1447个节点的系统,传统方法1形成雅可比矩阵有效行数组和有效列数组的时间非常长,无法接受,不得不中断运行,故没有得到计算时间;从表3可见,传统方法2形成雅可比矩阵有效行数组和有效列数组的时间虽然较少,但占整个潮流计算的计算时间也高达12.27%,对形成雅可比矩阵有效行数组和有效列数组这一简单功能的计算时间也显得太长了。本发明形成雅可比矩阵有效行数组和有效列数组的时间非常短,仅占整个潮流计算的计算时间的0.03%。
表4几种三相潮流算法计算算例4时的时间比较
潮流计算方法 | 有效数组形成时间(s) | 潮流计算时间(s) | 时间占比 |
方法2 | 5.86933 | 51.42213 | 11.41% |
方法3 | 0.00735 | 45.55681 | 0.016% |
对于具有2661个节点的系统,传统方法1形成雅可比矩阵有效行数组和有效列数组的时间非常长,无法接受,不得不中断运行,故没有得到计算时间;从表4可见,传统方法2形成雅可比矩阵有效行数组和有效列数组的时间较长,占整个潮流计算的计算时间高达11.41%,对形成雅可比矩阵有效行数组和有效列数组这一简单功能的计算时间也显得太长了。本发明形成雅可比矩阵有效行数组和有效列数组的时间仍然非常短,仅占整个潮流计算的计算时间的0.016%。本发明不仅计算时间短,程序也进一步简化。
本发明可以在任何版本的MATLAB编程语言实现,但建议使用较新版本的MATLAB语言。
本发明不局限于本实施例,任何在本发明披露的技术范围内的等同构思或者改变,均列为本发明的保护范围。
Claims (1)
1.一种基于雅可比元素提取的配电系统三相潮流计算方法,包括以下步骤:
A、输入原始数据和初始化电压;
根据电力系统节点的特点,潮流计算把电力系统节点分成3类:节点有功功率和无功功率已知,节点电压幅值和电压相角未知的节点称为PQ节点;节点有功功率和电压幅值已知,节点无功功率和电压相角未知的节点称为PV节点;节点电压幅值和电压相角已知,节点有功功率和无功功率未知的节点称为平衡节点;
B、形成三相导纳矩阵;
C、设置迭代计数t=0;
D、形成偏差量列向量ΔW并计算最大偏差量ΔWmax;
PQ节点和PV节点的电流相量偏差公式为:
式中,为节点i的A相注入电流相量偏差,为节点i的B相注入电流相量偏差,为节点i的C相注入电流相量偏差;为节点i的中性点注入电流相量偏差;为节点i的A相复功率,为节点i的B相复功率,为节点i的C相复功率;为节点i的A相电压相量,为节点i的B相电压相量,为节点i的C相电压相量;导纳矩阵元素为节点i的A相和节点k的A相之间的互导纳,当k=i时,为节点i的A相的自导纳;为节点i的A相和节点k的B相之间的互导纳;为节点i的A相和节点k的C相之间的互导纳;为节点i的B相和节点k的B相之间的互导纳,当k=i时,为节点i的B相的自导纳;为节点i的B相和节点k的A相之间的互导纳;为节点i的B相和节点k的C相之间的互导纳;为节点i的C相和节点k的C相之间的互导纳,当k=i时,为节点i的C相的自导纳;为节点i的C相和节点k的A相之间的互导纳;为节点i的C相和节点k的B相之间的互导纳;上标*表示共轭;k=1、2、…、n;
PV节点的电压幅值偏差公式如下:
式中,为节点i的A相电压幅值偏差,为节点i的B相电压幅值偏差,为节点i的C相电压幅值偏差;为节点i的A相给定电压幅值,为节点i的B相给定电压幅值,为节点i的C相给定电压幅值;fi A分别为的实部和虚部,fi B分别为的实部和虚部,fi C分别为的实部和虚部;
偏差量列向量ΔW为:
式中,n为节点数;m为PV节点数;设PV节点的编号为1~m;
式(4)中ΔWi P为:
式(4)中ΔWi PV0为:
在偏差量列向量ΔW中查找绝对值最大的值,得到最大偏差量ΔWmax;
平衡节点不参与迭代计算,不需要计算电流相量偏差或电压幅值偏差;
E、判断|ΔWmax|是否满足收敛精度ε,如果满足,转至步骤L;否则,执行步骤F;
F、形成基本雅可比子矩阵J0;
基本雅可比子矩阵J0为:
式中,Jik、JDii为分块子矩阵,diag表示对角矩阵;
式(7)中子矩阵Jik为:
式中,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部;分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部;分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部,分别为导纳矩阵元素的实部和虚部;
式(7)中子矩阵JDii为:
式(9)中元素分别表示如下:
如果节点i的中性点接地,则式(9)的第4行和第8行、第4列和第8列元素清零;
形成基本雅可比子矩阵J0后,平衡节点对应的雅可比矩阵行和列元素清零;
G、形成雅可比矩阵J;
在基本雅可比子矩阵J0的基础上追加与PV节点有关的附加雅可比矩阵元素形成完整的雅可比矩阵J如下:
式中,J0为(8n)×(8n)阶基本雅可比子矩阵;D为(8m)×(3m)阶子矩阵;F为(3m)×(8m)阶子矩阵;O1为(8n-8m)×(3m)阶零矩阵;O2为(3m)×(8n-8m)阶零矩阵;O3为(3m)×(3m)阶零矩阵;
式(14)中D为分块对角阵:
D=diag(D1,D2,…,Dm) (15)
式(15)中元素分别表示如下:
如果节点i的中性点接地,则式(16)的第4行和第8行元素清零;
式(14)中F为分块对角阵:
F=diag(F1,F2,…,Fm) (17)
式(17)中元素分别表示如下:
由基本雅可比子矩阵J0、对角矩阵D、对角矩阵F以及零矩阵O1、O2、O3形成完整的雅可比矩阵J;
其特征在于:还包括以下步骤:
I、形成雅可比矩阵有效行标号数组JR和有效列标号数组JC;
雅可比矩阵J中有许多整行元素全部为0的行或整列元素全部为0的列,为无效的行或列,进行计算时,需要去掉这些行和列,只取有效的行和列;用Matlab的sum函数和abs函数计算雅可比矩阵各行元素或各列元素绝对值之和;如果各行元素或各列元素绝对值之和为0,则雅可比矩阵对应行或对应列的元素全部为0,从而找到了矩阵J中元素都为0的行或元素都为0的列;形成雅可比矩阵有效行标号数组和有效列标号数组的具体步骤如下:
I1、令矩阵A为由雅可比矩阵J元素的绝对值组成的矩阵,即:
A=abs(J) (19)
式中,abs为Matlab中对矩阵元素取绝对值的函数;
I2、利用Matlab的函数sum求矩阵A各行或各列的和,得数组R和C:
R=sum(A,2) (20)
C=sum(A,1) (21)
式中,sum为Matlab中矩阵元素求和函数,sum的第2个参数为1时表示对各列元素求和,第2个参数为2时表示对各行元素求和;
I3、分别用式(22)和式(23)提取数组R和C中不为0元素的标号;
JR=find(R) (22)
JC=find(C) (23)
式中,find为Matlab中查找不为0的数组元素的标号;
J、解修正方程并修正状态变量;
潮流计算的修正方程为:
ΔW=-JΔX (24)
式中,J为雅可比矩阵;ΔW为偏差量列向量;ΔX为状态变量修正量列向量;
式(24)中ΔX为:
考虑去掉雅可比矩阵中元素全部为0的行和列后的修正方程为:
ΔWJR=-JJR,JCΔXJC (28)
式中,ΔWJR为按数组JR记录的标号提取偏差量列向量ΔW的元素形成的新列向量;JJR,JC为按数组JR记录的行标号和JC记录的列标号提取雅可比矩阵J的元素形成的新矩阵;ΔXJC为按数组JC记录的标号提取状态变量修正量ΔX的元素形成的新列向量;
利用Matlab的除法运算求解式(28)所示的修正方程,得到状态变量修正量ΔXJC:
ΔXJC=-JJR,JC\ΔWJR (29)
式中,\为Matlab的除法运算符号,其左侧是分母,右侧是分子;
按式(30)修正状态变量:
K、令t=t+1,转至步骤D;
L、输出节点及支路数据,结束。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910695403.5A CN110336287B (zh) | 2019-07-30 | 2019-07-30 | 一种基于雅可比元素提取的配电系统三相潮流计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910695403.5A CN110336287B (zh) | 2019-07-30 | 2019-07-30 | 一种基于雅可比元素提取的配电系统三相潮流计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110336287A CN110336287A (zh) | 2019-10-15 |
CN110336287B true CN110336287B (zh) | 2022-08-12 |
Family
ID=68148056
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910695403.5A Active CN110336287B (zh) | 2019-07-30 | 2019-07-30 | 一种基于雅可比元素提取的配电系统三相潮流计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110336287B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111355236B (zh) * | 2020-04-10 | 2023-05-30 | 大连海事大学 | 一种计及中性点电压变量的配电网三相潮流计算方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106356859A (zh) * | 2016-09-29 | 2017-01-25 | 大连海事大学 | 一种基于Matlab的直角坐标牛顿法潮流计算方法 |
-
2019
- 2019-07-30 CN CN201910695403.5A patent/CN110336287B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106356859A (zh) * | 2016-09-29 | 2017-01-25 | 大连海事大学 | 一种基于Matlab的直角坐标牛顿法潮流计算方法 |
Also Published As
Publication number | Publication date |
---|---|
CN110336287A (zh) | 2019-10-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Ghatak et al. | A fast and efficient load flow technique for unbalanced distribution system | |
CN106356859B (zh) | 一种基于Matlab的直角坐标牛顿法潮流计算方法 | |
Wei et al. | An innovative fixed-pole numerical approximation for fractional order systems | |
CN106532711B (zh) | 随迭代和节点类型改变雅可比矩阵的牛顿法潮流计算方法 | |
CN107204617B (zh) | 基于线性规划的直角坐标形式的区间潮流计算方法 | |
Wang et al. | Three-phase distribution power flow calculation for loop-based microgrids | |
US11853384B2 (en) | Methods of patel loadflow computation for electrical power system | |
CN107196306B (zh) | 基于Matlab稀疏矩阵的快速分解法潮流计算方法 | |
CN111355236B (zh) | 一种计及中性点电压变量的配电网三相潮流计算方法 | |
CN106229988B (zh) | 一种基于Matlab的极坐标牛顿法潮流计算方法 | |
CN106709243B (zh) | 含小阻抗支路电网的补偿法极坐标牛顿法潮流计算方法 | |
CN106532712B (zh) | 含小阻抗支路电网的补偿法直角坐标牛顿法潮流计算方法 | |
CN110336287B (zh) | 一种基于雅可比元素提取的配电系统三相潮流计算方法 | |
CN109617080A (zh) | 基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法 | |
CN110336288B (zh) | 基于矩阵运算提取雅可比元素的配电网三相潮流计算方法 | |
CN110417022A (zh) | 矩阵乘法运算提取雅可比元素的配电网三相潮流计算方法 | |
CN106410811B (zh) | 首次迭代小阻抗支路端点改变雅可比矩阵的潮流计算方法 | |
CN109494748A (zh) | 基于节点类型和修正的雅可比矩阵的牛顿法潮流计算方法 | |
Wang et al. | Load flow analysis | |
CN106712029B (zh) | 小阻抗支路pq端点变雅可比矩阵的牛顿法潮流计算方法 | |
CN113452028B (zh) | 低压配电网概率潮流计算方法、系统、终端和存储介质 | |
CN107181260B (zh) | 基于Matlab稀疏矩阵直角坐标牛顿法潮流计算方法 | |
CN114566969A (zh) | 一种适合研究目的使用的直角坐标牛顿法潮流计算方法 | |
CN106529089B (zh) | 用于含小阻抗支路电网的补偿法快速分解法潮流计算方法 | |
CN110417021A (zh) | 一种基于雅可比元素快速提取的配电网三相潮流计算方法 |
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 |