牛拉法潮流计算.docx
《牛拉法潮流计算.docx》由会员分享,可在线阅读,更多相关《牛拉法潮流计算.docx(17页珍藏版)》请在冰豆网上搜索。
牛拉法潮流计算
%本程序的功能是用牛拉法进行潮流计算
%原理介绍详见鞠平著《电气工程》
%默认数据为鞠平著《电气工程》例8.4所示数据
%B1是支路参数矩阵
%第一列和第二列是节点编号。
节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点编号%第三列为支路的串列阻抗参数,含变压器支路此值为变压器短路电抗%第四列为支路的对地导纳参数,含变压器支路此值不代入计算%第五烈为含变压器支路的变压器的变比,变压器非标准电压比
%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,“0”为不含有变压器
%B2为节点参数矩阵
%第一列为节点注入发电功率参数%第二列为节点负荷功率参数%第三列为节点电压参数
%第四列
%第五列%第六列为节点类型参数,“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数
%X为节点号和对地参数矩阵
%第一列为节点编号
%第二列为节点对地参数
clear;
clc;
num=input('是否采用默认数据?
(1-默认数据;2-手动输入)');
ifnum==1
n=4;
n1=4;
isb=4;
pr=0.00001;
B1=[120.1667i00.88641;130.1302+0.2479i0.0258i10;140.1736+0.3306i
0.0344i10;340.2603+0.4959i0.0518i10];
B2=[001002;0-0.5-0.3i1002;0.201.05003;0-0.15-0.1i1.05001];
X=[10;20.05i;30;40];
else
n=input('请输入节点数:
n=');n1=input('请输入支路数:
n1=');isb=input('请输入平衡节点号:
isb=');pr=input('请输入误差精度:
pr=');B1=input('请输入支路参数:
B1=');B2=input('请输入节点参数:
B2=');X=input('节点号和对地参数:
X=');
end
Times=1;%迭代次数
%创建节点导纳矩阵
Y=zeros(n);
fori=1:
n1
ifB1(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)-B1(i,5)/B1(i,3);
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3);
end
end
fori=1:
n1
Y(i,i)=Y(i,i)+X(i,2);%计及补偿电容电纳
end
disp('导纳矩阵为:
');
disp(Y);%显示导纳矩阵
%初始化OrgS、DetaSOrgS=zeros(2*n-2,1);
DetaS=zeros(2*n-2,1);
%创建OrgS,用于存储初始功率参数
h=0;
j=0;
fori=1:
n%对PQ节点的处理
ifi~=isb&B2(i,6)==2%不是平衡点&是PQ点
h=h+1;
forj=1:
n
%公式8-74
%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)
%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)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)));
j,3))-
OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(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
fori=1:
n%对PV节点的处理,注意这时不可再将h初始化为0
ifi~=isb&B2(i,6)==3%不是平衡点&是PV点
h=h+1;
forj=1:
n
%公式8-75-a
%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)
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
%创建PVU用于存储PV节点的初始电压
PVU=zeros(n-h-1,1);
t=0;
fori=1:
n
ifB2(i,6)==3
t=t+1;
PVU(t,1)=B2(i,3);
end
end
%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;
fori=1:
n%对PQ节点的处理
ifi~=isb&B2(i,6)==2
h=h+1;
DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);%delPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);%delQiend
end
t=0;
fori=1:
n%对PV节点的处理,注意这时不可再将h初始化为0
ifi~=isb&B2(i,6)==3
h=h+1;
t=t+1;
DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);%delPi
DetaS(2*h,1)=real(PVU(t,1))A2+imag(PVU(t,1))A2-real(B2(i,3))A2-imag(B2(i,3))A2;%delUiend
end
%DetaS
%创建I,用于存储节点电流参数
i=zeros(n-1,1);
h=0;
fori=1:
n
ifi~=isb
h=h+1;
I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭
end
end
%创建Jacbi(雅可比矩阵)
Jacbi=zeros(2*n-2);
h=0;
k=0;
fori=1:
n%对PQ节点的处理
ifB2(i,6)==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==(n-1)%将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;
end
end
end
end
end
k=0;
fori=1:
n%对PV节点的处理
ifB2(i,6)==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==(n-1)%将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;
end
end
end
end
end
disp('初始雅可比矩阵为:
');disp(Jacbi);
%求解修正方程,获取节点电压的不平衡量
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;%inv矩阵求逆
%DetaU
%修正节点电压
j=0;
fori=1:
n%对PQ节点处理
ifB2(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
fori=1:
n%对PV节点的处理
ifB2(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
%开始循环**********************************************************************whileabs(max(DetaU))>pr
OrgS=zeros(2*n-2,1);
h=0;
j=0;
fori=1:
n
ifi~=isb&B2(i,6)==2h=h+1;
forj=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(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
fori=1:
n
ifi~=isb&B2(i,6)==3
h=h+1;
forj=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(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;
fori=1:
n
ifi~=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;
fori=1:
n
j,3))-
j,3))-
ifi~=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,1)=real(B2(i,1))-OrgS(2*h-1,1);
DetaS(2*h,1)=real(PVU(t,1))A2+imag(PVU(t,1))A2-real(B2(i,3))A2-imag(B2(i,3)F2;
end
end
%DetaS
%创建I
i=zeros(n-1,1);
h=0;
fori=1:
n
ifi~=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;
fori=1:
n
ifB2(i,6)==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==(n-1)
k=0;
end
end
end
end
end
k=0;
fori=1:
n
ifB2(i,6)==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==(n-1)
k=0;
end
end
end
end
end
%Jacbi
DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS;
%DetaU
%修正节点电压
j=0;
fori=1:
n
ifB2(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
fori=1:
n
ifB2(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
disp('迭代次数为:
');
disp(Times);
disp('收敛时电压修正量为:
:
');disp(DetaU);
fork=1:
n
E(k)=B2(k,3);e(k)=real(E(k));
f(k)=imag(E(k));
V(k)=sqrt(e(kF2+f(k)A2);
sida(k)=atan(f(k)./e(k))*180./pi;
end
%===============计算各输出量===========================disp('各节点的实际电压标幺值E为(节点号从小到大排列):
');
disp(E);%显示各节点的实际电压标幺值E用复数表示
disp('');
disp('各节点的电压大小V为(节点号从小到大排列):
');
disp(V);%显示各节点的电压大小V的模值
disp('');
disp('各节点的电压相角sida为(节点号从小到大排列):
');disp(sida);%显示各节点的电压相角
forp=1:
n
C(p)=0;
forq=1:
n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);%计算各节点的功率S=电压X注入电流的共轭值enddisp('各节点的功率S为(节点号从小到大排列):
');
disp(S);%显示各节点的注入功率
Sline=zeros(n1,5);
disp('');
disp('各条支路的首端功率Si为(顺序同您输入B1时一致):
');fori=1:
n1
p=B1(i,1);
q=B1(i,2);
Sline(i,1)=B1(i,1);
Sline(i,2)=B1(i,2);
ifB1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3)));
Siz(i)=Si(p,q);
end
SSi(p,q)=Si(p,q);
Sline(i,3)=Siz(i);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);
end
disp('');
disp('各条支路的末端功率Sj为(顺序同您输入B1时一致):
');
fori=1:
n1
p=B1(i,1);q=B1(i,2);
ifB1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3)));
Sjy(i)=Sj(q,p);
end
SSj(q,p)=Sj(q,p);
Sline(i,4)=Sjy(i);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);
end
disp('');
disp('各条支路的功率损耗DS为(顺序同您输入B1时一致):
');fori=1:
n1
p=B1(i,1);
q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
DDS(i)=DS(