潮流计算程序参考

基于MATLAB的电力系统潮流计算

%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：

%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号，将变压器的串联阻抗置于低压侧处理。

%第三列为支路的串列阻抗参数。

%第四列为支路的对地导纳参数。

%第五烈为含变压器支路的变压器的变比

%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。

%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中

%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地%参数。

n=input('请输入节点数:n=');

n1=input('请输入支路数:n1=');

isb=input('请输入平衡节点号:isb=');

pr=input('请输入误差精度:pr=');

B1=input('请输入支路参数:B1=');

B2=input('请输入节点参数:B2=');

X=input('节点号和对地参数:X=');

Y=zeros(n);

Times=1; %置迭代次数为初始值

%创建节点导纳矩阵

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)-1/(B1(i,3)*B1(i,5));

Y(q,p)=Y(p,q);

Y(p,p)=Y(p,p)+1/B1(i,3);

Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));

end

end

OrgS=zeros(2*n-2,1);

DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化

%创建OrgS，用于存储初始功率参数

h=0;

j=0;

for i=1:n %对PQ节点的处理

if i~=isb&B2(i,6)==2

h=h+1;

for j=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)))+im ag(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(B 2(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

h=h+1;

for j=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)))+im ag(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(B 2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));

end

end

end

OrgS

%创建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

PVU

%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量

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);

DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);

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,2))-OrgS(2*h-1,1);

DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;

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));

end

end

I

%创建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));