直角坐标系下牛顿法潮流计算

基于极坐标的牛顿拉夫逊潮流计算
极坐标牛顿拉夫逊潮流(PNF)计算也称为极坐标矢量法，是一种利用
极坐标变换计算复杂电力系统潮流的数值计算方法。

极坐标牛顿拉夫逊潮
流计算把复杂的三相系统潮流变换为一种简单的极坐标系统，把潮流计算
变换为水平和垂直运动，而每一个潮流方程只需要解决一个联立方程，从
而大大简化了潮流计算的复杂性。

极坐标牛顿拉夫逊潮流计算首先把潮流计算分为两部分，第一部分是
极坐标系变换，把复杂的三相系统潮流变换为一种简单的极坐标系统，这
部分计算主要包括把直流潮流变换为极坐标系统、设置极坐标参数等步骤；第二部分是牛顿拉夫逊潮流计算，即基于极坐标系统解决潮流方程。

牛顿
拉夫逊潮流计算的核心部分是极坐标潮流模型，包括水平潮流模型和垂直
潮流模型。

极坐标牛顿拉夫逊潮流计算的优势主要体现在潮流计算过程简单、求
解效率高。

直角坐标系下牛顿法潮流计算1电力系统潮流计算潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。

在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。

此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。

2节点导纳矩阵的形成在图1（a ）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b ）所示。

将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c ）的等值网络，其中1101I y E =和4404I y E =分别称为节点1和4的注入电流源。

(a)24İİ4(c)图1 电力系统及其网络以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下：1011212112212022323242423323434244234434044()()()()0()()0()()y U y U U I y U U y U y U U y U U y U U y U U y U U y U U y U I ⎫+-=⎪-++-+-=⎪⎬-+-=⎪⎪-+-+=⎭ （2-1） 上述方程组经过整理可以写成1111221211222233244322333344422433444400Y U Y U I Y U Y U Y U Y U Y U Y U Y U Y U Y U Y U I ⎫+ =⎪+++=⎪⎬++=⎪⎪ ++=⎭ （2-2）式中，111012Y y y =+；2220232412Y y y y y =+++；332334Y y y =+；44402434Y y y y =++；122112Y Y y ==-；233223Y Y y ==-；244224Y Y y ==-；344334Y Y y ==-。

课程设计说明书题目电力系统分析系 ( 部)专业( 班级 )姓名学号指导教师起止日期电力系统分析课程设计任务书系(部): 专业：指导教师：目录一、潮流计算基本原理1.1 潮流方程的基本模型1.2 潮流方程的讨论和节点类型的划分1.3、潮流计算的意义二、牛顿一拉夫逊法2.1 牛顿-拉夫逊法基本原理2.2节点功率方程2.3修正方程2.4 牛顿法潮流计算主要流程三、收敛性分析四、算例分析总结参考文献电力系统分析潮流计算一、潮流计算基本原理1.1潮流方程的基本模型电力系统是由发电机、变压器、输电线路及负荷等组成，其中发电机及负荷是非线性元件，但在进行潮流计算时，一般可以用接在相应节点上的一个电流注入量来代表。

因此潮流计算所用的电力网络系由变压器、输电线路、电容器、电抗器等静止线性元件所构成，并用集中参数表示的串联或并联等值支路来模拟。

结合电力系统的特点，对这样的线性网络进行分析，普通采用的是节点法，节点电压与节点电流之间的关系I=YV (1—1)其展开式为(i=1,2,3, …,n) (1—2)在工程实际中，已经的节点注入量往往不是节点电流而是节点功率，为此必须应用联系节点电流和节点功率的关系式 (i=1,2,3, …,n) (1—3)将 式 ( 1 - 3 ) 代 入 式 ( 1 - 2 ) 得 到 (i=1,2,3, …,n) (1-4)交流电力系统中的复数电压变量可以用两种极坐标来表示V =Vei8. (1-5)或 V=e+jf (1-6)而复数导纳为Y=G+jB (1-7)将式(1-6)、式(1- 7)代入以导纳矩阵为基础的式(1-4),并将实部与虚部分开，可以得到以下两种形式的潮流方程。

潮流方程的直角坐标形式为潮流方程的极坐标形式为(1—10)(1-11)以上各式中，j∈i表示乙号后的标号j的节点必须直接和节点i相联，并包括j=i的情况。

这两种形式的潮流方程通常称为节点功率方程，实牛顿一拉夫逊等潮流算法所采用的主要数学模型。

( k 1)fi(k ) f i
(k )
计算平衡节点功率 S s 和线路功率
~
停止
1 / 20
程序代码如下：
#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #include<conio.h> struct linetype // 线路参数 { int jiedian[2]; // 若为变压器，则 左端为“低”压侧（换算成Π型等值电路）,变压器的阻抗在“低”压侧 double R,X,K,B0; }line[30]; struct Nodetype // 节点功率 { int lei; // PQ 定义 1，PV 定义 2，平衡节点 定义 3 int jie; // 节点编号 double P,Q; double y1,y2; // 初始电压 }poin[30]; int point,road; // 节点数 point 支路数 road int p1,p2; // PQ PV 节点数目 //************************************************* 自定义 函数 *************************************************************** void chargescreen() // 调节屏幕 { int mode; printf("\t 请选择界面模式： ①. 106*45 ②. 134*45\n\t>>"); a: scanf("%d",&mode); if(mode!=1 && mode!=2) { printf("\n\t 错误，请重新输入...\n\t>>"); goto a; } printf("\n\t"); system("pause"); if(mode==1) system("mode con:cols=106 lines=45"); // 调整屏幕大小 else system("mode con:cols=134 lines=45"); } void pqpv() // 统计 PQ、PV 节点 数目

关于直角坐标牛拉法系统潮流分布计算的答辩关于直角坐标牛拉法系统潮流分布计算的答辩一、引言直角坐标牛拉法是一种常用的电力系统潮流计算方法，用于计算电力系统中各节点的电压和功率分布。

本文将从算法原理、计算步骤、应用场景等方面进行阐述和答辩。

二、算法原理直角坐标牛拉法基于功率平衡方程和节点电压方程，通过迭代求解的方式，逐步逼近系统的潮流分布。

其核心思想是将电压和功率分别表示为实部和虚部，通过复数运算来求解未知量。

具体而言，直角坐标牛拉法将电流和导纳分别表示为复数形式，利用复数的乘法和除法运算，将节点电流和导纳联系起来，从而得到节点电压和功率的计算结果。

三、计算步骤直角坐标牛拉法的计算步骤包括以下几个部分：1. 初始化：给定电网拓扑结构、节点导纳和负荷信息，初始化节点电压和功率。

2. 潮流计算：根据功率平衡方程和节点电压方程，通过迭代计算节点电压和功率。

具体而言，每次迭代中，首先根据节点电压和导纳计算节点电流；然后，根据节点电流和导纳计算节点电压；再根据节点电压和导纳计算节点功率。

通过多次迭代，直到收敛为止。

3. 收敛判断：判断节点电压和功率的迭代计算是否收敛。

一般来说，可以通过判断节点电压和功率的变化量是否小于设定的收敛阈值来进行判断。

若满足收敛条件，则停止迭代；否则，继续迭代。

4. 输出结果：输出最终的节点电压和功率分布结果。

根据需要，还可以输出其他相关信息，如潮流方向、线路功率损耗等。

四、应用场景直角坐标牛拉法广泛应用于电力系统潮流计算和分析。

具体而言，它可以用于以下几个方面：1. 网络规划：通过潮流计算，可以评估电力系统的稳定性和可靠性，为电网规划提供依据。

例如，可以通过潮流计算来确定新建变电站的容量和位置，优化电网结构。

2. 运行调度：在电力系统的日常运行中，潮流计算可以用于实时监测和调度。

通过潮流计算，可以了解各节点的电压和功率情况，及时发现问题并采取措施，确保电力系统的安全稳定运行。

3. 短路分析：在电力系统发生短路故障时，潮流计算可以用于分析故障电流的分布情况，确定故障点和故障线路，为故障处理和保护调整提供参考。

例3-5利用牛顿-拉夫逊法直角坐标方式计算例3-3所示网络潮流分布情况。

解：确定例3-3系统雅可比矩阵的维数。

系统有n = 5条母线（节点），采用直角坐标方法求解时组成2(n -1) =8个方程，J(i )维数为8×8。

按题意要求，该系统中，节点1为平衡节点，保持U 1=1+j0为定值，2，4，5为PQ 节点，3为PU 节点，U 3=1.05+j0。

（1）赋初值由已知可知平衡节点:111.0,0e f == 对PQ、PU节点赋电压初值：(0)(0)(0)(0)(0)(0)(0)(0)245245331.0,0, 1.05,0e e e f f f e f ========（2）求PQ 节点有功、无功不平衡量，PU 节点有功、电压不平衡量()()(){}55(0)(0)(0)(0)(0)(0)(0)(0)222222222211()()8.0 1.00 2.6783 1.0000.8928 1.00 1.7855 1.0008.0s s j jj j jj j j j j P P P P e GeB f f Gf B e ==∆=-=---+=--⨯+⨯-++-⨯-+-⨯-+=-⎡⎤⎣⎦∑∑()()(){}55(0)(0)(0)(0)(0)(0)(0)(0)222222222211(0)(0)(0)(0)(0)(0)333333333()()2.80 1.00028.4590 1.0009.9197 1.0019.8393 1.0 1.5()(s s j jj j jj j j j j s s j jj j jQ Q Q Q f GeB f e Gf B e P P P P e GeB f f G==∆=-=--++=---⨯+-⨯+++⨯++⨯=-⎡⎤⎣⎦∆=-=---∑∑()(){}()()55(0)(0)311(0)22(0)22(0)2(0)222333333(0)(0)(0)(0)(0)(0)(0)4444444444)4.4 1.05007.4580 1.0507.4580 1.0000 4.00851.05 1.0500()(j j j j j s s s s j jj j jj j f B e U U U U e f P P P P e GeB f f Gf B e ==+=-⨯++⨯-+-⨯-++=⎡⎤⎣⎦∆=-=-+=-+=∆=-=---+∑∑()()()(){}()55(0)1155(0)(0)(0)(0)(0)(0)(0)(0)444444444411)0 1.000.8928 1.007.4580 1.05011.9219 1.00 3.57111.0000.3729()()00 1.0009.9197 1.009j j j s s j jj j jj j j j j Q Q Q Q f GeB f e Gf B e =====-⨯+-⨯-+-⨯-+⨯-+-⨯-+=⎡⎤⎣⎦∆=-=--++=--⨯++⨯++∑∑∑∑()()(){}()()()(){}55(0)(0)(0)(0)(0)(0)(0)(0)5555555555119.4406 1.050147.9589 1.0039.6768 1.0 6.052()()0 1.0 3.7290 1.00 1.7855 1.000 3.57111.009.0856 1.000s s j jj j jj j j j j P P P P e GeB f f Gf B e ==⨯+-⨯++⨯=⎡⎤⎣⎦∆=-=---+=-⨯-⨯-+-⨯-++-⨯-+⨯-+=⎡⎤⎣⎦∑∑()()()(){}55(0)(0)(0)(0)(0)(0)(0)(0)5555555555110()()00 1.0049.7203 1.0019.8393 1.00039.6786 1.00108.5782 1.00.66s s j jj j jj j j j j Q Q Q Q f GeB f e Gf B e ==∆=-=--++=--⨯+⨯++⨯+++⨯+-⨯=⎡⎤⎣⎦∑∑（3）计算雅可比矩阵以节点2（PQ ）有功、无功功率和节点3（PU ）电压幅值分别对各节点电压实部、虚部求导为例，其他节点的求解过程略。

稳态分析课设目录1.任务书 (2)2.模型简介及等值电路 (3)3.修正方程的建立 (4)4程序流程图 (9)5. MATLAB程序编写 (10)6.结果分析 (16)7.设计总结 (18)8.参考文献 (19)1.任务书课程设计任务书2.、模型简介及等值电路电力网络接线如下图所示，各支路阻抗标幺值参数如下：Z 12=0.02+j0.06，Z 13=0.08+j0.24， Z 23=0.06+j0.18， Z 24=0.06+j0.12， Z 25=0.04+j0.12， Z 34=0.01+j0.03， Z 45=0.08+j0.24， k=1.1。

该系统中，节点1为平衡节点，保持1 1.060V j =+为定值；节点2、3、4都是PQ 节点，节点5为PV 节点，给定的注入功率分别为:20.200.20S j =+， 3-0.45-0.15S j =，40.400.05S j =--，50.500.00S j =-+，5 1.10V =。

各节点电压（初值）标幺值参数如下：节点12345U i （0）=e i （0）+j f i （0）1.06+j0.0 1.0+j0.0 1.0+j0.0 1.0+j0.0 1.1+j0.0 计算该系统的潮流分布。

计算精度要求各节点电压修正量不大于10-5。

3.修正方程的建立当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量1212,,,...,nnf f fe e e 由于平衡节点的电压向量是给定的，因此待求两共2(1)n 需要2(n-1)个方程式。

事实上，除了平衡节点的功率方程式在迭代过程中没有约束作用以外，其余每个节点都可以列出两个方程式。

对PQ 节点来说，is isQ P 和是给定的，因而可以写出()()0()()0i ijij i ijjijj isjj jj ij iijij ijjj ijj iisi jjj ij ip ff fe G e G e P B B Q Qf f fG e e G e B B ∈∈∈∈⎫∆=---+=⎪⎬∆=--++=⎪⎭∑∑∑∑ （3-2-1）对PV 节点来说，给定量是is is V P 和，因此可以列出2222()()0()0i is ijij i ij j ijj ji jj i j ii is i iff fe G e G e P P B B fV V e ∈∈⎫∆=---+=⎪⎬⎪∆=-+=⎭∑∑ （3-2-2）以直角坐标系形式表示3.1迭代推算式采用直角坐标时,节点电压相量及复数导纳可表示为:i i i ij ij ijV e jf Y G jB =+=+ (3-2-3)将以上二关系式代入上式中,展开并分开实部和虚部;假定系统中的第1,2,,m 号为P —Q 节点,第m+1,m+2,,n-1为P —V 节点,根据节点性质的不同,得到如下迭代推算式:⑴对于PQ 节点1111()()()()nni i i ij j ij j i ij j ij j j j n ni i i ij j ij j i ij j ij j j j P P e G e B f f G f B e Q Q f G e B f e G f B e ====⎫∆=---+⎪⎪⎬⎪∆=--++⎪⎭∑∑∑∑ (3-2-4) 1,2,,i m =⑵对于PV 节点112222()()()n ni i i ij j ij j i ij j ij j j j I i i i P P e G e B f f G f B e V V e f ==⎫∆=---+⎪⎬⎪∆=-+⎭∑∑ (3-2-5)1,2,,1i m m n =++-⑶对于平衡节点平衡节点只设一个,电压为已知,不参见迭代,其电压为:n n n V e jf =+ (3-2-6)3.2修正方程式(2-3-5)和(2-3-6)两组迭代式工包括2(n-1)个方程.选定电压初值及变量修正量符号之后代入式(2-3-5)和(2-3-6),并将其按泰勒级数展开,略去,i i e f ∆∆二次方程及以后各项,得到修正方程如下:x J f ∆=∆ (3-2-7)(3-2-8)3.3雅可比矩阵各元素的算式式(3-2-8)中, 雅可比矩阵中的各元素可通过对式(3-2-4)和(3-2-5)进行偏导而求得.当j i ≠时, 雅可比矩阵中非对角元素为22()0i iij i ij i j j i i ij i ij i j jj j P Q G e B f e f P Q B e G f f e U U e f ⎫∂∆∂∆=-=-+⎪∂∆∂∆⎪⎪∂∆∂∆⎪==-⎬∂∆∂∆⎪⎪∂∆∂∆⎪==∂∂⎪⎭(3-2-9)当j i =时,雅可比矩阵中对角元素为:111122()()()()22ni ij j ij j ii i ii i j i n iij j ij j ii i ii i j j n iij j ij j ii i ii i j i niij j ij j ii i ii i j j i iji i i P G e B f G e B f e P G f B e G f B e f Q G f B e G f B e e Q G e B f G e B f f U e e U f f ====∂∆⎫=----⎪∂⎪⎪∂∆=-+-+⎪∂⎪⎪∂∆⎪=+-+∂⎪⎪⎬∂∆⎪=-∆-++⎪∂⎪∂∆⎪=-⎪∂⎪∂∆=-∂⎭∑∑∑∑⎪⎪⎪ (3-2-10)雅可比矩阵各元素的表示如下:()()()()ij i ij i i ij ij j ij j ii i ii i j j iG e B f j i P H G e B f G e B f j i e ∈-+⎧≠∂∆⎪==⎨----=∂⎪⎩∑)()()()ij i ij i i ij ij j ij j ii i ii i j j iB e G f j i P N G f B e B e G f j i f ∈-⎧≠∂∆⎪==⎨-++-=∂⎪⎩∑)()()()ij i ij i i ij ij j ij j ii i ii i j j i B e G f j i Q M G f B e B e G f j i e ∈-⎧≠∂∆⎪==⎨++-=∂⎪⎩∑)()()()ij i ij i i ij ij j ij j ii i ii i j j i G e B f j i Q L G e B f G e B f j i f ∈+⎧≠∂∆⎪==⎨--++=∂⎪⎩∑20()2()i ij ij j i U R e j i e ≠⎧∂∆==⎨-=∂⎩4.程序流程图启动输入数据形成节点导纳矩阵设定节点起始计算电压迭代次数k=0应用式(7-16)和(7-17)计算置节点号i=1雅克比矩阵J 是否已经全部形成，i>n?按式(7-21)和(7-22)计算雅克比矩阵元素J 增大节点号i=i+1解修正方程式,由计算电压修正量求出迭代是否收敛 ,计算平衡节点功率 和线路功率结束S s S ijV i (k)maxi(k)max3?<V i (k)V i (k)V i (k)V i (0)V i(k)V i (k+1)i(k)i(0)i(k)i(k)i(k)i(k+1),和P i (k)(k)(k)i iQ J 和(k)iP i(k)(k)iQ ,=+=+k=k+1否是是否5.MATLAB程序编写程序编写如下：clc;clear;g(4,1)=2.7500;b(4,1)=-8.2500;g(4,3)=1.2500;b(4,3)=-3.7500;g(1,4)=2.7500;b(1,4)=-8.2500;g(1,2)=1.6667;b(1,2)=-5.0000;g(1,3)=3.3334;b(1,3)=-6.6667;g(1,5)=5.0000;b(1,5)=-15.0000;g(2,1)=1.6667;b(2,1)=-5.0000;g(2,3)=10.0000;b(2,3)=-30.0000;g(2,5)=1.2500;b(2,5)=-3.7500;g(3,4)=1.2500;b(3,4)=-3.7500;g(3,1)=3.3334;b(3,1)=-6.6667;g(3,2)=10.0000;b(3,2)=-30.0000;g(5,1)=5.0000;b(5,1)=-15.0000;g(5,2)=1.2500;b(5,2)=-3.7500;b(1,1)=-0.8250;b(4,4)=0.7500;g(1,1)=0.275;g(4,4)=-0.25;N1=4;for m=1:N1+1;for n=1:N1+1;if m==nG(m,m)=g(m,1)+g(m,2)+g(m,3)+g(m,4)+g(m,5);B(m,m)=b(m,1)+b(m,2)+b(m,3)+b(m,4)+b(m,5);elseG(m,n)=-g(m,n);B(m,n)=-b(m,n);endendendY=G+j*B;e(1)=1.0;e(2)=1.0;e(3)=1.0;e(4)=1.10;f(1)=0;f(2)=0;f(3)=0;f(4)=0;u(1)=1.0;u(2)=1.0;u(3)=1.0;u(4)=1.1;uu(1)=0;uu(2)=0;uu(3)=0;uu(4)=0;P(1)=0.20;Q(1)=0.20;P(2)=-0.45;Q(2)=-0.15;P(3)=-0.40;Q(3)=-0.05;P(4)=-0.50;Q(4)=0.00;k=0;disp('迭代前的参数：')e,f,uuN1=4;precision=1;%除平衡节点外节点数while precision>0.00001;e(5)=1.06;f(5)=0.00;for m=1:N1for n=1:N1+1pt(n)=e(m)*(G(m,n)*e(n)-B(m,n)*f(n))+f(m)*(G(m,n)*f(n)+B(m,n)*e(n));endpp(m)=P(m)-sum(pt);endfor m=1:N1-1for n=1:N1+1qt(n)=f(m)*(G(m,n)*e(n)-B(m,n)*f(n))-e(m)*(G(m,n)*f(n)+B(m,n)*e(n));endqq(m)=Q(m)-sum(qt);endfor m=N1uu2(m)=u(m)*u(m)-(e(m)*e(m)+f(m)*f(m))enda=1;for m=1:4dW(a)=pp(m);a=a+2;enda=2;for m=1:3dW(a)=qq(m);a=a+2;enda=3*2+2;for m=4,dW(a)=uu2(m);a=a+2;if max(dW)< precisionfprintf('\n 迭代是收敛的。

Sb SG STc S0c jQB 2 jQB3
1 b Tb 2 c Tc 3
A
d Td
SLDb
G
SG
SL D d
14
二、两级电压的开式电力网计算 计算方法一：包含理想变压器，计算时，经过理 想变压器功率保持不变，两侧电压之比等于实际 变比k。 T b d c L-1 L-2 SLD A
V1 arctg V1 V1
4
网络元件的功率损耗
功率损耗包括：电阻和等值电抗上的损耗 对地等值导纳上产生的损耗
V1S1 , I1 S ' I
jQB1
B j 2
R jX
S '', I S 2 , I 2 V2
jQB 2
B j 2
线路
VS1 , I1
线路
S0 (GT jBT )V 2
I0% S0 P0 jQ0 P0 j SN 100
开式网络的电压和功率分布计算
一、已知供电点电压和负荷点功率时的计算方法 已知末端的功率和电压：从末端开始依次计算出 电压降落和功率损耗。 已知电源点的电压和负荷的功率：采取近似的方 法通过叠代计算求得满足一定精度的结果
X2 k2 X2
T
A
A
B2 B2 / k 2 d c L-2 SLD
R'2+ j X'2 j B'2/2
16
R1+ jX1
j B1/2 j B1/2
b ΔS0
Z'T
c' j B'2/2
d'
SLD
二、两级电压的开式电力网计算 计算方法三：用π型等值电路代表变压器

极坐标系下的潮流计算
潮流计算
在电力系统中，潮流计算是一种常用的计算方法，用于确定在给定网络结构和参数下，各节点的电压 、电流和功率分布。在极坐标系下进行潮流计算，可以更好地描述和分析电力系统的电磁场分布和变 化。
极坐标系下的潮流计算特点
在极坐标系下进行潮流计算，可以更直观地描述电力线路的走向和角度变化，更好地反映电力系统的 复杂性和实际情况。此外，极坐标系下的潮流计算还可以方便地处理电力系统的非对称性和不对称故 障等问题。
03
CATALOGUE
极坐标系下的牛顿-拉夫逊算法
极坐标系简介
极坐标系
一种二维坐标系统，由一个原点（称为极点）和一条从极点出发的射线（称为 极轴）组成。在极坐标系中，点P的位置由一个角度θ和一个距离r确定。
极坐标系的应用
极坐标系广泛应用于物理学、工程学、经济学等领域，特别是在电力系统和通 信网络中，用于描述电场、磁场、电流和电压等物理量的分布和变化。
极坐标形式
将电力系统的节点和支路参数以极坐 标形式表示，将实数问题转化为复数 问题，简化计算过程并提高计算效率 。
02
CATALOGUE
牛顿-拉夫逊算法原理
算法概述
牛顿-拉夫逊算法是一种迭代算法，用于求解非线性方程组。在电力系统中，它 被广泛应用于潮流计算，以求解电力网络中的电压、电流和功率等参数。
准确的结果。
通过极坐标系的处理，算法 能够更好地处理电力系统的 复杂结构和不对称性，提高 了计算的准确性和适应性。
算例分析表明，该算法在处理 大规模电力系统时仍具有较好 的性能，能够满足实际应用的
需求。
展望
进一步研究牛顿-拉夫逊算法在极坐标 系下的收敛性分析，探讨收敛速度与电 力系统规模、结构和参数之间的关系， 为算法的优后的电压、电流和功 率等参数。

直角坐标系的计算机潮流算法实验程序：%直角坐标法求解潮流计算Branch=[1 2 0.1+0.4j 0.01528j 1 03 1 0.3j inf 1.1 11 4 0.12+0.50j 0.01920j 1 02 4 0.08+0.40j 0.01413j 1 0]%Branch矩阵：1、支路首端号；2、支路末端号；3、支路阻抗；4、支路对地导纳；5、支路的变化；6、支路首端处于K侧为1,1侧为0Y=zeros(4); %节点导纳矩阵for i=1:4if Branch(i,6)==0 %不含变压器的支路j=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-1/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3)+Branch(i,4);Y(k,k)=Y(k,k)+1/Branch(i,3)+Branch(i,4);elsej=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-Branch(i,5)/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3);Y(k,k)=Y(k,k)+Branch(i,5)^2/Branch(i,3);endenddisp('节点导纳矩阵');YG=real(Y);B=imag(Y);V=[1;1;1.1;1.05];%给定V的初始计算值e=real(V);f=imag(V);disp('节点电压的实部：')edisp('节点电压的虚部：')fdisp('节点注入有功功率：')Ps=[-0.3;-0.55;0.5;0]disp('节点注入无功功率：')Qs=[-0.18;-0.13;0;0]%由各节点电压向量（状态变量）可得各节点注入功率：P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];n=0;%辅助循环计数变量while (any(del_W>1e-5))&(n<5)n=n+1;disp('迭代次数：');disp(n)%--------------------------------------------------------%修正方程式：% [-0.30-P(1) ] [del_e(1)]% [-0.18-Q(1) ] [del_f(1)]% [-0.55-P(2) ] [del_e(2)]% [-0.13-Q(2) ] + J * [del_f(2)] = 0 % [0.5-P(3) ] [del_e(3)]% [1.1^2-e(3)^2-f(3)^2] [del_f(3)]%雅克比矩阵为：J=[P1_e1 P1_f1 P1_e2 P1_f2 P1_e3 P1_e3;% Q1_e1 Q1_f1 Q1_e2 Q1_f2 Q1_e3 Q1_f3% P2_e1 P2_f1 P2_e2 P2_f2 P2_e3 P2_f3% Q2_e1 Q2_f1 Q2_e2 Q2_f2 Q2_e3 Q2_f3% P3_e1 P3_f1 P3_e2 P3_f2 P3_e3 P3_f3% V3_e1 V3_f1 V3_e2 V3_f2 V3_e3 V3_f3]%--------------------------------------------------------%-------------------求雅可比矩阵的参数-----------------------%GeBf=G*e-B*f;%辅助计算函数GfBe=G*f+B*e;%辅助计算函数for i=1:2for j=1:3if i==jP_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);elseP_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);Q_f(i,j)=-P_e(i,j);P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);Q_e(i,j)=P_f(i,j);endendendfor j=1:2P_e(3,j)=-G(3,j)*e(3)-B(3,j)*f(3);P_f(3,j)=B(3,j)*e(3)-G(3,j)*f(3);V_e(3,j)=0;V_f(3,j)=0;endP_e(3,3)=-GeBf(3)-G(3,3)*e(3)-B(3,3)*f(3);P_f(3,3)=-GfBe(3)-G(3,3)*f(3)+B(3,3)*e(3);V_e(3,3)=-2*e(3);V_f(3,3)=-2*f(3);%----------------------------------------------------%%--------------------求出雅克比矩阵-------------------%J=zeros(6,6);for i=1:3for j=1:3J(2*i-1,2*j-1)=P_e(i,j);J(2*i-1,2*j)=P_f(i,j);endendfor i=1:2for j=1:3J(2*i,2*j-1)=Q_e(i,j);J(2*i,2*j)=Q_f(i,j);endendfor j=1:3J(6,2*j-1)=V_e(3,j);J(6,2*j)=V_f(3,j);enddisp('雅可比矩阵为：')J%-----------------------------------------------------%%------------解修正方程并得出修正后的电压向量-----------% del_V=-J\del_W;clear i;for k=1:3V(k)=V(k)+del_V(2*k-1)+i*del_V(2*k);enddisp('节点电压为：')disp(V)%-----------------------------------------------------%%------------------计算节点不平衡量--------------------%e=real(V);f=imag(V);P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];disp('节点不平衡量为：')disp(del_W)end%------------------------最终结论-------------------------%disp('最终各节点电压幅值为：')disp(abs(V))disp('最终各节点电压相角(度)为：')disp(180*angle(V)/pi)disp('最终各节点注入功率：')S=P+i*Qdisp('最终各节点注入有功功率为：')Pdisp('最终各节点注入无功功率为：')Q实验运行结果如下：Branch =1.00002.0000 0.1000 + 0.4000i 0 + 0.0153i 1.0000 03.0000 1.0000 0 + 0.3000i Inf 1.1000 1.00001.0000 4.0000 0.1200 + 0.5000i 0 + 0.0192i1.0000 02.0000 4.0000 0.0800 + 0.4000i 0 + 0.0141i 1.0000 0节点导纳矩阵Y =1.0421 - 8.2429i -0.5882 +2.3529i 0 +3.6667i -0.4539 + 1.8911i-0.5882 + 2.3529i 1.0690 - 4.7274i 0 -0.4808 + 2.4038i0 + 3.6667i 0 0 - 3.3333i 0-0.4539 + 1.8911i -0.4808 + 2.4038i 0 0.9346 - 4.2616i节点电压的实部：e =1.00001.10001.0500节点电压的虚部：f =节点注入有功功率：Ps =-0.3000-0.55000.5000节点注入无功功率：Qs =-0.1800-0.1300迭代次数：1雅可比矩阵为：J =-1.0194 -8.3719 0.5882 2.3529 0 3.6667 -8.1138 1.0648 2.3529 -0.5882 3.6667 00.5882 2.3529 -1.0450 -4.8770 0 02.3529 -0.5882 -4.5778 1.0930 0 00 4.0333 0 0 0 -3.66670 0 0 0 -2.2000 0节点电压为：0.9935 - 0.0088i0.9763 - 0.1078i1.1000 + 0.1267i1.0500节点不平衡量为：-0.0013-0.0028-0.0135-0.0547-0.0160迭代次数：2雅可比矩阵为：J =-0.8091 -8.3613 0.6052 2.3325 0.0324 3.6429 -7.9992 1.4071 2.3325 -0.6052 3.6429 -0.03240.8280 2.2338 -1.0190 -4.6364 0 02.2338 -0.8280 -4.3641 2.0878 0 0-0.4644 4.0333 0 0 -0.0324 -3.64290 0 0 0 -2.2000 -0.2533节点电压为：0.9847 - 0.0086i0.9590 - 0.1084i1.0924 + 0.1289i1.0500节点不平衡量为：-0.0000-0.0000-0.0003-0.00110.0001-0.0001迭代次数：3雅可比矩阵为：J =-0.7940 -8.2936 0.5995 2.3120 0.0315 3.6107 -7.9228 1.4000 2.3120 -0.5995 3.6107 -0.03150.8191 2.1927 -0.9865 -4.6144 0 02.1927 -0.8191 -4.2210 2.0885 0 0-0.4728 4.0056 0 0 -0.0315 -3.61070 0 0 0 -2.1849 -0.2579节点电压为：0.9846 - 0.0086i0.9587 - 0.1084i1.0924 + 0.1290i1.0500节点不平衡量为：1.0e-006 *-0.0069-0.0077-0.38310.0099-0.0014最终各节点电压幅值为：0.98470.96481.10001.0500最终各节点电压相角(度)为：-0.5002-6.45036.7323最终各节点注入功率：S =-0.3000 - 0.1800i-0.5500 - 0.1300i0.5000 + 0.0934i0.3679 + 0.2647i最终各节点注入有功功率为：P =-0.3000-0.55000.50000.3679最终各节点注入无功功率为：Q =-0.1800-0.13000.09340.2647。

%主程序"PowerFlow_NR.m"function [bus_res,S_res] = PowerFlow_NR_2 % 牛顿-拉夫逊法解潮流方程的主程序[bus,line] = OpDF_; % 打开数据文件的子程序，返回bus（节点数据）和line（线路数据）回主程序[nb,mb]=size(bus);[nl,ml]=size(line); % 计算bus和line矩阵的行数和列数[bus,line,nPQ,nPV,nodenum] = Num_(bus,line); % 对节点重新排序的子程序Y = y_(bus,line); % 计算节点导纳矩阵的子程序myf = fopen('Result.m','w');fprintf(myf,'\n\n');fclose(myf); % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵format longEPS = 1.0e-10; % 设定误差精度for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出[dP,dQ] = dPQ_(Y,bus,nPQ,nPV); % 计算功率偏差dP和dQ的子程序J = Jac_(bus,Y,nPQ); % 计算雅克比矩阵的子程序UD = zeros(nPQ,nPQ);for i = 1:nPQUD(i,i) = bus(i,2); % 生成电压对角矩阵enddAngU = J \ [dP;dQ];dAng = dAngU(1:nb-1,1); % 计算相角修正量dU = UD*(dAngU(nb:nb+nPQ-1,1)); % 计算电压修正量bus(1:nPQ,2) = bus(1:nPQ,2) - dU; % 修正电压bus(1:nb-1,3) = bus(1:nb-1,3) - dAng; % 修正相角if (max(abs(dU))<EPS)&(max(abs(dAng))<EPS)breakend % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代endbus = PQ_(bus,Y,nPQ,nPV); % 计算每个节点的有功和无功注入的子程序[bus,line] = ReNum_(bus,line,nodenum); % 对节点恢复编号的子程序YtYm = YtYm_(line); % 计算线路的等效Yt和Ym的子程序，以计算线路潮流bus_res = bus_res_(bus); % 计算节点数据结果的子程序S_res = S_res_(bus,line,YtYm); % 计算线路潮流的子程序myf = fopen('Result.m','a');fprintf(myf,'---------------牛顿－拉夫逊法潮流计算结果----------\n\n 节点计算结果：\n节点节点电压节点相角（角度）节点注入功率\n');for i = 1:nbfprintf(myf,'%2.0f ',bus_res(i,1));fprintf(myf,'.6f ',bus_res(i,2));fprintf(myf,'.6f ',bus_res(i,3));fprintf(myf,'.6f + j .6f\n',real(bus_res(i,4)),imag(bus_res(i,4)));endfprintf(myf,'\n线路计算结果：\n节点I 节点J 线路功率S(I,J)线路功率S(J,I) 线路损耗dS(I,J)\n');for i = 1:nlfprintf(myf,'%2.0f ',S_res(i,1));fprintf(myf,'%2.0f ',S_res(i,2));fprintf(myf,'.6f + j .6f ',real(S_res(i,3)),imag(S_res(i,3)));fprintf(myf,'.6f + j .6f ',real(S_res(i,4)),imag(S_res(i,4)));fprintf(myf,'.6f + j .6f\n',real(S_res(i,5)),imag(S_res(i,5)));endfclose(myf); % 迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束bus = [1 1.00 0.00 -0.30 -0.18 1;2 1.00 0.00 -0.55 -0.13 1;3 1.10 0.00 0.50 0.00 2;4 1.05 0.00 0.00 0.00 3]line = [1 2 0.10 0.40 0.0 0.01528 0;4 2 0.08 0.40 0.0 0.01413 0;1 4 0.12 0.50 0.0 0.0192 0;3 1 0.00 0.3 0.0 0.0 -1.1]----------------------------------------牛顿－拉夫逊法潮流计算结果-------------------------------------------- 节点计算结果：节点节点电压节点相角（角度）节点注入功率1 0.984674906330845 -0.500170385513657 -0.300000000000000 - 0.180000000000000i2 0.964797665550885 -6.450305258622626 -0.550000000000000 - 0.130000000000000i3 1.100000000000000 6.732349388989963 0.500000000000000 + 0.093411003244513i4 1.050000000000000 0 0.367882692523292 + 0.264698252215732i 线路计算结果：节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)1 2 0.246244-j0.014651 -0.239990+j0.010627 0.006254- j0.0040241 3 -0.500001-j0.029264 0.500000+j0.093409 -0.000001+j0.0641451 4 -0.046244-j0.136088 0.048216+j0.104522 0.001972-j0.0315662 4 -0.310010-j0.140627 0.319666+j0.160176 0.009656+j0.019549%子程序1 "OpDF_.m" 作用为打开数据文件function [bus,line] = OpDF_[dfile,pathname]=uigetfile('*.m','Select Data File'); % 数据文件类型为m文件，窗口标题为“Select Data File”if pathname == 0error(' you must select a valid data file') % 如果没有选择有效文件，则出现错误提示elselfile =length(dfile); % 读取文件字符串长度eval(dfile(1:lfile-2)); % 去除后缀，打开文件！注意！新浪博客中"eval"函数会自动加上"_r"后缀，此处的函数名是"eval"而不是"eval_r"，拷贝后请去除"_r"后缀end%子程序2 "Num.m" 作用为对节点重排序，并修改相应的线路数据function [bus,line,nPQ,nPV,nodenum] = Num_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);nSW = 0; % nSW为平衡节点个数nPV = 0; % nPV为PV节点个数nPQ = 0; % nPQ为PQ节点个数for i = 1:nb, % nb为总节点数type= bus(i,6);if type == 3,nSW = nSW + 1;SW(nSW,:)=bus(i,:); % 计算并储存平衡节点elseif type == 2,nPV = nPV +1;PV(nPV,:)=bus(i,:); % 计算并储存PV节点elsenPQ = nPQ + 1;PQ(nPQ,:)=bus(i,:); % 计算并储存PQ节点endendbus=[PQ;PV;SW]; % 对bus矩阵按PQ、PV、平衡节点的顺序重新排序nodenum=[[1:nb]' bus(:,1)];% 生成新旧节点对照表for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendend % 按排序以后的节点顺序对line矩阵重新编号%子程序3 "y_.m" 作用为计算节点导纳矩阵function Y = y_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb); % 对导纳矩阵赋初值0for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4); % 读入线路参数if Zt~=0Yt=1/Zt; % 接地支路不计算YtendYm=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0Y(I,I)=Y(I,I)+Ym;endif K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/K/K;Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endif K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endend%子程序4 "dPQ_.m" 作用为计算功率偏差function [dP,dQ] =dPQ_(Y,bus,nPQ,nPV) % nPQ、nPV为相应节点个数n = nPQ + nPV +1; % 总节点个数dP = bus(1:n-1,4);dQ = bus(1:nPQ,5); % 对dP和dQ赋初值PV节点不需计算dQ 平衡节点不参与计算for i = 1:n-1for j = 1:ndP(i,1) =dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));if i<nPQ+1dQ(i,1) = dQ(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endendend % 利用循环计算求取dP和dQ%子程序5 "Jac_.m" 作用为计算雅克比矩阵function J = Jac_(bus,Y,nPQ)[nb,mb]=size(bus);H = zeros(nb-1,nb-1);N = zeros(nb-1,nPQ);K = zeros(nPQ,nb-1);L = zeros(nPQ,nPQ); % 将雅克比矩阵分块，即：J = [H N; K L]，并初始化Qi = zeros(nb-1,1);Pi = zeros(nb-1,1);for i = 1:nb-1for j = 1:nbPi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3) ));Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3) ));endend % 初始化并计算Qi和Pifor i = 1:nb-1for j = 1:nb-1if i~=jH(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseH(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);end % 分别计算H矩阵的对角及非对角元素if j < nPQ+1if i~=jN(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseN(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);endend % 分别计算N矩阵的对角及非对角元素if i < nPQ+1if i~=jK(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseK(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);end % 分别计算K矩阵的对角及非对角元素if j < nPQ+1if i~=jL(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseL(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);endend % 分别计算L矩阵的对角及非对角元素endendendJ = [H N; K L]; % 生成雅克比矩阵%子程序6 "PQ_.m" 作用为计算每个节点的功率注入function bus = PQ_(bus,Y,nPQ,nPV)n = nPQ+nPV+1; % n为总节点数for i = nPQ+1:n-1for j = 1:nbus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j, 3)));endend % 利用公式计算PV节点的无功注入for j =1:nbus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-b us(j,3)));bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-imag(Y(n,j))*cos(bus(n,3)-b us(j,3)));end % 计算平衡节点的无功及有功注入%子程序7 "ReNum_.m" 作用为对节点和线路数据恢复编号function [bus,line] = ReNum_(bus,line,nodenum)[nb,mb]=size(bus);[nl,ml]=size(line);bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus矩阵的数据k = 1;for i = 1 :nbfor j = 1 : nbif bus(j,1) == kbus_temp(k,:) = bus(j,:);k = k + 1;endendend % 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,1)line(i,j)=nodenum(k,2);breakendendendend % 恢复line的编号%子程序8 "YtYm_.m" 作用为计算线路的等效Yt和Ym，以计算线路潮流function YtYm = YtYm_(line)[nl,ml]=size(line);YtYm = zeros(nl,5); % 对YtYm矩阵赋初值0YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ymj = sqrt(-1);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);if Zt~=0Yt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=0YtYm(k,3) = Yt;YtYm(k,4) = Ym;YtYm(k,5) = Ym;endif (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0YtYm(k,4) = Ym;endif K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧YtYm(k,3) = Yt/K;YtYm(k,4) = Ym+Yt*(K-1)/K;YtYm(k,5) = Yt*(1-K)/K/K;endif K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧YtYm(k,3) = -Yt*K;YtYm(k,4) = Ym+Yt*(1+K);YtYm(k,5) = Yt*(K^2+K);endend%子程序9 "bus_res_.m" 计算并返回节点数据结果function bus_res = bus_res_(bus);[nb,mb]=size(bus);bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果bus_res(:,1:2) = bus(:,1:2);bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率%子程序10 "S_res_.m" 计算并返回线路潮流function S_res = S_res_(bus,line,YtYm)[nl,ml]=size(line);S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果S_res(:,1:2) = line(:,1:2); % 前两列为节点编号for k=1:nlI = S_res(k,1);J = S_res(k,2);if (J~=0)&(I~=0)S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(cos(bus(I,3))+j*sin(bu s(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtYm(k,3));S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*s in(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3));S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流elseif J==0S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));S_res(k,5)=S_res(k,3)+S_res(k,4);elseS_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流endend。

况下 ，尽管极 坐标 牛顿法 中的方程数和迭代次数可能均少 于直角 坐标牛顿法 ，但无论 系统 中 ＰＶ节点数 多少 ，直角
坐标 牛顿法的潮流计算速度均快于极坐标 牛顿法。对 ＩＥＥＥ一３０、一５７、一１１８节点 系统进行 编程计算结 果表 明 ，直
角坐标牛顿法 的计算效率远高 于极 坐标 牛顿法 ，尤其 在运 用稀 疏矩阵技术后 。
第 ３８卷第 １期 ２０１６年 ３月
南昌大学学报 （工科版 ） Ｊｏｕｒ ｎａｌ ｏｆ Ｎａｎｃｈａｎｇ Ｕｎｉｖｅｒｓｉｔｙ（Ｅｎｇｉｎｅｅｒ ｉｎｇ＆ Ｔｅｃｈｎｏｌｏｇｙ）
文 章 编 号 ：１００６—０４５６（２０１６）０１—００９８—０５
Ｖｏ１．３８ Ｎｏ．１ Ｍ ａｒ．２０１６
ＳＨＡＯ Ｗｅｉｚｈｅ，ＷＡＮＧ Ｙｕｊｕｎ，ＷＡＮ Ｘｉｎｒｕ，ＧＯＮＧ Ｊｉａｗｅｉ，ＣＨＥＮ Ｋｅｎ
（Ｓｃｈｏｏｌ ｏｆ Ｉｎｆｏｒｍａｔｉｏｎ Ｅｎｇｉｎｅｅｒ ｉｎｇ，Ｎａｎｃｈａｎｇ Ｕｎｉｖｅｒｓｉｔｙ，Ｎａｎｃｈａｎｇ ３３００３ １，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｔｈｅ ｃｏｍｐａｒｉｓｏｎ ｏｆ ｔｈｅ ｃａｌｃｕｌａｔｉｎｇ ｓｐｅｅｄｓ ｏｎ Ｎｅｗｔｏｎ—Ｒａｐｈｓｏｎ ｐｏｗｅｒ ｆｌｏｗ ｍｅｔｈｏｄ ｉｎ ｒｅｃｔａｎｇｕｌａｒ ａｎｄ ｐｏｌａｒ ｃｏｏｒｄｉｎａｔｅｓ ｗａｓ ｐｒｅｓｅｎｔｅｄ ｉｎ ｔｈｅ ｐａｐｅｒ，ｗｈｉｃｈ ｉｎｃｌｕｄｅｄ ｔｈｅ ａｍｏｕｎｔ ｏｆ ｔｈｅ ｅｑｕａｔｉｏｎｓ ｓｏｌｖｅｄ，ｔｈｅ ｃａｌｃｕｌａｔｉｎｇ ａ— ｍｏｕｎｔ ｏｆ ｔｈｅ ｅｌｅｍｅｎｔｓ ｉｎ Ｊａｃｏｂｉａｎ ｍａｔｒｉｘ，ｔｈｅ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｔｒ ｉａｎｇｌｅ ｆｕｎｃｔｉｏｎ，ｔｈｅ ｉｔｅｒａｔｉｏｎ ｎｕｍｂｅｒ，ｔｈｅ ａｍｏｕｎｔ ｏｆ ＰＶ ｂｕｓｅｓ ｉｎ ｔｈｅ ｓｙｓｔｅｍ ａｎｄ ＳＯ ｏｎ．Ｉｔ ｗａｓ ｏｂｖｉｏｕｓ ｔｈａｔ ｔｈｅ ｎｕｍｂｅｒ ｏｆ ＰＶ ｂｕｓｅｓ ｉｎ ｔｈｅ ｓｙｓｔｅｍ ｐｌａｙｅｄ ａ ｖｉｔａｌ ｒｏｌｅ ｉｎ ｄｅｃｉ－ ｄｉｎｇ ｔｈｅ ｃａｌｃｕｌａｔｉｎｇ ｓｐｅｅｄｓ ｗｉｔｈ Ｎｅｗｔｏｎ ｍｅｔｈｏｄ ｉｎ ｒｅｃｔａｎｇｕｌａｒ ａｎｄ ｐｏｌａｒ ｃｏｏｒｄｉｎａｔｅｓ．Ｗ ｉｔｈｏｕｔ ｃｏｎｓｉｄｅｒｉｎｇ ｔｈｅ ｅｌｅ— ｍｅｎｔｓ’ ｓｐａｒｓｉｔｙ，ｉｆ ｔｈｅ ｓｕｐｐｏｓｉｎｇ ｔｈａｔ ｔｈｅ ａｍｏｕｎｔ ｏｆ ＰＶ ｂｕｓｅｓ ｉｎ ｔｈｅ ｓｙｓｔｅｍ ｗａｓ ｌｏｗ ，ｔｈｅ ｃａｌｃｕｌａｔｉｎｇ ｓｐｅｅｄ ｏｆ Ｎｅｗｔｏｎ ｍｅｔｈｏｄ ｉｎ ｒｅｃｔａｎｇｕ ｌａｒ ｃｏｏｒｄｉｎａｔｅｓ ｗｏｕｌｄ ｂｅ ｆａｓｔｅｒ ｔｈａｎ ｔｈｅ ｓｐｅｅｄ ｉｎ ｐｏｌａｒ ｃｏｏｒｄｉｎａｔｅｓ．Ｏｎ ｔｈｅ ｃｏｎｔｒａｒｙ，ｓｕｐｐｏｓｅ ｔｈａｔ ｔｈｅ ａｍ ｏｕｎｔ ｏｆ ＰＶ ｂｕｓｅｓ ｉｎ ｔｈｅ ｓｙｓｔｅｍ ｗａｓ ｈｉｇｈ，ｔｈｅｎ ｉｔ ｗｏｕｌｄ ｌｅａｄ ｔｏ ａ ｏｐｐｏｓｉｔｅ ｒｅｓｕｌｔ．Ｗ ｉｔｈ ｔｈｅ ｅｌｅｍｅｎｔｓ’ ｓｐａｒｓｉｔｙ ｃｏｎｓｉｄｅｒｅｄ，ｅｖｅｎ ｔｈｏｕｇｈ ｔｈｅ ａｍｏｕｎｔ ｏｆ ｔｈｅ ｅｑｕａｔｉｏｎｓ ｓｏｌｖｅｄ ａｎｄ ｔｈｅ ｉｔｅｒａｔｉｏｎ ｎｕｍ ｂｅｒ ｗｉｔｈ Ｎｅｗｔｏｎ ｍ ｅｔｈｏｄ ｉｎ ｐｏｌａｒ ｃｏｏｒｄｉｎａｔｅｓ ｗａｓ ｌｅｓｓ ｔｈａｎ ｗｉｔｈ Ｎｅｗｔｏｎ ｍｅｔｈｏｄ ｉｎ ｒｅｃｔａｎｇｕｌａｒ ｃｏｏｒｄｉｎａｔｅｓ．Ｂｕｔ ｗｈａｔｅｖｅｒ ｔｈｅ ａｍｏｕｎｔ ｏｆ ＰＶ ｂｕｓｅｓ ｉｎ ｔｈｅ ｓｙｓｔｅｍ ｗａｓ ｍｏｒｅ ｏｒ ｌｅｓｓ，ｔｈｅ ｃａｌｃｕｌａｔｉｎｇ ｓｐｅｅｄ ｗｉｔｈ Ｎｅｗｔｏｎ ｍ ｅｔｈｏｄ ｉｎ ｒｅｃｔａｎｇｕｌａｒ ｃｏｏｒｄｉｎａｔｅｓ ｗａｓ ａｌｗａｙｓ ｆａｓｔｅｒ ｔｈａｎ ｔｈｅ ｓｐｅｅｄ ｉｎ ｐｏｌａｒ ｃｏｏｒｄｉｎａｔｅｓ．Ｔｈｅ ｒｅｓｕｌｔｓ ｏｆ ｐｒｏｇｒａｍｍｉｎｇ ａｎｄ ｃａｌｃｕｌａｔｉｎｇ ｔｈｅ ｓｙｓｔｅｍｓ ｏｆ ＩＥＥＥ 一３０、一５７、

function PowerFlowNRRec(linedatabusdataerror)%节点标号无特殊要求clcif nargin < 1% linedata = [% 1 2 0.1 0.4 0.01528 1;% 1 3 0 0.3 0 1.1;% 1 4 0.12 0.5 0.01920 1;% 2 4 0.08 0.40 0.01413 1];% busdata = [% 1 0 -0.30 -0.18 1 0;% 2 0 -0.55 -0.13 1 0;% 3 2 0.5 0 1.10 0;% 4 1 0 0 1.05 0;];linedata = [1 4 0.1 0.4 0.01528 1;1 3 0 0.3 0 1.1;1 2 0.12 0.5 0.01920 1;4 2 0.08 0.40 0.01413 1];busdata = [1 0 -0.30 -0.18 1 0;4 0 -0.55 -0.13 1 0;3 2 0.5 0 1.10 0;2 1 0 0 1.05 0;];endif nargin < 3error = 10^-5;endclc%优化节点排序pqbus = find(busdata(:2)==0);pvbus = find(busdata(:2)==2);baxxxxsebus = find(busdata(:2)==1);bussort = [pqbus' pvbus' baxxxxsebus'];Num = busdata(bussort1);%对节点进行重新排列，使其排列顺序为 PQ节点->PV节点->平衡节点linedata = NumResort(linedataNum);P = busdata(:3);Q = busdata(:4);U = busdata(:5);Us = U;deta = busdata(:6);e = U.*cos(deta);f = U.*sin(deta);busnum = length(busdata(bussort1));pqnum = length(pqbus);%生成节点导纳矩阵Y = BuildY(linedata);for k = 1:100[dPdQdU2] = CalcCollection(PQUsefYbusnumpqnum); if max(abs([dP dQ]))<errorbreak;endJ = CalcJ(efYbusnumpqnum)dPQU = zeros(busnum-11);for ii = 1:pqnumdPQU(2*ii-1) = dP(ii);dPQU(2*ii) = dQ(ii);endfor ii = pqnum+1:busnum-1dPQU(2*ii-1) = dP(ii);dPQU(2*ii) = dU2(ii);enddef = Jordan(JdPQU);for ii = 1:busnum-1f(ii) = f(ii)+def(2*ii-1);e(ii) = e(ii)+def(2*ii);endendif k >= 100disp('不收敛!');return;end[SbusSijSvc] = CalcS(linedataefYbusnum);disp('编号:');disp(Num');disp('节点电压:');disp(e'+1j*f');disp('平衡节点功率:');disp(Svc);disp('线路功率');for n = 1:length(Sij(:1))fprintf('%o %o 'Num(Sij(n1))Num(Sij(n2)))disp(Sij(n3))enddisp('节点功率');for n = 1:busnumfprintf('%o 'Num(n))disp(Sbus(n))endendfunction [SbusSijSvc] = CalcS(linedataefYbusnum)Sbus = zeros(1busnum);Qi = zeros(1busnum);Ud = e+1j*f;%公式11-57for ii = 1:busnumfor jj = 1:busnumSbus(ii) = Sbus(ii)+Ud(ii)*conj(Y(iijj))*conj(Ud(jj));endendSvc = Sbus(busnum);Sij = zeros(length(linedata(:1))3);for n=1:length(linedata(:1))ii=linedata(n1);jj=linedata(n2);Sij(2*n-11) = ii;Sij(2*n-12) = jj;Sij(2*n-13) = Ud(ii)*(conj(Ud(ii))...*conj(1j*linedata(n5)+linedata(n6)*(linedata(n6)-1)/(linedata(n3)+1j*linedata(n4)))... +(conj(Ud(ii))-conj(Ud(jj)))*conj(-Y(iijj)));ii=linedata(n2);jj=linedata(n1);Sij(2*n1) = ii;Sij(2*n2) = jj;Sij(2*n3) = Ud(ii)*(conj(Ud(ii))...*conj(1j*linedata(n5)+(1-linedata(n6))/(linedata(n3)+1j*linedata(n4)))...+(conj(Ud(ii))-conj(Ud(jj)))*conj(-Y(iijj)));endendfunction linedata1 = NumResort(linedataNum) %对节点进行重排linedata1 = linedata;for n = 1:length(Num)n1 = find(linedata(:1)==Num(n));n2 = find(linedata(:2)==Num(n));linedata1(n11) = n*ones(length(n1)1);linedata1(n22) = n*ones(length(n2)1);endendfunction J = CalcJ(efYbusnumpqnum)%计算雅可比矩阵G = real(Y);B = imag(Y);J = zeros(2*(busnum-1));H = zeros(pqnumpqnum);N = H;M = H;L = H;for ii = 1:pqnumfor jj = 1:pqnumif ii~=jjH(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii); N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii); M(iijj) = -N(iijj);L(iijj) = H(iijj);J(2*ii-12*jj-1) = H(iijj);J(2*ii-12*jj) = N(iijj);J(2*ii2*jj-1) = M(iijj);J(2*ii2*jj) = L(iijj);endendIii = G(iiii)*e(ii)-B(iiii)*f(ii)...+1j*(G(iiii)*f(ii)+B(iiii)*e(ii));for jj = 1:busnumif ii~=jjIii = Iii+G(iijj)*e(jj)-B(iijj)*f(jj)... +1j*(G(iijj)*f(jj)+B(iijj)*e(jj));endendH(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)+imag(Iii); N(iiii) = G(iiii)*e(ii)+B(iiii)*f(ii)+real(Iii); M(iiii) = -G(iiii)*e(ii)-B(iiii)*f(ii)+real(Iii); L(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)-imag(Iii); J(2*ii-12*ii-1) = H(iiii);J(2*ii-12*ii) = N(iiii);J(2*ii2*ii-1) = M(iiii);J(2*ii2*ii) = L(iiii);endH = zeros(pqnumbusnum-pqnum-2);N = H;M = H;L = H;for ii = 1:pqnumfor jj = pqnum+1:busnum-1if ii~=jjH(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);M(iijj) = -N(iijj);L(iijj) = H(iijj);J(2*ii-12*jj-1) = H(iijj);J(2*ii-12*jj) = N(iijj);J(2*ii2*jj-1) = M(iijj);J(2*ii2*jj) = L(iijj);endendendH = zeros(busnum-pqnum-2pqnum);N = H;R = H;S = H;for ii = pqnum+1:busnum-1for jj = 1:pqnumif ii~=jjH(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);R(iijj) = 0;S(iijj) = 0;J(2*ii-12*jj-1) = H(iijj);J(2*ii-12*jj) = N(iijj);J(2*ii2*jj-1) = R(iijj);J(2*ii2*jj) = S(iijj);endendendH = zeros(busnum-pqnum-2);N = H;R = H;S = H;for ii = pqnum+1:busnum-1for jj = pqnum+1:busnum-1if ii~=jjH(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);R(iijj) = 0;S(iijj) = 0;J(2*ii-12*jj-1) = H(iijj);J(2*ii-12*jj) = N(iijj);J(2*ii2*jj-1) = R(iijj);J(2*ii2*jj) = S(iijj);endendIii = G(iiii)*e(ii)-B(iiii)*f(ii)...+1j*(G(iiii)*f(ii)+B(iiii)*e(ii));for jj = 1:busnumif ii~=jjIii = Iii+G(iijj)*e(jj)-B(iijj)*f(jj)...+1j*(G(iijj)*f(jj)+B(iijj)*e(jj));endendH(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)+imag(Iii);N(iiii) = G(iiii)*e(ii)+B(iiii)*f(ii)+real(Iii);R(iiii) = 2*f(ii);S(iiii) = 2*e(ii);J(2*ii-12*ii-1) = H(iiii);J(2*ii-12*ii) = N(iiii);J(2*ii2*ii-1) = R(iiii);J(2*ii2*ii) = S(iiii);endendfunction [dPdQdU2] = CalcCollection(PQUsefYbusnumpqnum) %计算修正值G = real(Y);B = imag(Y);Pt = P;Qt = Q;dP = zeros(1busnum);dQ = zeros(1busnum);dU2 = zeros(1busnum);for ii = 1:pqnumPi = 0;Qi = 0;for jj = 1:busnumPi = Pi+e(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))... +f(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));Qi = Qi+f(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))... -e(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));enddP(ii) = P(ii)-Pi;dQ(ii) = Q(ii)-Qi;Pt(ii) = Pi;Qt(ii) = Qi;endfor ii = pqnum+1:busnum-1Pi = 0;for jj = 1:busnumPi = Pi+e(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))... +f(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));enddP(ii) = P(ii)-Pi;Pt(ii) = Pi;dU2(ii) = Us(ii)^2-(e(ii)^2+f(ii)^2);endendfunction Y = BuildY(linedata)if nargin < 1linedata = [1 4 0.1 0.4 0.01528 1;1 3 0 0.3 0 1.1;1 2 0.12 0.5 0.01920 1;4 2 0.08 0.40 0.01413 1];endnf = linedata(:1);nt = linedata(:2);r = linedata(:3);x = linedata(:4);b = linedata(:5);k = linedata(:6);branchnum = length(nf);busnum = max([nf'nt']);y = ones(branchnum1)./(r+1j*x);Y = zeros(busnum);for n = 1:branchnumY(nf(n)nt(n)) = Y(nf(n)nt(n))-k(n)*y(n);Y(nt(n)nf(n)) = Y(nf(n)nt(n));Y(nf(n)nf(n)) = Y(nf(n)nf(n))+k(n)^2*y(n)+1j*b(n);Y(nt(n)nt(n)) = Y(nt(n)nt(n))+y(n)+1j*b(n);endendfunction s = Jordan(Ab)%约当消元if nargin < 1A=[1 2 3;2 7 5;1 4 9];b=[1 6 -3]';endif rank(A)~=rank([Ab])error('A矩阵的秩和增广矩阵的秩不相同方程不存在唯一解'); end[~n]=size(A);A(:n+1) = b;for k = 1:nA(kk:end) = A(kk:end)/A(kk);% 规格化for q = 1:k-1% 消上三角if A(qk)~=0A(qk:end) = A(qk:end)-A(kk:end).*A(qk);endendfor p = k+1:n% 消下三角if A(pk)~=0A(pk:end) = A(pk:end)-A(kk:end).*A(pk);endendends = A(:end);end。

第 33 卷 第9期 2009 年 5 月 10 日
Vol. 33 No . 9 May 10 , 2009
带最优乘子牛顿法在交直流系统潮流计算中的应用
胡泽春 , 严 正
( 上海交通大学电子信息与电气工程学院 , 上海市 200240)
摘要 : 交直流系统潮流计算的统一迭代法的收敛性较交替迭代法更优 , 但对一些病态运行条件和 潮流无解的情况 ,统一迭代法将发散 ,得不到有用信息 。文中将交流系统中广泛应用的最优乘子潮 流方法推广到交直流系统的潮流计算中 。根据直流系统稳态方程的特点引入辅助变量 ,使得直流 系统潮流方程全部转换为二次或线性方程 。给出了直角坐标系和极坐标系下交直流系统带最优乘 子牛顿法潮流的实现方法 。通过算例验证了所述方法的有效性 ,并对 2 种坐标系下的结果进行了 比较和分析 。 关键词 : 潮流计算 ; 牛顿法 ; 最优乘子 ; 直角坐标和极坐标系 ; 交直流系统 中图分类号 : TM744
1 ( i) μ = arg min μ 2 式中 : n 为潮流方程数 。 定义 :
( i) ai = f i ( x )
n i =1
PdR = V dR Id PdI = V dI Id V dR = V dI + Rd Id
式中 : V d , Id 分别为直流电压和电流 ; V t 为交流侧电 压 ; a 为换流变压器变比 ; nb 为换流桥个数 ;α,γ分 别为触发角和熄弧角 ; Pd , Qd 分别为换流器吸收的 有功和无功功率 ( 包括变压器) ; S d 为换流器和换流 变压器吸收的复功率 ; X c 为换相电抗 ; Rd 为直流电 阻 ; k≈0. 995 ; 下标 R ,I 分别表示整流侧和逆变侧 。 以上 9 个方程中共有 13 个待求变量 。根据整 流站和逆变站的控制方式可确定 4 个方程 , 这样直 流系统的方程数正好等于待求变量数 。 1. 2 方程变换 在计算时 ,将 co s α R 和 co s γ I 当做待求量 。对 不同的直流系统控制方式 , 式 ( 1 ) ～式 ( 4 ) 可能是三 次、 二次或一次方程 。为了易于实现最优乘子潮流 计算方法 ,可对直流系统方程进行变换 。首先引入 2 个辅助变量 V sR 和 V sI ,定义为 :

Pgi min Pgi Pgi max
Qgi min Qgi Qgi max
ijmin ij ijmax
21
三、极坐标下的牛拉法潮流 计算
n m nm1 1 m
PQ
n
Pi Vi Vj (Gij cosij Bij sin ij ) j 1
n
Qi Vi Vj (Gij sinij Bij cosij ) j 1
⑦ 引入修正系数； ⑧ 初值、平值电压启动。
19
计算步骤
输入原始数据
形成节点导纳矩阵
给定节点电压初值
e(0) i
,
f (0) i
k 0
用公式计算 Pi(k ) , Qi(k )及Vi2(k )
是
max{| Pi(k) , Qi(k) , Vi2(k) |} ?
否
按公式计算雅可比矩阵各元素
解修正方程式，求ei(
已知
x(0) 1
,
x(0) 2
，与真解的差为
x1(0) , x2(0)
f1 (
x(0) 1
x1(0)
,
x(0) 2
x2(0)
)
0
f2(
x(0) 2
x2(0) ,
x(0) 2
x2(0)
)
0
8
一、牛顿一拉夫逊法的基本原理
展开：
f1
(
x(0) 1
,
x(0) 2
)
f1 x1
x1( 0 )
0
f1 x2
0
x2( 0 )
10.954526
x4
x3
f ( x3 ) f ( x3 )
10.954526 0.00163988 2 10.954526

基于MATLAB的直角坐标下牛顿拉夫逊法潮流计算基于MATLAB的直角坐标下牛顿-拉夫逊法潮流计算摘要潮流计算，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。

潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。

通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。

待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。

它是基于配电网络特有的层次结构特性，论文提出了一种新颖的分层前推回代算法。

该算法将网络支路按层次进行分类，并分层并行计算各层次的支路功率损耗和电压损耗，因而可大幅度提高配电网潮流的计算速度。

论文在MATLAB环境下，利用其快速的复数矩阵运算功能，实现了文中所提的分层前推回代算法，并取得了非常明显的速度效益。

另外，论文还讨论发现，当变压器支路阻抗过小时，利用Π型模型会产生数值巨大的对地导纳，由此会导致潮流不收敛。

为此，论文根据理想变压器对功率和电压的变换原理，提出了一种有效的电压变换模型来处理变压器支路，从而改善了潮流算法的收敛特性。

关键词：电力系统；潮流分析；MATLABAbstractFlow calculation is an important analysis function of power system and is the necessary facility of fault analysis, relay protection setting and security analysis. In addition, the traditional design method is a structured program design method based on functional decomposition, the entire software engineering as a combination of objects, as the domain of a particular issue, the composition of the object will remain basically unchanged Therefore, this decomposition methodbased on object design software structure relatively stable, easy to maintain and expand. . Combine the characteristics of power systems, software running on the use of MATLAB language WINDOWS OS graphical flow calculation software. The main features of the system are simple and intuitive graphical interface and stable operation. Calculated accurately Calculations, the algorithm has done a number of improvements to enhance the computing speed, the various types of effective package makes the procedure has good modularity maintainability and reusability. The MATLAB language is used to calculate flow distribution of power system in this paper. The typical examples explain that the method has the characteristics of simple programming high calculation efficiency and matching people habit the calculation result can satisfy the engineering calculation needs and at the same time verify the usefulness of the method.Key words: Electric power system; flow calculation; MATLAB 第一章电力系统潮流计算概述1.1电力系统潮流概述潮流计算是电力系统分析中的一种最基本的计算，它的任务是在给定的接线方式和运行条件下，确定系统的运行状态，如各母线上的电压（幅值和相角）、网络中的功率分布及功率损耗等，是电力系统的稳态计算。