基于MATLAB的潮流计算源程序代码

%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算********** clear clc

load E:\\data\\IEEE014_Node.txt Node=IEEE014_Node; weishu=size(Node);

nnum=weishu(1,1); %节点总数

load E:\\data\\IEEE014_Branch.txt branch=IEEE014_Branch;

bwei=size(branch); bnum=bwei(1,1); %支路总数

Y=(zeros(nnum)); Sj=100;

%********************************节点导纳矩阵******************************* for m=1:bnum;

s=branch(m,1); %首节点

e=branch(m,2); %末节点

R=branch(m,3); %支路电阻

X=branch(m,4); %支路电抗

B=branch(m,5); %支路对地电纳

k=branch(m,6);

if k==0 %无变压器支路情形

Y(s,e)=-1/(R+j*X); %互导纳

Y(e,s)=Y(s,e); end

if k~=0 %有变压器支路情形

Y(s,e)=-(1/((R+j*X)*k)); Y(e,s)=Y(s,e);

Y(s,s)=-(1-k)/((R+j*X)*k^2); Y(e,e)=-(k-1)/((R+j*X)*k); %对

地导纳

end

Y(s,s)=Y(s,s)-j*B/2;

Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形 end

for t=1:nnum;

Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13); %求支路自导纳 end

G=real(Y); %电导 B=imag(Y); %电纳

%******************节点分类**************************************

pq=0; pv=0; blancenode=0; pqnode=zeros(1,nnum); pvnode=zeros(1,nnum); for m=1:nnum;

if Node(m,2)==3

blancenode=m; %平衡节点编号

else if Node(m,2)==0 pq=pq+1;

pqnode(1,pq)=m; %PQ节点编号

else if Node(m,2)==2 pv=pv+1;

pvnode(1,pv)=m; %PV节点编号

end end end end

%*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化 for n=1:nnum

Uoriginal(1,n)=Node(n,9); %对各点电压赋初值

if Node(n,9)==0;

Uoriginal(1,n)=1; %该节点为非PV节点时，将电压值赋为1 end end

Presion=input('请输入误差精度要求：Presion=');

disp('该电力系统节点数:'); disp(nnum); xiumax=0.1; counter=0;

while xiumax>Presion

%****************************计算不平衡量*********************************** e=real(Uoriginal); %取初始电压的实部

f=imag(Uoriginal); %取初始电压的虚部

deta=zeros(2*pq+2*pv,1); %构造储存功率变化量的列矩阵 n=1; for m=1:pq; Pi=0;Qi=0;

for t=1:nnum;

Pi=Pi+e(1,pqnode(1,m))*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t))+f(1,pqnode(1,m))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t));%计算该PQ节点的负荷有功

Qi=Qi+f(1,pqnode(1,m))*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1,m))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t));%计算该PQ节点的负荷无功 end

S1(1,pqnode(1,m))=Pi+j*Qi;

P=(Node(pqnode(1,m),7)-Node(pqnode(1,m),5))/Sj-Pi;%计算该PQ节点的实际有功功率 deta(n,1)=P;%在该列向量中储存有功功率

n=n+1;

Q=(Node(pqnode(1,m),8)-Node(pqnode(1,m),

6))/Sj-Qi;%计算该PQ节点的实际无功功率 deta(n,1)=Q;%在该列向量中储存无功功率

n=n+1; end

for m=1:pv; Pv=0; Qv=0; for t=1:nnum;

Pv=Pv+e(1,pvnode(1,m))*(G(pvnode(1,m),t)*e(1,t)-B(pvnode(1,m),t)*f(1,t))+f(1,pvnode(1,m))*(G(pvnode(1,m),t)*f(1,t)+B(pvnode(1,m),t)*e(1,t));%计算该PV节点的负荷有功

Ui=e(1,pvnode(1,m))^2+f(1,pvnode(1,m))^2;%计算该节点的负荷电压值

Qv=Qv+f(1,pqnode(1,m))*(G(pqnode(1,m),t)*e(1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1,m))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t)*e(1,t)); end

S1(1,pvnode(1,m))=Pv+j*Qv;

P=(Node(pvnode(1,m),7)-Node(pvnode(1,m),5))/Sj-Pv; %计算该节点的实际有功功率 deta(n,1)=P; %储存该有功功率 n=n+1;

U=Node(pvnode(1,m),3)^2-Ui; %计算电压变化量

deta(n,1)=U; %储存该电压变化量 n=n+1; end deta;

cerate=zeros(pq+pv,1); for k=1:pq

cerate(k,1)=pqnode(1,k); end

for v=1:pv

cerate(pq+v,1)=pvnode(1,v);

end

%******************************雅克比矩阵****************************** Jacob=ones(2*nnum-2);

L=0;J=0;H=0;N=0; R=0;S=0;n=1;k=1; for

m=1:pq; %m表示雅克比矩阵中pq节点的行数

for u=1:pq+pv; %u表示雅克比矩阵中pq节点的列数

t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的个数

if pqnode(1,m)~=t %非对角元素的情况

H=G(pqnode(1,m),t)*f(1,pqnode(1,m))-B(pqnode(1,m),t)*e(1,pqnode(1,m));

N=G(pqnode(1,m),t)*e(1,pqnode(1,m))+B(pqnode(1,m),t)*f(1,pqnode(1,m)); L=H; J=-N; else if

pqnode(1,m)==t %对角线元素时的情况

I=0;

for g=1:nnum

I=Y(t,g)*Uoriginal(1,g)+I; %计算节点的注入电流

end aii=real(I); bii=imag(I);

H=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode(1,m))+bii;

N=G(t,t)*e(1,pqnode(1,m))+B(t,t)*f(1,pqnode(1,m))+aii;

L=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode(1,m))-bii;

J=-G(t,t)*e(1,pqnode(1,m))-B(t,t)*f(1,pqnode(1,m))+aii;

end end

Jacob(n,k)=H; k=k+1;

Jacob(n,k)=N;

k=k-1;n=n+1; Jacob(n,k)=J; k=k+1;

Jacob(n,k)=L;

n=n-1; k=k+1; %按照雅克比矩阵的排列规则排列pq节点的雅克比元素 end

k=1; n=2*m+1; %将光标定位于下一个待排列PQ节点元素的第一个位置 end

n=2*pq+1; k=1; %定位于PV节点的第一个位置处 for m=1:pv;

for u=1:pq+pv;

t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的位置

if pvnode(1,m)~=t %非对角线元素情况

H=G(pvnode(1,m),t)*f(1,pvnode(1,m))-B(pvnode(1,m),t)*e(1,pvnode(1,m));

N=G(pvnode(1,m),t)*e(1,pvnode(1,m))+B(pvnode(1,m),t)*f(1,pvnode(1,m)); R=0; S=0; end

if pvnode(1,m)==t %对角线元素情况 I=0;

for g=1:nnum

I=Y(t,g)*Uoriginal(1,g)+I; %计算PV节点的注入电流 end