基于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(pqn ode(1,m),t)*e(1,pqnode(1,m));
N=G(pqnode(1,m),t)*e(1,pqnode(1,m))+B(pq node(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(pvno de(1,m),t)*e(1,pvnode(1,m));
N=G(pvnode(1,m),t)*e(1,pvnode(1,m))+B(pvn ode(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
aii=real(I);
bii=imag(I);
H=-B(t,t)*e(1,pvnode(1,m))+G(t,t)*f(1,pvnode (1,m))+bii;
N=G(t,t)*e(1,pvnode(1,m))+B(t,t)*f(1,pvnode( 1,m))+aii;
R=2*f(1,pvnode(1,m));
S=2*e(1,pvnode(1,m));
end
Jacob(n,k)=H; k=k+1;
Jacob(n,k)=N; k=k-1;n=n+1;
Jacob(n,k)=R; k=k+1;
Jacob(n,k)=S;
n=n-1;k=k+1; %按照雅克比矩阵的排列规则排列PV节点的雅克比元素end
k=1;
n=n+2; %定位于下一个待排列PV节点的雅克比元素第一个位置
end
%*************************电压变化量化的计算与存储************************************ Detau=inv(Jacob)*deta; %构建电压的变化量的列向量
f=zeros(1,nnum); %给电压实部赋初值0
e=zeros(1,nnum); %给电压虚部赋初值0
for p=1:(pq+pv);
f(1,cerate(p,1))=j*Detau(2*p-1,1);%将电压变量的奇数行赋值给f
e(1,cerate(p,1))=Detau(2*p,1); %将电压变量的偶数行赋值给e
end
t=e+f;
xiumax=abs(Detau(1,1)); %将电压变化量的第一个元素赋值给最大允许误差
for n=2:2*nnum-2;
if abs(Detau(n,1))>xiumax
xiumax=abs(Detau(n,1)); %找出最大的电压误差
end
end
Uoriginal=Uoriginal+t; %迭代修正后的电压值
counter=1+counter; %统计迭代次数
end
disp('迭代次数counter:');
disp(counter);
%**************************平衡节点功率及显示**********************************
m=blancenode;
t=0;
for n=1:nnum;
t=t+(G(m,n)-j*B(m,n))*(real(Uoriginal(1,n))-j*i mag(Uoriginal(1,n)));
end
S1(1,m)=Uoriginal(1,m)*t;
%**************************直角坐标下各节点电压及显示****************************
U=zeros(1,nnum);
for n=1:nnum
Ui(n,1)=Node(n,1);
U(1,n)=real(Uoriginal(1,n))+i*imag(Uoriginal(1 ,n)); %将电压值由极坐标转化为直角坐标形式
Ui(n,2)=U(1,n);
Ui(n,3)=S1(1,n);
end
disp('各节点电压直角坐标形式及节点注入功率:');
disp(' 节点号节点电压值节点注入功率');
disp(Ui);
disp('修正电压的最大误差:')
disp(xiumax);
%**************************功率损耗************************************* **
for m=1:bnum %支路功率及损耗
startnode=branch(m,1);
endnode=branch(m,2); %终止节点
y=sum(Y,2);
Sij=Uoriginal(1,startnode)*((conj(Uoriginal(1,s tartnode))*conj(y(startnode,1)))+(conj(Uorigi nal(1,startnode))-conj(Uoriginal(1,endnode))) *conj(-Y(startnode,endnode)));
Sji=Uoriginal(1,endnode)*((conj(Uoriginal(1,e ndnode))*conj(y(endnode,1)))+(conj(Uoriginal (1,endnode))-conj(Uoriginal(1,startnode)))*co nj(-Y(endnode,startnode)));
S(m,1)=startnode; %始节点
S(m,2)=endnode; %末节点
S(m,3)=Sij; %始节点流入的功率
S(m,4)=Sji; %末节点流入的功率
S(m,5)=Sij+Sji; %线路损耗
end
disp('支路功率及损耗：');
disp('节点号(I) 节点号(J) 支路功率（I-J）支路功率（J-I）线路损耗（delta_S)')
disp(S); 最新文件---------------- 仅供参考--------------------已改成word文本--------------------- 方便更改。