电力系统分析潮流计算matlab.doc

目录:
一、 软件需求说明书。。。。。...。。.。。.....。。....。...。。...。。.。...。.。。。.。.。。..。.。。。。.。..。.。。。..。.。.。。。。。.....3
二、 概要设计说明书。。..。。。。。..。.。。.。。.。..。.。。。..。。.。。。....。。.。.。..。.。。。.。..。.。..。。。。.。。.。.。。。...。。..。。4
1、编写潮流计算程序。..。。。..。..。。。..。...。。..。.。。..。。....。.。。.。......。。。。.....。...。..。.。。。。。。。4
2、数据的输入测试。.。.。...。。。。。.。.。。....。.。。.。。。...。...。。...。。。。.。.。...。.。。。.。。.。..。。...。.。....4
3、运行得出结果...。。。.。。..。。。。。。。..。。。。。...。.。.。。..。..。。..。.。..。.。。.。。...。.。。.。.。。.。。.。.。...。。..4
4、进行实验结果验证.。。。.。....。。。。。。...。.。.。。。...。。。.。。。。....。。。....。.。.。...。。...。。..。。.。。。.。4
三、 详细设计说明书。。.....。......。。.。。。.。。.。.。.。.。..。.。。.。。.。..。。..。.。。.....。..。。.。..。..。...。....。。.。。5
1、数据导入模块。。..。..。...。.。。。。..。。.。.。。。。...。.。。.。。。.。。。。。。。........。..。。....。。。...。.。。.。.。.。。5
2、节点导纳矩阵模块..。.。。。...。。...。.。.。..。.。.....。..。。.....。。。.。。。.。。。。。。。......。...。。。.....5
3、编号判断模块..。.。......。。.。.。。.....。。。。。.。.....。.。.。.。.。。..。。.。。...。.。。.。。。。。。.....。..。.。..。.5
4、收敛条件判定模块......。。。.。。。。.。。。.。。.。.。。.。..。。..。...。。。.。.。。。..。。.。.。。..。。..。。...。...。。。5
5、雅可比矩阵模块。.。.。。.。。。。。。。....。。。.。.。。.。。。。。.....。。。。.。。..。。。。。..。.。。。。。.。.。.。。..。..。。.。。。5
6、迭代计算模块.。.。....。。。。。.。。。。...。。..。..。。。。....。。。。。..。.....。。。.。。..。.。。。。。...。.。。.。.。.。。。...5
7、计算输出参数模块。。..。.。.。....。.。。.。..。。...。.。..。。。...。。.。。.....。。....。。.。.。..。。...。。。.。。.。5
四、 程序代码。..。.。。.。..。。..。。..。。......。。.。.。。。。......。.。....。。。。。.。。。.。。.。.。。。。..。。.。。。..。..。。。。。。。。...。。。6
五、 最测试例..。.。。...。。。.。。...。。.。.....。。.。.。。......。。。..。。。。。。...。。.。。.。.。。。。..。。。..。。.。.。。。....。。。。。.。.15
1、输入结果...。....。。...。。..。.......。...。....。.。...。...。。。。。..。.。.。。.。。..。。。。。....。.。。.。...。。..。。。15
2、输出结果。...。...。..。.。..。..。..。。....。。.。。.。.。。。....。。.。.。.。。...。。.。。。。。.....。。..。。..。。。..。.。..。15
3、结果验证......。...。..。。。。。。..。.....。..。..。。..。...。..。.。。.。..。。。...。。。。..。....。.。.。.。。。。.。.。。。..15
一、 软件需求说明书
本次设计利用MATLAB/C++/C(使用MATLAB)编程工具编写潮流计算，实现对节点电压和功率分布的求取。
潮流方程的求解基本方法是迭代,包括牛顿-拉夫逊法,以及P-Q分解法，本次设计采用牛顿迭代法。
牛顿迭代法（Newton’s method)又称为牛顿-拉夫逊方法,它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法.多数方程不存在求根公式，因此求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f（x) = 0的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根。 
牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代,从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说,各个母线所供负荷的功率是已知的，各个节点电压是未知的(平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的,电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组,需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导,得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量(未知的）构成了误差方程,解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
二、 概要设计说明书
1、编写潮流计算程序
本程序主要分为七个模块:数据导入模块、节点导纳矩阵模块、编号判断模块、收敛条件判定模块、雅可比矩阵模块、迭代计算模块、计算输出参数模块。
下图为潮流迭代框图.
2、 数据的输入测试
本次设计是将《电力系统分析》课本P88例题进行潮流计算。
3、 运行得出结果
得出电压、有功功率、无功功率、角度的数据。
4、进行实验结果验证
三、 详细设计说明书
1、数据导入模块：
利用Excel输入已知节点、支路数据,通过“读取”将数据导入MATLAB中.
2、节点导纳矩阵模块：
利用已知的电阻电抗及导纳的值，根据导纳的计算公式,计算出节点的自导纳及节点间的互导纳的值，按照节点编号组成导纳矩阵,利用MATLAB“real”和“imag”调用导纳矩阵中的实部和虚部,分别形成电导和电纳的矩阵。
3、编号判断模块：
当首节点在变压器左侧，设为1，位于变压器右侧，设为2，既非1也非0为不含变压器;
节点类型为PQ时，为1,节点类型为PV时，为2，节点类型为VƟ时，为3。
4、收敛条件判定模块:
根据节点的类型赋初值，并进行失配功率的初步计算,判断是否符合收敛条件max{ΔPi，ΔQi}<.如不符合，则进行后续的计算。
5、雅可比矩阵模块：
根据节点类型确定雅克比矩阵的阶数，然后根据n维非线性方程组的修正方程求出雅克比矩阵.
6、迭代计算模块：
解修正方程，并进行收敛条件判断,如不符合条件则进行下一次迭代，以一直到符合条件为止.
7、计算输出参数模块:
当满足收敛条件max{ΔPi，ΔQi｝〈时，结束迭代计算,计算电压、有功功率、无功功率、电压角度等参数.
四、 程序代码
clear %清除变量；
filename=’E：\MATLAB\zhilu。xlsx’;
a=xlsread(filename)； %读取支路信息；
filename=’E:\MATLAB\jiedian.xlsx’;
c=xlsread（filename)； ％读取节点信息
b=zeros(a（1,7)）; ％定义节点导纳矩阵
G=zeros(a(1，7));
B=zeros（a(1,7)）；
for i1=1：（a（1,7)）
if （a(i1，5）==0) ％%％首节点在变压器左侧
b(a(i1，1)，a（i1，1)）=b(a（i1，1），a（i1，1））+1/(a（i1,3)+1j*a（i1，4));
b(a(i1，2）,a(i1,2)）=b(a（i1,2）,a(i1，2））+1/（(a(i1，3）+1j＊a（i1，4)）＊a(i1,6）^2）；
b(a（i1，1),a(i1，2)）=b(a（i1，1）,a（i1，2）)+（-1/((a(i1,3)+1j*a（i1，4））＊a（i1，6））);
b（a(i1,2),a(i1,1)）=b（a(i1,2），a(i1,1））+（—1/((a（i1,3）+1j＊a(i1，4）)＊a（i1，6))）；％％％进行导纳计算
G（a(i1,1),a(i1,1）)=real（b(a（i1，1），a（i1，1)）);
G(a（i1,2），a（i1，2））=real(b(a(i1，2)，a(i1，2）)）;
G(a（i1,1）,a(i1,2））=real（b（a(i1，1)，a（i1,2）））；
G(a（i1，2），a(i1，1）)=real(b(a（i1，2），a（i1，1）））；
B(a(i1,1)，a（i1，1）)=imag(b（a(i1,1)，a（i1，1）))；
B(a（i1，2）,a(i1，2））=imag(b（a(i1,2）,a（i1，2））);
B(a(i1，1），a(i1,2）)=imag(b（a（i1，1），a(i1,2）))；
B（a(i1，2),a(i1，1))=imag(b(a（i1,2）,a(i1，1)))；
elseif （a(i1,5）==1) ％％%首节点在变压器右侧
b(a(i1,1）,a(i1,1)）=b（a（i1,1)，a（i1,1)）+1/((a(i1,3）+1j*a(i1，4))*a(i1，6)^2）;
b（a(i1，2),a（i1,2））=b(a（i1,2),a（i1,2))+1/(a（i1，3）+1j*a（i1，4));
b(a（i1，1)，a（i1，2)）=b（a(i1,1）,a(i1,2）)+（—1/（(a（i1，3）+1j＊a（i1，4))*a(i1,6）））;
b(a（i1，2）,a(i1，1))=b(a（i1,2)，a(i1，1））+(-1/（（a（i1，3)+1j＊a(i1，4)）＊a(i1,6))）; %％％进行导纳计算
G（a（i1,1)，a(i1,1））=real（b(a（i1,1),a（i1，1)））;
G（a(i1,2），a（i1,2））=real（b(a(i1，2)，a(i1,2）));
G（a（i1，1）,a（i1,2)）=real（b（a（i1,1)，a(i1,2)）);
G(a(i1,2),a（i1,1））=real（b(a（i1，2)，a（i1，1）));
B（a（i1,1），a（i1,1))=imag（b（a(i1，1)，a（i1,1）))；
B（a（i1,2)，a（i1，2）)=imag(b(a（i1,2),a（i1,2))）;
B（a（i1，1）,a（i1,2）)=imag(b（a(i1,1），a(i1,2）)）；
B（a(i1，2）,a(i1,1））=imag(b(a（i1,2),a(i1,1）)）；
else %%％不含变压器支路
b(a(i1，1)，a（i1，1）)=b（a（i1,1），a(i1,1）)+1/（a（i1,3)+1j*a（i1,4))+1j＊a(i1,5）;
b(a(i1,2），a(i1,2））=b(a(i1，2),a（i1，2））+1/（a（i1，3）+1j*a(i1,4）)+1j*a(i1,5）;
b（a(i1，1）,a(i1,2）)=b（a（i1，1),a(i1，2))+（-1/((a(i1，3）+1j*a（i1,4））)）；
b(a(i1,2)，a（i1，1）)=b(a(i1，2）,a（i1，1））+(-1/((a(i1，3)+1j*a（i1，4）））)； %％%进行导纳计算
G(a（i1，1），a(i1,1)）=real（b（a（i1，1），a（i1，1)));
G(a(i1，2）,a（i1，2))=real(b(a（i1,2）,a(i1，2)）)；
G(a(i1，1）,a（i1，2))=real(b（a（i1，1)，a(i1，2))）；
G(a（i1,2),a（i1，1))=real（b（a（i1,2),a（i1,1）)）;
B(a（i1,1)，a(i1,1）)=imag（b(a（i1,1），a（i1,1)））;
B（a（i1,2),a(i1,2）)=imag（b（a（i1，2）,a（i1,2)））；
B（a(i1，1）,a（i1,2））=imag(b（a（i1,1),a(i1，2)））;
B(a(i1,2）,a（i1,1)）=imag（b(a（i1,2)，a(i1，1）））;
end
end
％％％％计数各个节点个数
PQjd=0;
PVjd=0;
VOjd=0;
for i2=1：a（1,7）
if c(i2，6）==1
PQjd=PQjd+1;
end
if c(i2，6）==2
PVjd=PVjd+1;
end
if c(i2,6)==3
VOjd=VOjd+1;
end
end
%注入功率的计算
zhuru_PQ=zeros（PQjd＊2+PVjd,1）;
k=1；
for i3=1：a（1，7）
if c(i3，6）==1
zhuru_PQ（k，1)=c(i3,2）-c（i3，4);
zhuru_PQ（k+a(1,7）—1,1)=0;
for i10=1：a(1,7)
zhuru_PQ（k+a(1，7）-1,1)=zhuru_PQ（k+a（1，7）-1,1）+B(i3，i10)；
end
zhuru_PQ(k+a(1，7)—1,1)=zhuru_PQ（k+a(1,7）—1,1)-c(i3,5）
k=k+1；
end
if c(i3，6)==2
zhuru_PQ(k,1）=c（i3，2）-c（i3，4);
k=k+1；
end
end
n=a(1,7)；
%%%%％后续矩阵变量定义
angle_u=zeros(PQjd*2+PVjd，1）； %迭代角度、电压矩阵
angle_u（n：PQjd*2+PVjd,1）=1; %电压初值设为1,角度初值设为0
delta_au=zeros(PQjd＊2+PVjd，1）；
delta_au（n：PQjd*2+PVjd，1）=0；
%角度、电压修正量向量
%％％%%失配功率的计算
k=1；
k_=PQjd+1；
delta_PQ= zeros(PQjd＊2+PVjd，1）;
for i4=1：n
％%%%％％%％％％%%%％％％%％％PQ节点的计算
if c（i4，6)==1
delta_PQ(k，1)=zhuru_PQ（k，1)；
delta_PQ(k+n—1,1）=zhuru_PQ(k+n—1,1);
i6=1；
for i5=1:n
if c（i5,6)==1
delta_PQ(k,1）=delta_PQ（k，1)-angle_u(n-1+k）＊angle_u(n-1+i6)＊（G(i4,i5）*cos(angle_u(k,1)-angle_u(i6,1)）+B(i4,i5）＊sin(angle_u(k,1)-angle_u（i6,1）)）；
delta_PQ（k+n—1,1）=delta_PQ(k+n-1,1）-angle_u(n—1+k)*angle_u（n—1+i6）＊（G（i4,i5)*sin（angle_u(k,1)-angle_u（i6,1)）—B(i4,i5）*cos（angle_u（k,1)—angle_u(i6，1）））;
i6=i6+1;
end
if c(i5，6）==2
delta_PQ（k,1)=delta_PQ（k，1）-angle_u（n-1+k）＊c（i5,1）＊（G（i4,i5）＊cos(angle_u（k，1)-angle_u（i6,1）)+B（i4，i5）＊sin（angle_u（k，1）-angle_u(i6,1）)）；
delta_PQ（k+n—1,1)=delta_PQ(k+n—1,1）—angle_u(n-1+k)＊c（i5,1)*（G（i4，i5）*sin（angle_u（k,1）-angle_u（i6,1））-B（i4,i5）*cos(angle_u（k,1)—angle_u（i6，1）));
i6=i6+1;
end
if c(i5，6)==3
delta_PQ(k，1)=delta_PQ(k,1）-angle_u(n—1+k）*c(i5,1)＊(G（i4,i5）*cos(angle_u(k，1))+B(i4，i5)*sin(angle_u(k,1）））；
delta_PQ(k+n—1，1)=delta_PQ(k+n-1,1)—angle_u(n-1+k)*c（i5，1)＊(G（i4，i5)*sin(angle_u（k,1））-B（i4,i5)＊cos(angle_u(k,1)）)；
end
end
k=k+1;
end
%%%％%％％％%％％%%%%%％％%％%PV节点的计算
if c(i4，6)==2
delta_PQ(k_，1)=zhuru_PQ（k_,1）；
i6=1；
for i5=1:n
if c(i5，6）==1
delta_PQ（k_，1）=delta_PQ（k_，1）-c(i4,1)*angle_u(n-1+i6)＊(G（i4,i5）＊cos(angle_u（k_,1)—angle_u(i6,1))+B(i4,i5）＊sin(angle_u（k_,1）—angle_u(i6,1)））;
i6=i6+1;
end
if c(i5,6)==2
delta_PQ(k_,1）=delta_PQ(k_,1)-c（i4,1)＊c(i5,1）＊（G(i4,i5）*cos（angle_u(k_,1)—angle_u(i6，1）)+B（i4,i5)＊sin(angle_u(k_,1）—angle_u(i6，1）)）；
i6=i6+1;
end
if c(i5，6)==3
delta_PQ(k_，1）=delta_PQ（k_,1）-c（i4,1)*c（i5,1)*（G（i4,i5)＊cos（angle_u（k_,1)）+B（i4,i5）*sin（angle_u（k_,1）))；
end
end
k_=k_+1;
end
end
zhongjian_delta_PQ=abs（delta_PQ）； ％％%％取绝对值
jingdu=max(zhongjian_delta_PQ）; %%％%为迭代精度的判断做赋值准备
while （jingdu>0。0001)
％％％％%％%雅各比矩阵的形成
J=zeros(2＊PQjd+PVjd); %％%%矩阵定义
U=ones(n,1)； %%%％为方便运算将电压取出放入一个新的矩阵
for i=1:PQjd
U（i，1）=angle_u（i+n-1，1)；
end
angle=zeros(n,1）; ％%％％为方便运算将角度取出放入一个新的矩阵
for i=1：n-1
angle(i，1）=angle_u(i,1）；
end
%％％％%H%％％％% %％％％%H部分生成
for i1=1:n-1
for i2=1：n-1
if i1～=i2
J（i1，i2）=—U（i1，1)*U(i2,1）*(G（i1，i2)＊sin（angle(i1，1）-angle(i2,1）)-B（i1，i2）＊cos(angle(i1，1)-angle（i2，1))）；
end
if i1==i2
for i3=1：n
if i3～=i1
J（i1，i2）=J（i1,i2）+U(i1，1)＊U(i3，1)*(G（i1，i3)＊sin(angle(i1,1）-angle（i3,1））-B（i1，i3)*cos(angle(i1，1）-angle（i3,1))）；
end
end
end
end
end
％%%%%N%％％%％ ％％%％％N部分生成
for i1=1:n—1
for i2=1：PQjd
if i1～=i2
J(i1,i2+n-1)=-U(i1，1)*(G(i1，i2)*cos（angle(i1，1）—angle（i2,1))-B（i1,i2）＊sin（angle（i1，1)-angle（i2，1)）)；
end
if i1==i2
J(i1，i2+n—1)=(—2)*U（i1，1）＊U(i1，1）*G（i1,i1）；
for i3=1：n
if i3~=i1
J(i1,i2+n-1）=J(i1,i2+n—1)-U(i3，1)*(G(i1，i3)*cos（angle(i1，1）—angle（i3，1））+B（i1,i3)＊sin(angle（i1，1)—angle(i3，1))）;
end
end
end
end
end
%%%％%K%％%%％ %％％%％K部分生成
for i1=1：PQjd
for i2=1:n-1
if i1~=i2
J(i1+n—1，i2）=U（i1，1)＊U（i2，1）＊（G（i1,i2）*cos（angle(i1,1）—angle(i2,1)）-B（i1,i2）*sin(angle（i1，1）-angle（i2，1））);
end
if i1==i2
for i3=1:n
if i3～=i1
J(i1+n—1，i2）=J（i1+n-1,i2)-U（i1，1）＊U（i3，1）＊(G（i1，i3)*cos(angle（i1，1）-angle（i3,1))+B(i1，i3）＊sin（angle(i1,1）—angle(i3，1)））；
end
end
end
end
end
%%%％％L%%%％% %％%％%L部分生成
for i1=1:PQjd
for i2=1:PQjd
if i1~=i2
J(i1+n—1,i2+n—1)=—U（i1,1）*(G（i1，i2)＊sin（angle(i1,1)-angle(i2，1）)-B（i1,i2)*cos(angle(i1，1）—angle(i2，1））);
end
if i1==i2
J（i1+n-1,i2+n-1）=2＊U(i1，1）＊U(i1,1）＊B(i1，i1）;
for i3=1：n
if i3～=i1
J（i1+n-1,i2+n-1)=J(i1+n-1，i2+n-1)—U（i3,1）*（G(i1，i3)＊sin(angle（i1,1）-angle(i3，1)）-B（i1,i3）＊cos（angle（i1,1)—angle（i3，1）));
end
end
end
end
end
delta_au=-inv（J)＊delta_PQ； ％%%解修正方程
angle_u=angle_u+delta_au； ％％％形成新的初始值
％%％％%％%%％%％%%%%%％％%％%％％％%%％％%%%%％％%％％％
％再次计算失配功率
k=1;
k_=PQjd+1；
delta_PQ= zeros(PQjd＊2+PVjd,1）；
for i4=1:n
％％%％％%％%％％％％%%％％％％%
if c(i4,6）==1
delta_PQ(k，1）=zhuru_PQ(k,1）;
delta_PQ(k+n-1，1)=zhuru_PQ（k+n-1，1)；
i6=1；
for i5=1:n
if c（i5，6）==1
delta_PQ(k,1）=delta_PQ(k,1）—angle_u（n-1+k)＊angle_u（n-1+i6)＊(G（i4，i5）*cos（angle_u(k,1)—angle_u(i6，1)）+B(i4，i5）*sin（angle_u(k，1）—angle_u（i6,1））)；
delta_PQ（k+n—1,1)=delta_PQ(k+n-1,1)—angle_u(n-1+k)＊angle_u(n-1+i6)*（G（i4,i5)*sin(angle_u（k,1）—angle_u（i6,1)）-B（i4,i5）＊cos(angle_u(k，1)-angle_u（i6,1）));
i6=i6+1;
end
if c（i5,6）==2
delta_PQ(k,1)=delta_PQ（k,1)-angle_u(n-1+k)*c（i5，1)*(G（i4，i5)*cos（angle_u（k，1)-angle_u(i6，1))+B(i4,i5)*sin(angle_u(k,1）-angle_u(i6,1)));
delta_PQ(k+n-1，1）=delta_PQ（k+n—1,1)—angle_u（n-1+k)*c（i5，1)＊（G（i4，i5）*sin(angle_u（k，1）—angle_u（i6,1))—B（i4，i5)*cos（angle_u(k,1）-angle_u（i6,1)）)；
i6=i6+1；
end
if c(i5,6)==3
delta_PQ(k,1）=delta_PQ（k，1）-angle_u(n-1+k)*c（i5,1)＊（G(i4，i5）＊cos(angle_u(k，1))+B（i4,i5)＊sin(angle_u（k,1）));
delta_PQ（k+n-1，1）=delta_PQ（k+n—1，1）-angle_u(n—1+k）＊c（i5，1)*（G(i4,i5）*sin(angle_u（k，1）)-B(i4，i5)*cos（angle_u(k，1）))；
end
end
k=k+1;
end
%%%％％%％%％％％％％%％%%%%%％
if c（i4，6)==2
delta_PQ(k_，1）=zhuru_PQ（k_,1）;
i6=1；
for i5=1:n
if c（i5，6）==1
delta_PQ（k_，1)=delta_PQ（k_，1）—c（i4，1）＊angle_u（n—1+i6）＊（G(i4，i5）*cos（angle_u(k_，1）—angle_u（i6,1))+B(i4,i5)*sin(angle_u（k_,1)—angle_u（i6，1）))；
i6=i6+1;
end
if c(i5，6）==2
delta_PQ（k_，1）=delta_PQ(k_,1)-c（i4,1)＊c(i5，1）*(G（i4,i5）＊cos（angle_u（k_，1）-angle_u(i6，1））+B(i4，i5)*sin（angle_u(k_，1）—angle_u（i6，1)）);
i6=i6+1;
end
if c（i5,6）==3
delta_PQ（k_,1）=delta_PQ（k_,1）—c（i4，1）*c(i5，1）＊(G（i4，i5）＊cos（angle_u(k_，1))+B（i4，i5)*sin（angle_u(k_，1）））；
end
end
k_=k_+1；
end
end
%％％％再次进行精度判断前的准备
m=0；
zhongjian_delta_PQ=delta_PQ;
zhongjian_delta_PQ=abs（delta_PQ）;
jingdu=max(zhongjian_delta_PQ);
m=m+1;
end
%%％％%％迭代完毕
%%%％％％电压与角度换算结果输出
U=ones(n,1）;
for i=1：PQjd
U（i,1)=angle_u(i+n-1，1）;
end
angle=zeros(n,1);
for i=1：n-1
angle(i,1)=angle_u（i，1）；
end
angle_jiaodu=angle/（2＊pi）＊360；
％%％%％%%％％计算各个节点功率
P=zeros（n,1)；
Q=zeros(n,1）;
for i1=1：n
for i2=1:n
P（i1，1)=P(i1，1)+U（i1，1)*U(i2，1）＊(G(i1，i2)＊cos(angle（i1，1）—angle(i2,1)）+B(i1,i2)*sin(angle(i1,1）-angle（i2，1）)）；
Q（i1,1）=Q(i1，1）+U(i1，1）＊U（i2，1）*(G（i1,i2)*sin（angle(i1，1）-angle（i2,1)）—B(i1，i2)*cos（angle（i1,1）—angle（i2,1))）；
end
end
五、最测试例
本次设计试例为《电力系统分析》课本P88例题3-3.
1、输入数据
支路数据如下：
节电数据如下:
2、 输出结果
J矩阵
k 电压U 角度Ɵ 有功功率P 无功功率Q
3、 结果验证
输出结果与课本P90结果一致，潮流计算程序正确.
20