大学毕设论文--复杂网络nr法潮流分析与计算的设计课程设计.doc

电气工程及其自动化专业课程设计
复杂网络N-R法潮流分析与计算的设计

学生学号：
学生姓名：
班 级：
指导教师：
起止日期：

哈尔滨工程大学自动化学院

课程设计报告撰写内容
一、设计要求（宋体，小四号字，加黑）
用matlab编程，N_R法计算潮流分布
具体要求为：
（1）给出程序，并给出注释
（2）输出迭代次数，各节点电压，各支路电流
（3）在图中标明功率流向

节点数据如下表所示（标幺值）

1
2
3
4
5
6
P

3
1.8
0.6
3.5
5
Q

1
0.5
0.8
1.3

V
1

1.05

支路及变压器数据
线路
T1
T2
L2
L3
L4
L5
阻抗
j0.04
j0.02
0.06+j0.025
0.01+j0.2
0.06+j0.5
0.05+j0.3
导纳/2

j0.25
j0.25
j0.25
j0.25
变比
1.05:1
1.05:1

精度要求：0.0001

二、设计方案（要求给出详细的设计思路及其必要的论证）
（1.）潮流计算的方法
（1）高斯雅克比迭代法
（2）高斯-塞得尔法（对初值要求底，迭代次数多）
（3）牛顿-拉夫逊法（使用广泛）
（4）PQ 快速分解法（提升运算速度）
目前广泛应用的潮流计算方法都是基于节点电压法的，以节点导纳矩阵Y
作为电力网络的数学模型。节点电压Ui 和节点注入电流Ii 由节点电压方程
YV=I （1）
根据S=VI﹡(I﹡为I 的共轭)可得非线性的节点方程
YV=I=(S/V)﹡ （2）
在实际的电力系统中，已知的运行条件不是节点的注入电流，而是负
荷和发电机的功率，而且这些功率一般不随节点电压的变化而变化。由于
各节点注入功率与注入电流的关系为Si＝Pi＋jQi =ViIi﹡，因此可将式（2）
改写为
Ii=Si/Vi=Pi+jQi/Vi (i= 1, 2,3 ⋯ n) （3）
式中，Pi 和Qi 分别为节点i 向网络注入的有功功率和无功功率，当
i 为发电机节点时Pi﹥0；当i 为负荷节点时Pi﹤0；当i 为无源节点Pi＝
0，Qi＝0；Vi 和Ii 分别为节点电压相量Vi 和节点注入电流相量Ii 的共
轭。
式（3）亦即潮流计算的基本方程式,它可以在直角坐标也可以在极坐
标上建立2n 个实数形式功率方程式。发电机Pi、Qi 为正，负荷Pi、Qi 为
负。
展开YV=I 为
Ii=ΣYijVj=YiiVi+ΣYijVi( i=1 2 3 ⋯n) (4)
将式（4）代入式（3），得n 维的非线性复数的电压方程组
潮流计算的基本方程为
(Pi-jQi)/ Vi= YiiVi+ΣYijVi （i=1, 2, ⋯ n） (5)
（2.）变量的分类
假设系统中有n 个节点，构成n 个复数方程，2n 个实数方程，变量总数为6n 个。
a)不可控变量（2n 个）：负荷消耗的有功功率Li P 和无功功率Li Q .由于该类变
量无法控制，取决于用户，而且出现事先没有预计的变动，使系统偏离原始运行
状态，因此又称为不可控变量或扰动变量。
b)控制变量（2n 个）:发电机发出的有功功率Gi P 和无功功率Gi Q ，因为该类
变量可控。也称独立变量。
c)状态变量（2n 个）：母线电压或节点电压的幅值大小i V 与相角大小i δ ，又
称依从变量或因变量。并且i V 受Gi P 控制， i δ 受Gi Q 控制。
其中2n 个扰动变量是给定的，2n 个控制变量和2n 个状态变量中给定两个，求
另外两个。
（3.）变量的约束条件
a)扰动变量没有约束条件。
b)控制变量约束条件：为满足发电机的技术经济特性指标。

c)状态变量的i V 的约束条件：保证良好的电能质量。

d) 状态变量的i δ 的约束条件：保证系统的稳定运行。

(4.)系统节点的分类，根据给定的控制变量和状态变量进行分类如下：
（1）PQ 节点（即负荷节点）：
给定Pi、Qi，求Vi 和i δ （ i i e , f ）。通常变电所都是这一类型的节点，由于
没有发电设备，因而发电功率为零电力系统中的绝大多数节点属于这一节点。其
包含变电站节点（即联络节点或浮游节点）。

（2）PV 节点（即调节节点、电压控制节点）：
给定Pi 和Vi，求Qi 和i δ （ i i e , f ）。这类节点必须有足够的可调无功容量，
用以维持给定的电压幅值。一般时选择有一顶武功储备的发电厂和具有可调无功
电源设备的变电所作为PV 节点。在电力系统中，这类节点数很少。
（3）平衡节点（即松弛节点、参考节点、基准节点）：
给定Vi 和i δ （ i δ =0），求Pi 和Qi。（只有一个）有功功率不能给定，这个节
点承担了系统的有功功率平衡。同时其电压幅值也是给定的，相位为零。

(5. )P-Q 分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思
想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率
误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，
把有功功率和无功功率迭代分开来进行。
牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系
统时，修正方程式展开为：
ΔP = HΔΘ + NΔV/ V
ΔQ = JΔΘ + LΔV /V
以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特
点。
电力系统中有功功率主要与各节点电压向量的角度有关，无功功率则主要受
各节点电压幅值的影响。大量运算经验也告诉我们，矩阵N 及J 中各元素的数
值相对是很小的，因此对牛顿法的第一步简化就是把有功功率和无功功率分开
来进行迭代，即将式(4)化简为：
ΔP = HΔΘ
ΔQ = LΔV /V (5)
这样，由于我们把2n 阶的线性方程组变成了二个n 阶的线性方程组，对牛
顿法的第二个化简，也是比较关键的一个化简，即把式(5)中的系数矩阵简化为
在迭代过程中不变的对称矩阵。众所周知，一般线路两端电压的相角差是不大
的(通常不超过10~20 度)，因此可以认为：
(6)
此外，与系统各节点无功功率相应的导纳Li B 必定远远小于该节点自导纳的
虚部，即：

因此， (7)
考虑到以上关系后，式(5)中系数矩阵中的元素表达式可以化简为：
（8）
这样，式(5)中系数矩阵可以表示为：
（9）
进一步可以把它们表示为以下矩阵的乘积：
（10）
将它代入(5)中，并利用乘法结合率，我们可以把修正方程式变为：

将以上两式的左右两侧用以下矩阵左乘

就可得到

以上两式就是P-Q 分解法达到修正方程式，其中系数矩阵只不过是系统导纳
矩阵的虚部，因而是对称矩阵，而且在迭代过程中维持不变。它们与功率误差
方程式

构成了P-Q 分解法迭代过程中基本计算公式，其迭代步骤大致是：
(1)给定各节点电压向量的电压初值V i (0) ,θi (0);
(2)根据(12)计算各节点有功功率误差ΔPi ，并求出;ΔPi/Vi
(3)解修正方程式(11)，并进而计算各节点电压向量角度的修正量i Δθ
(4)修正各节点电压向量角度θi ；

(5)根据式(16)计算各节点无功功率误差i ΔQ ，并求出 / ; i i ΔQ V
(6)解修正方程式(11)，求出各节点电压幅值的修正量i ΔV
(7)修正各节点电压幅值i V (i ) (i 1) (i 1)
i i i V =V − − ΔV − (18)
(8)返回(2)进行迭代，直到各节点功率误差及电压误差都满足收敛条件。
P-Q 分解法与牛顿法潮流程序的主要差别表现在它们的修正方程式上。P-Q
分解法通过对电力系统具体特点的分析，对牛顿法修正方程式的雅可比矩阵进
行了有效的简化和改进，有以下三个特点：
（1）在提高计算速度和减少内存方面的作用是明显的，不再叙述。
（2）使我们得到以下好处。首先，因为修正方程式的系数矩阵就是导
纳矩阵的虚部，因此在迭代过程中不必象牛顿法那样进行形成雅可比矩阵的
计算，这样不仅是仅减少了运算量，而且也大大简化了程序。其次，由于系数
矩阵在迭代过程中维持不变，因此在求解修正方程式时，可以迅速求得修正量，
从而显著提高了迭代速度。
（3）可以使我们减少形成因子表时的运算量，而且由于对称矩阵三角分解
后，其上三角矩阵和下三角矩阵有非常简单的关系，所以在计算机中可以只储
存上三角矩阵或下三角矩阵，从而也进一步节约了内存。
三、设计内容
%本程序的功能是用牛顿——拉夫逊法进行潮流计算
% B1矩阵：1、支路首端号； 2、末端号； 3、支路阻抗； 4、线路对地电纳 (或变压器导纳)；
% 5、支路的变比； 6、支路首端处于K侧为1，1侧为0；
% 7、线路/变压器标识（0/1）变压器参数当支路首端处于K侧标识为1时归算至末端侧，0归算至首端侧
% B2矩阵：1、该节点发电机功率； 2、该节点负荷功率；
% 3、PQ节点电压初始值，或平衡节点及PV节点电压的给定值
% 4、节点所接无功补偿并联电容（感）的电纳
% 5、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；3为PV节点；
clear;
isb=1; %input('请输入平衡母线节点号：isb=');
pr=0.0001; %input('请输入误差精度：pr=');
%---------------------------------------------------
n=6;%input('请输入节点数：n=');
nl=6;%input('请输入支路数：nl=');
B1=[1 2 0+0.04i 0 1.05 1 1;
2 3 0.06+0.025i 0+0.5i 1 0 0;
2 5 0.01+0.2i 0+0.5i 1 0 0;
3 4 0.06+0.50i 0+0.5i 1 0 0;
4 5 0.05+0.3i 0 1 0 0;
6 5 0+0.02i 0 1.05 1 1]
B2=[0 0 1 0 1;
0 3+1i 1.00 0 2;
0 1.8+0.50i 1.00 0 2;
0 0.6+0.8i 1.00 0 2;
0 3.5+1.3i 1.00 0 2;
0 -5+0i 1.05 0 3]
%input('请输入各节点参数形成的矩阵： B2=');
%X=[1 0;2 0;3 0;4 0;5 0;6 0]
%-------------------------------------------------------------
%n=4;%input('请输入节点数：n=');nl=4;%input('请输入支路数：nl=');
%B1=[1 2 4+16i 0 1 0 0;1 3 4+16i 0 1 0 0;2 3 2+8i 0 1 0 0;2 4 1.49+48.02i 0 11/110 0 1] %input('请输入由支路参数形成的矩阵： B1=');
%B2=[0 0 115 0 1;0 0 110 0 2;0 20+4i 110 0 2;0 10+6i 10 0 2] %input('请输入各节点参数形成的矩阵： B2=');
%-------------------------------------------------------------
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
% % %-----------求导纳矩阵------------------------
%for i=1:n
% if X(i,2)~=0;
% p=X(i,1);
% Y(p,p)=1/X(i,2);
%end
%end
for i=1:nl %从1到n1（总支路数）
if B1(i,7)==1 %-----------如果是变压器支路--------
if B1(i,6)==0 %左节点(首端)处于1侧
p=B1(i,1);q=B1(i,2);
else %左节点(首端)处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
Y(q,p)=Y(p,q); %非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2); %对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); %对角元1侧+励磁导纳
else %------------否则为线路支路--------------------
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(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./2.0000; %对角元j侧+线路电纳的一半
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2.0000; %对角元i侧+线路电纳的一半
end
end
disp('导纳矩阵 Y=');disp(Y);
%-----------给定各节点初始电压及给定各节点注入功率--------------------------
G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部
for i=1:n %给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=abs(B2(i,3)); %PV、平衡节点及PQ节点电压模值
end
for i=1:n %给定各节点注入功率
S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,4); %i节点无功补偿量(电纳值)
end
%==================用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）=======================
P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2
while IT2~=0 % N0=2*n 雅可比矩阵的阶数；N1=N0+1扩展列
IT2=0;a=a+1;
JZ=['Jacobi矩阵第(',num2str(a),')次消去运算'];JZ1=['Jacobi矩阵第(',num2str(a),')次回代运算'];JZ0=['功率方程第(',num2str(a),')次差值：'];
%----------------求取各个节点的功率及功率偏差及PV节点的电压偏差--------------------
for i=1:n %n个节点2n行(每节点两个方程P和Q或U)
p=2*i-1;m=p+1;C(i)=0;D(i)=0;
for j1=1:n %第i行共n列(n个节点间互导纳及节点电压相乘即电流)
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
%求i节点有功和无功功率P',Q'的计算值
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)
V2=e(i)^2+f(i)^2; %电压模平方
%===求取功率差及PV节点电压模平方差 =========
if i~=isb %非平衡节点(PQ或PV节点)
if B2(i,5)~=3 %非PV节点(只能是PQ节点)
J(m,N1)=P(i)-P1; %PQ节点有功功率差J(m,N1)扩展列△P
J(p,N1)=Q(i)-Q1; %PQ节点无功功率差J(p,N1)扩展列△Q
else %PV节点==================
J(m,N1)=P(i)-P1; %PV节点有功功率差J(m,N1)扩展列△P
J(p,N1)=V(i)^2-V2; %PV节点电压模平方差J(p,N1)扩展列△U
end
end %(if i~=isb) 非平衡节点(PQ或PV节点)
end %(for i=1:n) n个节点2n行(每节点两个方程P和Q或U)
for m=1:N0
JJN1(m)=J(m,N1);
end
disp(JZ0);disp(JJN1);
%-------------判断功率偏差量及PV节点的电压偏差量是否满足要求-----------------
for k=3:N0 %除去平衡节点1、2号以外的所有节点
DET=abs(J(k,N1));
if DET>=pr; %PQ节点的功率偏差量及PV节点的电压偏差量是否满足要求
IT2=IT2+1; %不满足要求的节点数加1

end
end
ICT2(a)=IT2; %不满足要求的节点数；a为迭代次数
ICT1=ICT1+1; %迭代次数
if ICT2(a)==0; %当前不满足要求的节点数为零
break %退出迭代运算
end
%--------------------以上为求取各个节点的功率及功率偏差及PV节点的电压偏差-------------
%================= 求取Jacobi矩阵形成修正方程 ===================
for i=2:n %n个节点2n行(每节点两个方程P和Q或U)
if i~=isb %非平衡节点(PQ或PV节点)
if B2(i,5)~=3 %下面是针对PQ节点来求取Jacobi矩阵的元素 ===========
C(i)=0;D(i)=0;
for j1=1:n %第i行共n列(n个节点间互导纳及节点电压相乘即电流)
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
for j1=2:n %第i行共n列(2n个Jacobi矩阵元素dP/de及dP/df或dQ/de及dQ/df)
if j1~=isb&j1~=i %非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % X1=dP/de=-dQ/df=-X4
X2=B(i,j1)*e(i)-G(i,j1)*f(i); % X2=dP/df=dQ/de=X3
X3=X2; % X2=dp/df X3=dQ/de
X4=-X1; % X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;m=p+1; % X3=dQ/de J(p,N)=DQ节点无功功率差 J(p,N)=DQ;
J(m,q)=X1;q=q+1; % X1=dP/de J(m,N)=DP节点有功功率差 J(m,N)=DP;
J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;%扩展列△Q J(p,N)=DQ;
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;%扩展列△P J(m,N)=DP;
J(m,q)=X2;
end
end
else %if B2(i,5)~=3 % 否则（即为PV节点）
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
for j1=1:n
if j1~=isb&j1~=i %非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5; % PV节点电压误差J(p,N)=DV;
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,N)=DP;
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5; % PV节点电压误差J(p,N)=DV;
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X6; % PV节点有功误差J(m,N)=DP;
J(m,q)=X2;
end
end
end %(if B2(i,5)~=3 else)
end %(if i~=isb)
end %(for i=1:n)n个节点2n行(每节点两个方程P和Q或U)
JZ0=['形成的第(',num2str(a),')次Jacobi矩阵：'];
disp(JZ0);disp(J);
%=============================== 以上为形成完整的Jacobi矩阵 ================================
%====下面用高斯消去法对由Jacobi矩阵形成的修正方程进行求解(按列消去、回代) ==========
for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）
for k1=k+1:N1 % 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U
J(k,k1)=J(k,k1)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1; % 对角元规格化K行K列对角元素赋1
%================== 按列消去运算 ==================================
for k2=k+1:N0 % 从k+1行到2*n最后一行
for k3=k+1:N1 %从k2+1列到扩展列消去k+1行后各行下三角元素
J(k2,k3)=J(k2,k3)-J(k2,k)*J(k,k3);%消去运算
end %用当前行K3列元素减去当前行K列元素乘以第k行K3列元素
J(k2,k)=0; %当前行第k列元素已消为0
end
end
JZ=['Jacobi矩阵第(',num2str(a),')次消去运算'];JZ1=['Jacobi矩阵第(',num2str(a),')次回代运算'];
disp(JZ);disp(J);
%==================== 按列回代运算 =======================================
for k=N0:-1:3
for k1=k-1:-1:3
J(k1,N1)=J(k1,N1)-J(k1,k)*J(k,N1);
J(k1,k)=0;
end
end
for m=1:N0
JJN1(m)=J(m,N1);
end
disp(JZ1);disp(JJN1);%disp(J);
%----------------------------------修改节点电压-------------------------------
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N1); %修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N1); %修改节点电压虚部
U(L)=sqrt(e(L)^2+f(L)^2);
end
disp('各个节点电压模');disp(U);
%============================== 结束一次迭代 ==============================
end
%********************** 下面为迭代计算结束后的有关输出过程 *****************
disp('迭代次数：');
disp(ICT1-1);
disp('没有达到精度要求的个数：');
disp(ICT2);
for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度
E(k)=e(k)+f(k)*j; %将各节点电压用复数表示
end
%=============== 计算各输出量 ===========================
disp('各节点的电压复数值E为(节点号从小到大排列)：');
disp(E); %显示各节点的实际电压值E用复数表示
disp('-----------------------------------------------------');
disp('各节点的电压模值大小V为(节点号从小到大排列)：');
disp(V); %显示各节点的电压大小V的模值
disp('-----------------------------------------------------');
disp('各节点的电压相角sida为(节点号从小到大排列)：');
disp(sida); %显示各节点的电压相角
for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点注入电流的共轭值
end
S(p)=E(p)*C(p); %计算各节点的功率 S = 电压 X 注入电流的共轭值
end
disp('各节点的功率S为(节点号从小到大排列)：');
disp(S); %显示各节点的注入功率
disp('-----------------------------------------------------');
disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,7)==0
Si(p,q)=E(p)*conj(E(p)*B1(i,4)./2+(E(p)-E(q))./B1(i,3));
Siz(i)=Si(p,q);
else
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p)*B1(i,4)...
+(E(p)*B1(i,5)-E(q))*(1./(B1(i,3)*B1(i,5)))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*conj((E(p)-E(q)*B1(i,5))*(1./(B1(i,3)*B1(i,5)^2)));
Siz(i)=Si(p,q