%本程序的功能是用牛拉法进行潮流计算 %原理介绍详见鞠平著《电气工程》
%默认数据为鞠平著《电气工程》例8.4所示数据
±是支路参数矩阵
%第一列和第二列是节点编号。节点编号由小到大编写
%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号 %第三列为支路的串列阻抗参数，含变压器支路此值为变压器短路电抗 %第四列为支路的对地导纳参数，含变压器支路此值不代入计算 %第五烈为含变压器支路的变压器的变比，变压器非标准电压比
%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器
2为节点参数矩阵
%第一列为节点注入发电功率参数 %第二列为节点负荷功率参数 %第三列为节点电压参数 %第四列 %第五列
%第六列为节点类型参数，“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数
%X为节点号和对地参数矩阵 %第一列为节点编号 %第二列为节点对地参数
clear; clc;
num=input('是否采用默认数据？(1-默认数据;2-手动输入)'); if num==1 n=4; n1=4; isb=4;
pr=0.00001;
B1=[1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0518i 1 0];
B2=[0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1]; X=[1 0;2 0.05i;3 0;4 0]; 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); 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)-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
for i=1:n1
Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳 end
disp('导纳矩阵为：');
disp(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 %不是平衡点&是PQ点 h=h+1; for j=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))); 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 %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3 %不是平衡点&是PV点 h=h+1; for j=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; for i=1:n
if B2(i,6)==3 t=t+1;
PVU(t,1)=B2(i,3); end end
%创建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); TlPi DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); TlQi 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,1))-OrgS(2*h-1,1); TlPi
DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; TlUi 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));%conj求共轭 end end
%创建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
disp('初始雅可比矩阵为：'); disp(Jacbi);
%求解修正方程，获取节点电压的不平衡量 DetaU=zeros(2*n-2,1);
DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆 % 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;