电力系统分析计算机算法PSDBPA潮流计算实验报告

第一篇：电力系统分析计算机算法PSDBPA潮流计算实验报告

电力系统分析的计算机算法

实验报告

学生姓名 课 程 电力系统分析的计算机算法 学 号

专 业 电气工程及其自动化 指导教师 邱晓燕

二Ο一四 年 六 月 二日

实验一

潮流计算

一、实验目的

1．了解并掌握电力系统计算机算法的相关原理。

2．了解和掌握PSD-BPA电力系统分析程序稳态分析方法（即潮流计算）。3．了解并掌握PSD-BPA电力系统分析程序单线图和地理接线图的使用。

二、实验背景

随着科学技术的飞速发展，电力系统也在不断地发展，电网通过互联变得越来越复杂，同时也使系统稳定问题越来越突出。无论是电力系统规划、设计还是运行，对其安全稳定进行分析都是极其重要的。

PSD-BPA软件包主要由潮流和暂稳程序构成，具有计算规模大、计算速度快、数值稳定性好、功能强等特点，已在我国电力系统规划、调度、生产运行及科研部门得到了广泛应用。

本实验课程基于PSD-BPA平台，结合《电力系统分析计算机算法》课程，旨在引导学生将理论知识和实际工程相结合，掌握电力系统稳态、暂态分析的原理、分析步骤以及结论分析。清晰认知电力系统分析的意义。

三、原理和说明

1.程序算法

PSD-BPA电力系统分析程序稳态分析主要是潮流计算，软件中潮流程序的计算方法有P_Q分解法，牛顿_拉夫逊法，改进的牛顿－拉夫逊算法。采用什么算法以及迭代的最大步数可以由用户指定。

注：采用P-Q分解法和牛顿－拉夫逊法相结合，以提高潮流计算的收敛性能，程序通常先采用P-Q分解法进行初始迭代，然后再转入牛顿－拉夫逊法求解潮流。

2.程序主要功能

可进行交流系统潮流计算，也可进行包括双端和多端直流系统的交直流混合潮流计算。除了潮流计算功能外，该软件还具有自动电压控制、联络线功率控制、系统事故分析（N-1开断模拟）、网络等值、灵敏度分析、节点P－V、Q－V和P－Q曲线、确定系统极限输送水平、负荷静特性模型、灵活多样的分析报告、详细的检错功能等功能。

3.输入、输出相关文件 *.dat

潮流计算数据文件

*.bse

潮流计算二进制结果文件（可用于潮流计算的输入或稳定计算）*.pfo

潮流计算结果文件

*.map 供单线图格式潮流图及地理接线图格式潮流图程序使用的二进制结果文件

*.pff，*.pfd 中间文件（正常计算结束后将自动删除。不正常时，将留在硬盘上，可随时删除）

pwrflo.dis 储存一个潮流作业计算时屏幕显示的信息。pfcard.def 定义潮流程序卡片格式文件，用户可更改及调整该文件。该文件安装时放在与潮流程序相同的目录中。打开TextEdit应用程序时先读入该文件。4.程序常用控制语句

常用的控制语句主要包括：

(1)指定潮流文件开始的一级控制语句“(POWERFLOW, CASEID=方式名, PROJECT=工程名)”

(2)指定计算方法和最大迭代次数的控制语句“/SOL_ITER, DECOUPLED=PQ法次数, NEWTON=牛拉法次数”；

(3)指定计算结果输出的控制语句“/P_OUTPUT_LIST, „”；(4)指定计算结果输出顺序的控制语句“/RPT_SORT= „”；

(5)指定计算结果分析列表的控制语句“/P_ANALYSIS, LEVEL= ？”；(6)指定潮流结果二进制文件名的控制语句“/NEW_BASE, FILE = 文件名”；

(7)指定潮流图和地理接线图使用的结果文件控制语句“/PF_MAP，FILE=文件名”；

(8)指定网络数据的控制语句“/NETWORK_DATA”；(9)指定潮流数据文件结束的控制语句“(END)”； 5.计算结果介绍（PFO文件）

潮流计算结果文件内容主要分下述几个方面： 1）程序控制语句列表。

2）输入、输出文件及输出的内容列表。

3）错误信息。如为致命性错误，则中断计算。4）误差控制参数列表。5）迭代过程。6）计算结果输出：

详细计算结果列表：按节点、与该节点相联接支路顺序，并根据用户的要求（通过控制语句控制）可按照字母、分区或区域排序输出潮流计算结果。

分析报告列表：并根据用户的要求（通过控制语句控制），输出各种潮流分析报告。

7）错误信息统计。6.算例

IEEE 9节点例题：

图1 IEEE9节点系统接线图

节点参数、线路参数及变压器参数分别见表1～表3。

表1 IEEE 9节点算例节点参数

表2 IEEE 9节点算例线路参数

表3 IEEE 9节点算例变压器参数

注：表1-表3中功率基准值为100MVA；电阻、电感值为标幺值。对应于上述系统及数据的潮流计算数据（IEEE90.DAT）见例1。例1：

(POWERFLOW,CASEID=IEEE9,PROJECT=IEEE_9BUS_TEST_SYSTEM)/SOL_ITER,DECOUPLED=2,NEWTON=15,OPITM=0./P_INPUT_LIST,ZONES=ALL /P_OUTPUT_LIST,ZONES=ALL /RPT_SORT=ZONE /NEW_BASE,FILE=IEEE90.BSE /PF_MAP,FILE = IEEE90.MAP /NETWORK_DATA BS GEN1

16.501 999.999.1.04 B

GEN1

230.01

B

STATIONA 230.01 125.50.0 0.B

STATIONB 230.01 90.30.0 0.B

STATIONC 230.01 100.35.0 0.000 B

GEN2

230.01

BE GEN2

18.001 163.999 10 25 B

GEN3

230.01 BE GEN3

13.801 85.999.1025

.L-----------------transmission lines----------------------------L

GEN1 230.STATIONA230..0100.0850.0440 L

GEN1 230.STATIONA230.2.0100.0850.0440 L

GEN1230.STATIONB230..0170.0920.0395 L

STATIONA230.GEN2230..0320.1610.0765 L

STATIONB230.GEN3230..0390.1700.0895 L

GEN2230.STATIONC230..0085.0720.03725 L

STATIONC230.GEN3230..0119.1008.05225.T-----transformers---------

T

GEN116.5 GEN1230..0576 16.5 230.T

GEN218.0 GEN2230..0625 18.0 230.T

GEN313.8 GEN3230..0586 13.8 230.(END)

四、实验过程及结果

（一）IEEE9节点算例： 1.系统接线图：

2.在BPA软件建立模型，并进行计算，结果如下： 1）系统数据 2）计算过程迭代信息及详细的输出列表：

小结

3.406

-60.2

0.000

28.2

0.000

0.0

3.406

-32.0

--------------

--------------

--------------

--------------

总结

3.406

-60.2

0.000

28.2

0.000

0.0

3.406

-32.0 * 并联无功补偿数据列表

/----------电容器（Mvar）-----------/

/-----------电抗器（Mvar）-------------/

区域/分区

最大容量

使用容量

备用

未安排容量

最大容量

使用容量

备用

未安排容量

01

73.4

73.4

0.0

0.0

0.0

0.0

0.0

0.0

-------

-------

-------

-------

-------

-------

-------

-------

总结

73.4

73.4

0.0

0.0

0.0

0.0

0.0

0.0

TRANSMISSION LINES CONTAINING COMPENSATION

OWN ZONE BUS1

BASE1 ZONE BUS2

BASE2

ID PERCENT

CASE CONTAINS NO TRANSMISSION LINES WITH SERIES COMPENSATION

* 节点相关数据列表

节点

电压

/--------发电--------/ /---负荷----/

/-----无功补偿-----/ 类型 拥有者 分区

电压/角度

kV

MW

MVAR 功率因数

MW

MVAR

使用的存在的未安排

PU/度

发电机1

16.5

16.5

105.4

23.1 0.98

0.0

0.0

0.0

0.0

0.0

S

01

1.000/

0.0

发电机2

18.0

18.0

180.0

40.6 0.98

17.0

8.0

0.0

0.0

0.0

E

01

1.000/

5.4

发电机3

13.8

13.8

85.0

13.8 0.99

0.0

0.0

0.0

0.0

0.0

E

01

1.000/

1.6

母线1

230.0

239.3

0.0

0.0

0.0

0.0

21.6

21.6

0.0

01

1.040/-3.5

母线2

230.0

238.3

0.0

0.0

35.0

10.0

0.0

0.0

0.0

01

1.036/-0.6

母线3

230.0

240.3

0.0

0.0

0.0

0.0

0.0

0.0

0.0

01

1.045/-1.3

母线A

230.0

232.6

0.0

0.0

125.0

70.0

20.5

20.5

0.0

01

1.011/-6.0

母线B

230.0

234.1

0.0

0.0

90.0

40.0

10.4

10.4

0.0

01

1.018/-5.7

母线C

230.0

235.6

0.0

0.0

100.0

55.0

21.0

21.0

0.0

01

1.024/-3.1

--------------

--------------------------------

整个系统

370.4

77.6

367.0

183.0

73.4

73.4

0.0

电容器总和

73.4

73.4

0.0

电抗器总和

0.0

0.0

0.0 * 旋转备用数据列表

------------有功功率-----------

------------------------无功功率-----------------------

区域/分区

最大值

实际出力

备用

最大值

最小值

已发无功

吸收无功

备用

(MW)

(MW)

(MW)

(MVAR)

(MVAR)

(MVAR)

(MVAR)

(MVAR)

01

370.4

370.4

0.0

2997.0

0.0

77.6

0.0

2919.4

-------

-------

------

-------

-------

-------

------

-------

总结

370.4

370.4

0.0

2997.0

0.0

77.6

0.0

2919.4

说明：

1.有功旋转备用不包含所有同步电动机的功率（如 抽水蓄能电机）。

有功出力为负值的发电机（包括电动机）作为负荷处理，不统计在内。

当最大出力值小于实际出力时，统计时最大出力值用实际出力值代替。

2.无功旋转备用不包含同步调相机的无功功率。

无功旋转备用只统计有功出力大于0并且基准电压小于30kV的发电机。

* 潮流计算迭代过程和平衡节点相关信息数据

计算结果收敛。牛顿-拉夫逊法迭代次数为 5次。

各区域平衡机出力数据列表

区域

平衡机

电压

额定有功

有功出力

无功出力

有功负荷

无功负荷

所属分区

SYSTEM

发电机1 16.5

1.000

0.00

105.41

23.11

0.00

0.00

01

* 没有遇到错误信息 23:03:48 3）单线图：

（二）课本习题：E2-5 1.网络接线图：

2.程序：

(POWERFLOW,CASEID=IEEE9,PROJECT=IEEE_9BUS_TEST_SYSTEM)/SOL_ITER,DECOUPLED=2,NEWTON=15,OPITM=0 /P_OUTPUT_LIST,ZONES=ALL /RPT_SORT=ZONE /NEW_BASE,FILE=IEEE90.BSE /PF_MAP,FILE = IEEE90.MAP /NETWORK_DATA.BUS-----------------节点数据-----BS

母线4

999

999

1.050

B

母线1

0.32 0.20

B

母线2

0.56 0.16

BE

母线3

0.5 999

1.10

.L-----------------支路数据-----L

母线1

母线2

0.11 0.40

0.015

L

母线2

母线4

0.08 0.40

0.014

L

母线4

母线1

0.12 0.51

0.019

.T--------------变压器数据，包括普通变压器、移相器、带调节的变压器等。

T

母线1

母线3

0.07 0.35

(END)

3.计算结果

4.系统单线图

五、总结及思考题

实验中遇到的问题及解决方法：

路径错误——————重设各个参数路径 卡片无法识别—————将参数规范化

本次实验使我初步掌握了PSD-BPA软件在电力系统潮流计算中的使用方法，收获良多，为今后的工作打下了基础。获益匪浅。

电力系统潮流计算是研究电力系统稳态运行情况的一种基本计算。它的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。电力系统潮流计算的结果是电力系统稳定计算和故障分析的基础。

潮流计算的方法有最基本的手算迭代方法，而利用电子计算机进行潮流计算从20世纪50年代中期就已经开始。此后，潮流计算曾采用了各种不同的方法，这些方法的发展主要是围绕着对潮流计算的一些基本要求进行的。从数学上说，潮流计算是求解一组由潮流方程描述的非线性代数方程组。牛顿-拉夫逊方法是解非线性代数方程组的一种基本方法，在潮流计算中也得到应用。PSD-BPA仿真软件中潮流计算模型建模的注意事项?

第二篇：电力系统分析潮流计算的计算机算法

潮流计算的计算机算法实验报告

姓名：学号：班级：

一、实验目的

掌握潮流计算的计算机算法。

熟悉MATLAB，并掌握MATLAB程序的基本调试方法。

二、实验准备

根据课程内容，熟悉MATLAB软件的使用方法，自行学习MATLAB程序的基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法(或PQ分解法)的计算程序，用相应的算例在MATLAB上进行计算、调试和验证。

三、实验要求

每人一组，在实验课时内，用MATLAB调试和修改运行程序，用算例计算输出潮流结果。

四、实验程序

clear;

%清空内存

n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');Y=zeros(n);

Times=1;

%一：创建节点导纳矩阵

for i=1:n1

if B1(i,6)==0

%不含变压器的支路

p=B1(i,1);

q=B1(i,2);

Y(p,q)=Y(p,q)-1/B1(i,3);

Y(q,p)=Y(p,q);

Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);

Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);

else

%含有变压器的支路

p=B1(i,1);

q=B1(i,2);

Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));

Y(q,p)=Y(p,q);

Y(p,p)=Y(p,p)+1/B1(i,3);

Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));

end end Y;

%将OrgS、DetaS初始化

OrgS=zeros(2*n-2,1);

DetaS=zeros(2*n-2,1);

%二：创建OrgS，用于存储初始功率参数

h=0;j=0;for i=1:n

%对PQ节点的处理

if i~=isb&B2(i,6)==2

h=h+1;

for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*Imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

end

end end %三：对PV节点的处理，注意这时不可再将h初始化为0 for i=1:n

if i~=isb&B2(i,6)==3

h=h+1;

for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

end

end end OrgS;%四：创建PVU 用于存储PV节点的初始电压

PVU=zeros(n-h-1,1);t=0;for i=1:n

if B2(i,6)==3

t=t+1;

PVU(t,1)=B2(i,3);

end end PVU;%五：创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量

h=0;for i=1:n

%对PQ节点的处理

if i~=isb&B2(i,6)==2

h=h+1;

DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);

DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);

end end t=0;for i=1:n

%六：对PV节点的处理，注意这时不可再将h初始化为0

if i~=isb&B2(i,6)==3

h=h+1;

t=t+1;

DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);

DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;

end end DetaS;%七：创建I，用于存储节点电流参数

i=zeros(n-1,1);h=0;for i=1:n

if i~=isb

h=h+1;

I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));

end end I;%八：创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n

%对PQ节点的处理

if B2(i,6)==2

h=h+1;

for j=1:n

if j~=isb

k=k+1;

if i==j

%对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));

else

%非对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);

end

if k==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0;

end

end

end

end end k=0;for i=1:n

%对PV节点的处理

if B2(i,6)==3

h=h+1;

for j=1:n

if j~=isb

k=k+1;

if i==j

%对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));

Jacbi(2*h,2*k-1)=2*imag(B2(i,3));

Jacbi(2*h,2*k)=2*real(B2(i,3));

else

%非对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));

Jacbi(2*h,2*k-1)=0;

Jacbi(2*h,2*k)=0;

end

if k==(n-1)

%将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0;

end

end

end

end end Jacbi;%九：求解修正方程，获取节点电压的不平衡量

DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正节点电压 j=0;for i=1:n

%对PQ节点处理

if B2(i,6)==2

j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);

end end for i=1:n

%对PV节点的处理

if B2(i,6)==3

j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);

end end B2;%十：开始循环********************************************************************** while abs(max(DetaU))>pr OrgS=zeros(2*n-2,1);

%！!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵

h=0;j=0;for i=1:n

if i~=isb&B2(i,6)==2

h=h+1;

for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

end

end end for i=1:n

if i~=isb&B2(i,6)==3

h=h+1;

for j=1:n

OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

end

end end OrgS;%创建DetaS h=0;for i=1:n

if i~=isb&B2(i,6)==2

h=h+1;

DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);

DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);

end end t=0;for i=1:n

if i~=isb&B2(i,6)==3

h=h+1;

t=t+1;

DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);

DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;

end end DetaS;%创建I i=zeros(n-1,1);h=0;for i=1:n

if i~=isb

h=h+1;

I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));

end end I;%创建Jacbi Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n

if B2(i,6)==2

h=h+1;

for j=1:n

if j~=isb

k=k+1;

if i==j

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));

else

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));

Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);

end

if k==(n-1)

k=0;

end

end

end

end end k=0;for i=1:n

if B2(i,6)==3

h=h+1;

for j=1:n

if j~=isb

k=k+1;

if i==j

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));

Jacbi(2*h,2*k-1)=2*imag(B2(i,3));

Jacbi(2*h,2*k)=2*real(B2(i,3));

else

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));

Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));

Jacbi(2*h,2*k-1)=0;

Jacbi(2*h,2*k)=0;

end

if k==(n-1)

k=0;

end

end

end

end end Jacbi;DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正节点电压 j=0;for i=1:n

if B2(i,6)==2

j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);

end end for i=1:n

if B2(i,6)==3

j=j+1;

B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);

end end B2;Times=Times+1;

%迭代次数加1 end Times;

五、实验流程

六、实验结果

参数输入：

运行结果：

七、实验体会

通过这次实验，让我第一次接触到了MATLAB,并深切体会到了它的强大之处；潮流计算的计算机算法的实现不仅巩固了我的学过的知识，还让我学到一些MATLAB的编程，虽然在实验的过程中出现了很多的错误，但在老师的细心指导下，问题都解决啦；计算机为我们省去了大量的人工计算，希望在以后的学习中能接触到更多的软件，学习到更多的知识。

第三篇：电力系统分析潮流实验报告

南昌大学实验报告

电力系统潮流计算实验 学生姓名： 学 号： 专业班级： 实验类型：□ 验证 □ 综合 ■ 设计 □ 创新 实验日期： 实验成绩：

一、实验目的：

本实验通过对电力系统潮流计算的计算机程序的编制与调试，获得对复杂电力系统进行潮流计算的计算机程序，使系统潮流计算能够由计算机自行完成，即根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，由计算程序运行完成该电力系统的潮流计算。通过实验教学加深学生对复杂电力系统潮流计算计算方法的理解，学会运用电力系统的数学模型，掌握潮流计算的过程及其特点，熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。

二、实验内容：

编制调试电力系统潮流计算的计算机程序。程序要求根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，完成该电力系统的潮流计算，要求计算出节点电压、功率等参数。

1、在各种潮流计算的算法中选择一种，按照计算方法编制程序。

2、将事先编制好的电力系统潮流计算的计算程序原代码由自备移动存储设备导入计算机。

3、在相应的编程环境下对程序进行组织调试。

4、应用计算例题验证程序的计算效果。

三、实验程序：

function [e,f,p,q]=flow_out(g,b,kind,e,f)%计算潮流后efpq的终值 s=flow(g,b,kind,e,f);k=0;

while max(abs(s))>10^-5 J=J_out(g,b,kind,e,f);J_ni=inv(J);dv=J_ni*s;l=length(dv)/2;

for i=1:l

e(i)=e(i)-dv(2*i-1);f(i)=f(i)-dv(2*i);

end

s=flow(g,b,kind,e,f);end

l=length(e);for i=1:l s1=0;s2=0;

for j=1:l

s1=s1+g(i,j)*e(j)-b(i,j)*f(j);s2=s2+g(i,j)*f(j)+b(i,j)*e(j);

end

p(i)=e(i)*s1+f(i)*s2;q(i)=f(i)*s1-e(i)*s2;end

function s=flow(g,b,kind,e,f)%计算当前ef与规定的pqv的差值 l=length(e);s=zeros(2*l-2,1);for i=1:(l-1)s1=0;s2=0;

for j=1:l

s1=s1+g(i,j)*e(j)-b(i,j)*f(j);s2=s2+g(i,j)*f(j)+b(i,j)*e(j);

end

s(2*i-1)=kind(2,i)-e(i)*s1-f(i)*s2;

if kind(1,i)==1

s(2*i)=kind(3,i)-f(i)*s1+e(i)*s2;

else

s(2*i)=kind(3,i)^2-f(i)^2-e(i)^2;

end end

function J=J_out(g,b,kind,e,f)%计算节点的雅克比矩阵 l=length(e);

J=zeros(2*l-2,2*l-2);for i=1:(l-1);

if kind(1,i)==1

s=PQ_out(g,b,e,f,i);

for j=1:(2*l-2)J(2*i-1,j)=s(1,j);J(2*i,j)=s(2,j);

end

else

s=PV_out(g,b,e,f,i);for j=1:(2*l-2)J(2*i-1,j)=s(1,j);J(2*i,j)=s(2,j);

end

end end

function pq=PQ_out(g,b,e,f,i)%计算pq节点的雅克比矩阵 l=length(e);pq=zeros(2,2*l-2);for j=1:(l-1)

if j==i s=0;

for k=1:l

s=s-(g(i,k)*e(k)-b(i,k)*f(k));

end

pq(1,2*i-1)=s-g(i,i)*e(i)-b(i,i)*f(i);s=0;

for k=1:l

s=s-(g(i,k)*f(k)+b(i,k)*e(k));

end

pq(1,2*i)=s+b(i,i)*e(i)-g(i,i)*f(i);s=0;

for k=1:l

s=s+(g(i,k)*f(k)+b(i,k)*e(k));

end

pq(2,2*i-1)=s+b(i,i)*e(i)-g(i,i)*f(i);s=0;

for k=1:l

s=s-(g(i,k)*e(k)-b(i,k)*f(k));

end

pq(2,2*i)=s+g(i,i)*e(i)+b(i,i)*f(i);

else

pq(1,2*j-1)=-(g(i,j)*e(i)+b(i,j)*f(i));pq(1,2*j)=b(i,j)*e(i)-g(i,j)*f(i);pq(2,2*j)=-pq(1,2*j-1);pq(2,2*j-1)=pq(1,2*j);

end end

function pv=PV_out(g,b,e,f,i)%计算pv节点的雅克比矩阵 l=length(e);pv=zeros(2,2*l-2);for j=1:(l-1)

if j==i s=0;

for k=1:l

s=s-(g(i,k)*e(k)-b(i,k)*f(k));

end

pv(1,2*i-1)=s-g(i,i)*e(i)-b(i,i)*f(i);s=0;for k=1:l

s=s-(g(i,k)*f(k)+b(i,k)*e(k));

end

pv(1,2*i)=s+b(i,i)*e(i)-g(i,i)*f(i);pv(2,2*i-1)=-2*e(i);pv(2,2*i)=-2*f(i);

else

pv(1,2*j-1)=-(g(i,j)*e(i)+b(i,j)*f(i));pv(1,2*j)=b(i,j)*e(i)-g(i,j)*f(i);

end end

%数据输入

g=[1.042093-0.588235 0-0.453858-0.588235 1.069005 0-0.480769 0 0 0 0

-0.453858-0.480769 0 0.9344627];

b=[-8.242876 2.352941 3.666667 1.891074 2.352941-4.727377 0 2.403846 3.666667 0-3.333333 0 1.891074 2.40385 0 4.26159];e=[1 1 1.1 1.05];f=[0 0 0 0];kind=[1 1 2 0-0.3-0.55 0.5 1.05-0.18-0.13 1.1 0];

[e,f,p,q]=flow_out(g,b,kind,e,f);e f

四、例题及运行结果

在上图所示的简单电力系统中，系统中节点1、2为PQ节点，节点3为PV节点，节点4为平衡节点，已给定 P1s+jQ1s=-0.30-j0.18 P2s+jQ2s=-0.55-j0.13 P3s=0.5 V3s=1.10 V4s=1.05∠0° 容许误差ε=10-5

节点导纳矩阵：

各节点电压：

节点

e

f

v

ζ

1.0.984637-0.008596 0.984675-0.500172 2.0.958690-0.108387 0.964798-6.450306 3.1.092415

0.128955 1.100000

6.732347 4.1.050000

0.000000 1.050000

0.000000

各节点功率：

节点

P

Q 1-0.300000-0.180000 2 –0.550000-0.130000

30.500000-0.551305

40.367883

0.264698 结果：

五、思考讨论题

1.潮流计算有几种方法？简述各种算法的优缺点。

答：高斯迭代法（高斯塞德尔法），牛顿拉夫逊法以及P-Q分解法。高斯迭代法是直接迭代，对初值要求比较低，程序简单，内存小，但收敛性差，速度慢，多用于配电网或辐射式网络中；牛顿拉夫逊法是将非线性方程线性化之后再迭代的，对初值要求比较高，收敛性好，速度快，迭代次数少，运行时间短，被广泛使用；P-Q分解法是在极坐标牛顿法的基础上进行三个简化所得，有功、无功分开迭代，迭代次数比牛顿多一倍但运算量小，整体速度更快，运行时间更短，多用于110KV以上的高压电网中

2.在潮流计算中，电力网络的节点分几类？各类节点的已知量和待求量是什么？

答: PQ节点：P、Q为已知量，V、为待求量；PV节点：给定P、V，求Q、；平衡节点：给定V、，求P、Q。

3.潮流计算中的雅可比矩阵在每次迭代时是一样的吗？为什么？

答:不一样,因为每次迭代的电压、有功、无功都是与前一次不同的新值，所以每次迭代过程中，雅可比矩阵都是变化的。

六、实验心得

这次实验是通过matlab编写出一个潮流计算的程序。我这用了牛顿法直角坐标系来编写程序的。通过编写这次程序可以更深一步的理解潮流计算的步骤，也明白了在潮流计算中要注意的一些细节。

第四篇：电力系统分析潮流计算例题

3．1 电网结构如图3—11所示，其额定电压为10KV。已知各节点的负荷功率及参数：

S2(0.3j0.2)MVA，S3(0.5j0.3)MVA，S4(0.2j0.15)MVA

Z12(1.2j2.4)，Z23(1.0j2.0)，Z24(1.5j3.0)

试求电压和功率分布。

解：（1）先假设各节点电压均为额定电压，求线路始端功率。

P32Q320.520.32S23(R23jX23)(1j2)0.0034j0.006822VN10P42Q420.220.152S24(R24jX24)(1.5j3)0.0009j0.001922VN10

则： S23S3S230.5034j0.3068

S24S4S240.2009j0.1519

'S23S24S21.0043j0.6587

S12又

''P12Q121.004320.65872S12(R12jX12)(1.2j2.4)22VN10

0.0173j0.0346'故： S12S12S121.0216j0.6933

（2）再用已知的线路始端电压V110.5kV及上述求得的线路始端功率S

12，求出线路各点电压。

(P12R12Q12X12)1.02161.20.69332.4V120.2752kVV110.5 V2V1V1210.2248kV(P24R24Q24X24)V240.0740kVV4V2V2410.1508kVV2

(P23R23Q23X23)V230.1092kVV3V2V2310.1156kVV2

（3）根据上述求得的线路各点电压，重新计算各线路的功率损耗和线路始端功率。

0.520.32(1j2)0.0033j0.0066 S23210.120.220.152(1.5j3)0.0009j0.0018 S24210.15故 S23S3S230.5033j0.3066 S24S4S240.2009j0.1518

'S23S24S21.0042j0.6584 则 S12又

1.004220.65842S12(1.2j2.4)0.0166j0.0331 210.22 从而可得线路始端功率 S121.0208j0.6915 这个结果与第（1）步所得计算结果之差小于0.3%，所以第（2）和第（3）的结果可作为最终计算结果；若相差较大，则应返回第（2）步重新计算，直道相差较小为止。

3.2 如图所示简单系统，额定电压为110KV 双回输电线路，长度为80km，采用LGJ-150导线，其单位长度的参数为：r=0.21Ω/km，x=0.416Ω/km,b=2.74106S/km。变电所中装有两台三相110/11kV的变压器，每台的容量为15MVA,其参数为：

P040.5kW，Ps128kW,Vs%10.5,Io%3.5。母线A的实际运行电压为117kV，负荷功率：

SLDb30j12MVA,SLDc20j15MVA。当变压器取主轴时，求母线c的电压。

解（1）计算参数并作出等值电路。

输电线路的等值电阻、电抗和电纳分别为 RL1800.218.4 21 XL800.41616.6 Bc2802.74106S4.38104S

由于线路电压未知，可用线路额定电压计算线路产生的充电功率，并将其等分为两部分，便得

112QBBcVN4.381041102Mvar2.65Mvar22将QB分别接于节点A 和b，作为节点负荷的一部分。

两台变压器并联运行时，它们的等值电阻、电抗及励磁功率分别为

1PsVN211281102RT3.4 2221000SN21000151Vs%VN2110.51102RT42.4 22100SN2100153.515PojQo2(0.0405j)MVA0.08j1.05MVA100

变压器的励磁功率也作为接于节点b的负荷，于是节点b的负荷

SbSLDbjQB(P0jQ0)30j120.08j1.05j2.65MVA30.08j10.4MVA点c的功率即是负荷功率 Sc20j15MVA 这样就得到图所示的等值电路

节

（2）计算母线A输出的功率。

先按电力网络的额定电压计算电力网络中的功率损耗。变压器绕组中的功率损耗为

Sc202152RTjXTST(3.4j42.4)MVA2

110VN0.18j2.19MVA2由图可知

Sc'ScPTjQT20j150.18j2.19MVA20.18j17.19MVASc''Sc'Sb20.18j17.1930.08j10.4MVA50.26j27.59MVA

线路中的功率损耗为

S1SLVRLjXLN

50.26227.592(8.4j16.6)MVA2.28j4.51MVA2110''2 于是可得 S1'S1''SL50.26j27.592.28j4.51MVA52.54j32.1MVA

由母线A输出的功率为

SAS1'jQB52.54j32.1j2.65MVA52.54j29.45MVA

（3）计算各节点电压。

线路中电压降落的纵分量和横分量分别为

P1'RLQ1'XL52.248.432.116.6VLkV8.3kVVA117

P1'XLQ1'RL52.2416.632.18.4VLkV5.2kV

VA117b点电压为

Vb VAVLVL221178.325.22kV108.8kV变压器中电压降落的纵，横分量分别为

Pc'RTQc'XT20.183.417.1942.4VTkV7.3kV

Vb108.8

'Pc'XTQCRT20.1842.417.193.4VTkV7.3kV

Vb108.8归算到高压侧的c点电压

Vc' VbVTVT22108.87.327.32kV101.7kV变电所低压母线c的实际电压

1111VcVc101.7kV10.17kV

110110'如果在上述计算中都不计电压降落的横分量，所得结果为

Vb108.7kV，Vc'101.4kV，Vc10.14kV 与计及电压降落横分量的计算结果相比，误差很小。

3.3 某一额定电压为10kV的两端供电网，如图所示。线路L1、L2和L3导线型号均为LJ-185，线路长度分别为10km，4km和3km，线路L4为2km长的LJ-70导线；各负荷点负荷A10.50kV、VB10.40kV时的初始如图所示。试求V功率分布，且找到电压最低点。（线路参数LJ-185：z=0.17+j0.38Ω/km；LJ-70：z=0.45+j0.4Ω/km）

解 线路等值阻抗

ZL110(0.17j0.38)1.7j3.8

ZL24(0.17j0.38)0.68j1.52

ZL33(0.17j0.38)0.51j1.14 ZL42(0.45j0.4)0.9j0.8 求C点和D点的运算负荷，为 SCE SC0.320.162(0.9j0.8)1.04j0.925kVA 2102600j1600300j1601.04j0.9252901.04j1760.925kVA

SD600j2001600j10002200j1200kVA

循环功率

ScVAVBVN10.510.410339.430.17j0.38kVA580j129kVA 170.17j0.38Z12901.04722003j1760.9257j12003Sc17 1582.78j936.85580j1292162.78j1065.85kVASAC12901.0410220014j1760.92510j120014Sc 173518.26j2024.07580j1292938.26j1895.07kVASBDSACSBD2162.78j1065.852938.26j1895.075101.04j2960.92kVA SCSD2901.04j1760.9252200j12005101.04j2960.92kVA SCDSBDSD2938.26j1895.072200j1200738.26j695.07kVA

C点为功率分点，可推算出E点为电压最低点。进一步可求得E点电压 SAC2.1621.072(1.7j3.8)MVA98.78j220.8kVA

102' SAC2162.78j1065.8598.78j220.82261.56j1286.65kVA

VAC2.261.71.293.80.8328kV10.5

VCVAVAC10.50.83289.6672kV VCE0.3010.90.1610.80.041kV9.6672

VEVCVCE9.66720.0419.6262kV

3.4 图所示110kV闭式电网，A点为某发电厂的高压母线，其运行电压为117kV。网络各组件参数为：

变电所b SNXT63.5

20MVA，S00.05j0.6MVA，RT4.84，变电所c SNXT127

10MVA，S00.03j0.35MVA，RT11.4，负荷功率 SLDb24j18MVA，SLDc12j9MVA

试求电力网络的功率分布及最大电压损耗。

解（1）计算网络参数及制定等值电路。

线路Ⅰ： Z(0.27j0.423)6016.2j25.38

B2.6910660S1.61104S

2QB1.611041102Mvar1.95Mvar

线路Ⅱ： Z(0.27j0.423)5013.5j21.15

B2.6910650S1.35104S

2QB1.351041102Mvar1.63Mvar

线路Ⅱ： Z B(0.45j0.44)4018j17.6

2.5810640S1.03104S 1.031041102Mvar1.25Mvar

14.84j63.52.42j31.75 2 2QB 变电所b：ZTb S0b20.05j0.6MVA0.1j1.2MVA

变电所b：ZTc111.4j1275.7j63.5 2 S0c20.03j0.35MVA0.06j0.7MVA 等值电路如图所示

（2）计算节点b和c的运算负荷。

STb2421822.24j31.75MVA0.18j2.36MVA 2110SbSLDbSTbSobjQBIjQB24j180.18j2.360.1j1.2j0.975j0.623M24.28j19.96MVASTc122925.7j63.5MVA0.106j1.18MVA 1102ScSLDcSTcSocjQBjQB12j90.106j1.180.06j0.7j0.623j0.815MVA12.17j9.44MVA（3）计算闭式网络的功率分布。

SbZZScZ24.28j19.9631.5j38.7512.17j9.4413.5j21.15MVS47.7j64.13ZZZ18.64j15.79MVASScZZSbZZZZ12.17j19.4434.2j42.9824.28j19.9616.2j25.38MVA47.7j64.1317.8j13.6MVASIS18.64j15.7917.8j13.6MVA36.44j29.39MVA

SbSc24.28j19.9612.17j9.44MVA36.45j29.4MVA

可见，计算结果误差很小，无需重算。取S续进行计算。

S18.64j15.79MVA继SbS24.28j19.9618.65j15.8MVA5.63j4.16MVA

由此得到功率初分布，如图所示。（4）计算电压损耗。

由于线路Ⅰ和Ⅲ的功率均流向节点b为功率分点，且有功功率分点和无公功功率分点都在b点，因此这点的电压最低。为了计算线路Ⅰ的电压损耗，要用A点的电压和功率SA1。SA1SSL18.64215.8216.2j25.38MVA19.45j17.05MVA18.65j15.81102VPA1RQAX19.4516.217.0525.386.39MVA VA117变电所b高压母线的实际电压为 VbVAV1176.39110.61MVA

3.5 变比分别为k1110/11和k2115.5/11的两台变压器并联运行，如图所示，两台变压器归算到低压侧的电抗均为1Ω，其电阻和导纳忽略不计。已知低压母线电压10kV，负荷功率为16+j12MVA，试求变压器的功率分布和高压侧电压。解（1）假定两台变压器变比相同，计算其功率分布。因两台变压器电抗相等，故

S1LDS2LD11SLD16j12MVA8j6MVA 22（2）求循环功率。因为阻抗已归算到低压侧，宜用低压侧的电压求环路电势。若取其假定正方向为顺时针方向，则可得

k2 EVB10.51101kV0.5kV 10k1故循环功率为 ScVBE100.5MVAj2.5MVA ZT1ZT2j1j1（3）计算两台变压器的实际功率分布。

ST1S1LDSc8j6j2.5MVA8j8.5MVA

ST2S2LDSc8j6j2.5MVA8j3.5MVA

（4）计算高压侧电压。不计电压降落的横分量时，按变压器T-1计算可得高压母线电压为

8.51 VA10k1100.8510kV108.5kV

10按变压器T-2计算可得

3.51 VA10k2100.3510.5kV108.68kV

10 计及电压降落的横分量，按T-1和T-2计算克分别得。

VA108.79kV，VA109kV

（5）计及从高压母线输入变压器T-1和T-2的功率。S S'T1828.528j8.5j1MVA8j9.86MVA

102823.528j3.5j1MVA8j4.26MVA 210'T2输入高压母线的总功率为

S'ST'1ST228j9.868j4.26MVA16j14.12MVA 计算所得功率分布，如图所示。

3．6 如图所示网络，变电所低压母线上的最大负荷为40MW，cos0.8，Tmax4500h。试求线路和变压器全年的电能损耗。线路和变压器的参数如下：

线路（每回）：r=0.17Ω/km, x=0.409Ω/km, b2.28106S/km

变压器（每台）:P086kW,Ps200kW,I0%2.7,Vs%10.5

解 最大负荷时变压器的绕组功率损耗为

V%SSTPTjQT2PsjsSN1002SN2

2

10.540/0.82200j31500kVAkVA252j4166100231.5变压器的铁芯损耗为

I%2.7S02P0j0SN286j31500kVA kVA172j1701100100线路末端充电功率

QB22blV22.821061001102Mvar3.412Mvar

2等值电路中流过线路等值阻抗的功率为

S1SSTS0jQB240j300.252j4.1660.172j1.701j3.412MVA40.424j32.455MVA线路上的有功功率损耗

S1240.424232.455210.17100MW1.8879MW PL2RL2V1102

已知cos0.8，Tmax4500h，从表中查得3150h，假定变压器全年投入运行，则变压器全年的电能损耗 WT2P08760PT315017287602523150kWh2300520kWh

线路全年的电能损耗

WLPL31501887.93150kWh5946885kWh 输电系统全年的总电能损耗

WTWL23005205946885kWh8247405kWh

第五篇：电力系统分析潮流计算大作业

电力系统分析潮流计算大作业

（源程序及实验报告）

源程序如下：

采用直角坐标系的牛顿-拉夫逊迭代 function chaoliujisuan()m=3;%m=PQ节点个数 v=1;%v=PV节点个数

P=[-0.8055-0.18 0];%P=PQ节点的P值 Q=[-0.5320-0.12 0];%Q=PQ节点的Q值 PP=[0.5];%PP=PV节点的P值 V=[1.0];%V=PV节点的U值

E=[1 1 1 1.0 1.0]';%E=PQ,PV,Vθ节点e的初值

F=[0 0 0 0 0]';%F=PQ,PV,Vθ节点f的初值 G=[

6.3110-3.5587-2.7523 0 0;

-3.5587 8.5587-5

0 0;

-2.7523-5

7.7523 0 0;

0

0

0

0 0;

0

0

0

0 0

];B=[

-20.4022 11.3879

9.1743

0

0;

11.3879-31.00937 15

4.9889 0;

9.1743

-28.7757 0

4.9889;

0

4.9889

0

5.2493 0;

0

0

4.9889

0

-5.2493

];Y=G+j*B;X=[];%X=△X n=m+v+1;%总的节点数

FX=ones(2*n-2,1);%F(x)矩阵

F1=zeros(n-1,n-1);%F(x)导数矩阵 a=0;%记录迭代次数

EF=zeros(n-1,n-1);%最后的节点电压矩阵 while max(FX)>=10^(-5)for i=1:m %PQ节点

FX(i)=P(i);%△P FX(n+i-1)=Q(i);%△Q

for w=1:n FX(i)= FX(i)-E(i)*G(i,w)*E(w)+E(i)*B(i,w)*F(w)-F(i)*G(i,w)*F(w)-F(i)*B(i,w)*E(w);%△P

FX(n+i-1)=FX(n+i-1)-F(i)*G(i,w)*E(w)+F(i)*B(i,w)*F(w)+E(i)*G(i,w)*F(w)+E(i)*B(i,w)*E(w);%△Q

end end for i=m+1:n-1 %PV节点 FX(i)=PP(i-m);%△P FX(n+i-1)=V(i-m)^2-E(i)^2-F(i)^2;%△Q

for w=1:n FX(i)= FX(i)-E(i)*G(i,w)*E(w)+E(i)*B(i,w)*F(w)-F(i)*G(i,w)*F(w)-F(i)*B(i,w)*E(w);%△P

end end

for i=1:m %PQ节点

for w=1:n-1

if i~=w

F1(i,w)=-(G(i,w)*E(i)+B(i,w)*F(i));

F1(i,n+w-1)=B(i,w)*E(i)-G(i,w)*F(i);

F1(n+i-1,w)=B(i,w)*E(i)-G(i,w)*F(i);

F1(n+i-1,n+w-1)=G(i,w)*E(i)+B(i,w)*F(i);

else

F1(i,w)=-G(i,i)*E(i)-B(i,i)*F(i);

F1(i,n+w-1)=B(i,i)*E(i)-G(i,i)*F(i);

F1(n+i-1,w)=B(i,i)*E(i)-G(i,i)*F(i);

F1(n+i-1,n+w-1)=G(i,i)*E(i)+B(i,i)*F(i);

for k=1:n

F1(i,w)=F1(i,w)-G(i,k)*E(k)+B(i,k)*F(k);

F1(i,n+w-1)= F1(i,n+w-1)-G(i,k)*F(k)-B(i,k)*E(k);

F1(n+i-1,w)=F1(n+i-1,w)+G(i,k)*F(k)+B(i,k)*E(k);

F1(n+i-1,n+w-1)=F1(n+i-1,n+w-1)-G(i,k)*E(k)+B(i,k)*F(k);

end

end

end end for i=m+1:n-1 %PV节点

for w=1:n-1

if i~=w

F1(i,w)=-(G(i,w)*E(i)+B(i,w)*F(i));

F1(i,n+w-1)=B(i,w)*E(i)-G(i,w)*F(i);

F1(n+i-1,w)=0;

F1(n+i-1,n+w-1)=0;

else

F1(i,w)=-G(i,i)*E(i)-B(i,i)*F(i);

F1(i,n+w-1)=B(i,i)*E(i)-G(i,i)*F(i);

F1(n+i-1,w)=-2*E(i);

F1(n+i-1,n+w-1)=-2*F(i);

for k=1:n

F1(i,w)=F1(i,w)-G(i,k)*E(k)+B(i,k)*F(k);

F1(i,n+w-1)= F1(i,n+w-1)-G(i,k)*F(k)-B(i,k)*E(k);

end

end

end end X=inv(F1)*(-FX);for i=1:n-1

E(i)=E(i)+X(i);

F(i)=F(i)+X(n+i-1);end a=a+1;fprintf('第%d次迭代后的节点电压分别为:n',a);disp(E+j*F);fprintf('第%d次迭代后功率偏差△P △Q电压偏差△V的平方分别为:n',a);disp(FX);end

disp('收敛后的节点电压用极坐标表示为:');EF=E+j*F;for i=1:n-1 fprintf('%d号节点电压的幅值为:',i)

disp(abs(EF(i)));fprintf('%d号节点电压的相角度数为',i)

disp(angle(EF(i))*180/pi);end PPH=0;for i=1:n

PPH=PPH+EF(n)*conj(Y(n,i))*conj(EF(i));end fprintf('平衡节点的功率');disp(PPH);

运行结果：

运行结果复制如下：

第1次迭代后的节点电压分别为:

1.00340.1019i

1.03390.0017i

1.0000

第1次迭代后功率偏差△P △Q电压偏差△V的平方分别为:

-0.8055

-0.1800

0

0.5000

-0.3720

0.2474

0.3875

0

第2次迭代后的节点电压分别为:

0.98360.1038i

1.01830.0035i

1.0000

第2次迭代后功率偏差△P △Q电压偏差△V的平方分别为:

0.0512

-0.0222

-0.0403

0.0002

-0.1012

-0.0219

-0.0099

-0.0000

第3次迭代后的节点电压分别为:

0.98310.1038i

1.01800.0035i

1.0000

第3次迭代后功率偏差△P △Q电压偏差△V的平方分别为:

0.0008

-0.0003

-0.0005

-0.0001

-0.0021

-0.0004

-0.0003

-0.0000

第4次迭代后的节点电压分别为:

0.98310.1038i

1.01800.0035i

1.0000

第4次迭代后功率偏差△P △Q电压偏差△V的平方分别为:

1.0e-005 *

0.0280

-0.0083

-0.0164

-0.0121

-0.1085

-0.0199

-0.0135

-0.0005

收敛后的节点电压用极坐标表示为: 1号节点电压的幅值为:

0.9916

1号节点电压的相角度数为-7.4748

2号节点电压的幅值为:

1.0175

2号节点电压的相角度数为-5.8548

3号节点电压的幅值为:

1.0229

3号节点电压的相角度数为-5.5864

4号节点电压的幅值为:

1.0000

4号节点电压的相角度数为-0.2021

平衡节点的功率 0.4968-10.3280i