电力系统潮流程序设计.doc

潮流程序的设计与开发
1 数据结构的设计
要求：将设备铭牌参数和有名值作为原始输入数据，潮流结果以有名值的形式输出。

支路数据与节点数据分别以一个矩阵的形式输入，矩阵的每一行表示每一个节点或每一条支路，矩阵的每一列表示不同的参数数据。

2 变量命名设计
变量名称
程序变量表示
变量名称
程序变量表示
节点导纳矩阵

无功功率

电导

视在功率

电钠

雅克比矩阵

电压幅值

平衡节点编号

电压相角

有功不平衡量

有功功率

无功不平衡量

电压相角修正量

电压幅值修正量

3 程序流程图

4 程序代码
4.1 主程序
n=input('请输入节点数:n=');
l=input('请输入支路数:l=');
%支路数不要包括三绕组变压器
sw=input('请输入平衡节点号:sw=');
ac=input('请输入误差精度:ac=');
SB=input('请输入基准功率:SB=');
B1=input('请输入支路参数:B1=')
%支路参数节点参数和对地支路参数均以矩阵形式储存
%第一列储存支路编号
%第二列与第三列分别储存支路的两个端点，分别为p,q
%第四列储存支路阻抗
%第五列储存支路对地导纳,注意对地导纳不要除以2
%第六列储存该支路是否含有变压器，有为1，无为0
%第七列储存变压器变比k,p指向q的变压器变比为k:1，且k大于等于1
%第八列储存变压器短路损耗
%第九列储存变压器短路电压百分数
%第十列储存变压器空载损耗
%第十一列储存变压器空载电流百分数
%第十二列储存变压器低压侧电压
%第十三列储存变压器额定功率
%第十四列储存归算所取基准电压
%注意，将三绕组变压器转换为双绕组变压器输入
A1=input('请输入节点参数:A1=');
%第一列为节点编号
%第二列为注入发电功率
%第三列为负荷功率
%第四列为节点电压幅值,为方便起见，以标幺值形式表示
%第五列为节点电压相角
%第六列储存节点对地导纳
%第七列为节点的类型，1为PQ节点，2为PV节点，3为平衡节点

%%%首先求解双绕组变压器参数
for i=1:l
if B1(i,6)==1
Zt(i)=B1(i,8)*B1(i,12)^2/(1000*B1(i,13)^2);
Xt(i)=B1(i,9)*B1(i,12)^2/(100*B1(i,13));
Gt(i)=B1(i,10)/(1000*B1(i,12)^2);
Bt(i)=B1(i,11)*B1(i,13)/(100*B1(i,12)^2);

end
end

B2(:,1:7)=B1(:,1:7);
for i=1:l
B2(i,8)=Zt(i);
B2(i,9)=Xt(i);
end

for i=1:l
if B1(i,6)~=1
ZB(i)=B1(i,14)^2/SB;
YB(i)=1/ZB(i);
B2(i,4)=B1(i,4)/ZB(i);
B2(i,5)=B1(i,5)/YB(i);
end
end

for i=1:l
if B1(i,6)==1
ZB(i)=B1(i,14)^2/SB;
YB(i)=1/ZB(i);
B2(i,8)=B2(i,8)/ZB(i);
B2(i,9)=B2(i,9)/ZB(i);
Gt(i)=Gt(i)/YB(i);
Bt(i)=Bt(i)/YB(i);
A1(B1(i,3),6)=Gt(i)-Bt(i)*(1i);
end
end

%%%下面求解节点导纳矩阵
Y=zeros(n);
for i=1:n
if A1(i,6)~=0
Y(i,i)=A1(i,6);
end
end
B3=B2;
for i=1:l
p=B3(i,2);
q=B3(i,3);
if B3(i,6)==1 %%%含有变压器支路支路
Y(p,p)=Y(p,p)+1./((B3(i,8)+B3(i,9)*(1i))*B3(i,7)^2);
Y(p,q)=Y(p,q)-1./((B3(i,8)+B3(i,9)*(1i))*B3(i,7));
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./((B3(i,8)+B3(i,9)*(1i)));
else%%% 无变压器支路
Y(p,p)=Y(p,p)+1./B3(i,4)+B3(i,5)/2;
Y(p,q)=Y(p,q)-1./B3(i,4);
Y(q,p)=Y(p,q);
Y(q,q)=Y(q,q)+1./B3(i,4)+B3(i,5)/2;
end
end
Y

A2=A1;
A2(:,2)=A1(:,2)/SB;
A2(:,3)=A1(:,3)/SB;%功率参数标幺化

G=real(Y);
B=imag(Y);
S=A2(:,2)-A2(:,3) ;
P=real(S);
Q=imag(S);
V=A2(:,4);
delta=A2(:,5);

DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q);
J = jcb( G,B,V,delta,n,B3,A2,sw);
[ ddelta, dV ] = xzl( n,J,DeltaS,A2 );
e=1;
while (max(ddelta)>ac | max(dV)>ac)

a=0;
b=0;
for i=1:n
switch A2(i,7)
case 1
delta(i)=delta(i)+ddelta(i-b);
V(i)=V(i)+dV(i-a-b);
case 2
delta(i)=delta(i)+ddelta(i-b);
V(i)=V(i);
a=a+1;

case 3
delta(i)=delta(i);
V(i)=V(i);
b=b+1;
end
end
DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q);
J = jcb( G,B,V,delta,n,B3,A2,sw);
[ ddelta, dV ] = xzl( n,J,DeltaS,A2 );
e=e+1;
end
V
delta
e
%%%%下面求平衡节点功率
v=V.*cos(delta)+V.*sin(delta)*(1i);
for j=1:n
yu(j)=conj(Y(sw,j))*conj(v(j));
end
S(sw)=sum(yu)*v(sw);
input('平衡节点的功率为');
S(sw)
B2(sw,1)=S(sw);

%%%%下面求解线路功率
for i=1:n
for j=1:n
Sl(i,j)=v(i)*(conj(v(i))*conj(A2(i,6))+conj(v(i)-v(j))*conj(-Y(i,j)));
end
end
input('线路功率为');
Sl

%%%%线路上损耗的功率
for i=1:n
for j=1:n
DertaS1(i,j)=(Sl(i,j)+Sl(j,i))/2;
end
end
input('线路上损耗的功率为');
DertaSz=sum(sum(DertaS1))
4.2 计算功率不平衡量程序
function DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q)
%计算功率不平衡量
for i=1:n
if A2(i,7)~=sw
EP(i)=0;
EQ(i)=0;
for j=1:n
EP(i)=EP(i)+V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
EQ(i)=EQ(i)+V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
end
P1=EP(i);
Q1=EQ(i);

if A2(i,7)==1 %PQ节点
p=2*i-1;
DeltaS(p)=P(i)-P1;
p=p+1;
DeltaS(p)=Q(i)-Q1;
else
p=2*i-1;
DeltaS(p)=P(i)-P1;
p=p+1;
DeltaS(p)=0;
end
end
end
DeltaS(2*sw-1)=[];
DeltaS(2*sw-1)=[];
for i=1:n
if A2(i,7)==2
DeltaS(2*i)=[];
end
end
DeltaS=DeltaS';
End

4.3 计算雅克比矩阵程序
function [ J ] = jcb( G,B,V,delta,n,B3,A2,sw )
%计算雅克比矩阵

for i=1:n
if A2(i,7)==1 %PQ节点
for j=1:n
if j~=i&j~=sw
H=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
J1=-V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
N=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
L=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
p=2*i-1;
q=2*j-1;
J(p,q)=H;
m=p+1;
J(m,q)=J1;
q=q+1;
J(p,q)=N;
J(m,q)=L;
else if j==i&j~=sw
H1=0;
for h=1:n
H1=H1+(-V(i)*V(h)*(G(i,h)*sin(delta(i)-delta(h))-B(i,h)*cos(delta(i)-delta(h))));
end
H=H1+V(i)*V(i)*(-B(i,i));
J2=0;
for h=1:n
J2=J2+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*sin(delta(i)-delta(h))));
end
J1=J2-V(i)*V(i)*G(i,i);
N1=0;
for h=1:n
N1=N1+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*sin(delta(i)-delta(h))));
end
N=N1-V(i)*V(i)*G(i,i)+2*V(i)*V(i)*G(i,i);
L1=0;
for h=1:n
L1=L1+(V(i)*V(h)*(G(i,h)*sin(delta(i)-delta(h))-B(i,h)*cos(delta(i)-delta(h))));
end
L=L1-V(i)*V(i)*(-B(i,i))-2*V(i)*V(i)*B(i,i);
p=2*i-1;
q=2*j-1;
J(p,q)=H;
m=p+1;
J(m,q)=J1;
q=q+1;
J(p,q)=N;
J(m,q)=L;
end
end
end
end
if A2(i,7)==2 %PV节点
for j=1:n
if j~=i&j~=sw
H=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
J1=0;
N=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
L=0;
p=2*i-1;
q=2*j-1;
J(p,q)=H;
m=p+1;
J(m,q)=J1;
q=q+1;
J(p,q)=N;
J(m,q)=L;
else if j==i&j~=sw
H1=0;
for h=1:n
H1=H1+(-V(i)*V(h)*(G(i,h)*sin(delta(i)-delta(h))-B(i,h)*cos(delta(i)-delta(h))));
end
H=H1+V(i)*V(i)*(-B(i,i));
J1=0;
N1=0;
for h=1:n
N1=N1+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*sin(delta(i)-delta(h))));
end
N=N1-V(i)*V(i)*G(i,i)+2*V(i)*V(i)*G(i,i);
L=0;
p=2*i-1;
q=2*j-1;
J(p,q)=H;
m=p+1;
J(m,q)=J1;
q=q+1;
J(p,q)=N;
J(m,q)=L;
end
end
end
end
end

a=0;
for i=1:n
switch A2(i,7)
case 1

case 2
a=a+1;
J(2*i,:)=[];
J(:,2*i)=[];

case 3
J(2*i-1-a,:)=[];
J(2*i-1-a,:)=[];
J(:,2*i-1-a)=[];
J(:,2*i-1-a)=[];
end
end

end

4.4 计算电压修正量的程序
function [ ddelta, dV ] = xzl( n,J,DeltaS ,A2)
%求取电压不平衡量
dX=J\DeltaS;
a=2;
for i=1:n
switch A2(i,7)

case 1
ddelta(i)=dX(a-1);
dV(i)=dX(a);
a=a-0;
a=a+2;

case 2
ddelta(i)=dX(a-1);
dV(i)=0;
a=a-1;
a=a+2;

case 3
ddelta(i)=0;
dV(i)=0;
a=a-2;
a=a+2;
end
end
a=0;
for i=1:n
switch A2(i,7)

case 1

case 2
a=a+1;
dV(i)=[];

case 3
ddelta(i)=[];
dV(i-a)=[];
end
end

end

5 程序调试
首先，在MATLAB命令窗口输入主函数名“shuangraozu”运行程序，结果如下：

进而根据提示输入相应参数如下：

输入支路参数矩阵如下：
[1 2 1 0 0 1 1 665.2 14 189.6 0.299 15.75 246 15.75;
2 2 3 2.12+17.2i 1.176e-4i 0 0 0 0 0 0 0 0 220;
3 2 6 2.65+21.5i 1.47e-4i 0 0 0 0 0 0 0 0 220;
4 3 6 1.59+12.9i 8.82e-5i 0 0 0 0 0 0 0 0 220;
5 3 7 0 0 1 1 395.65 13.85 74.75 0.0444 10.5 180 10.5;
6 4 7 0 0 1 1 243.25 -0.95 0 0 10.5 180 10.5 ;
7 5 7 0 0 1 1 256.75 9.15 0 0 10.5 180 10.5]
随后输入节点参数矩阵如下：
[1 150 0 1 0 0 2;
2 0 100+100i 1 0 0 1;
3 0 0 1 0 0 1;
4 0 100+80i 1 0 0 1;
5 0 20+15i 1 0 0 1;
6 0 0 1 0 0 3;
7 0 0 1 0 0 1]
程序运行计算可得节点导纳矩阵如下：
可见与手算潮流例题中所得的节点导纳矩阵基本一致。
计算初次功率不平衡量矩阵如下：

计算第一次迭代的雅克比矩阵如下：

计算第一次迭代的电压修正量矩阵如下：

通过收敛条件判断电压修正量是否满足精度要求，若满足，则输出结果，若不满足，继续迭代。
迭代4次后，满足精度要求，输入电压结果如下：

可见，机算潮流迭代4次后得到的电压结果与手算潮流迭代两次所得的潮流计算结果虽有一定误差，但是收敛的方向和趋势基本相同，进而验证了程序的正确性。
19