潮流计算matlab牛顿拉夫逊法法.docx
《潮流计算matlab牛顿拉夫逊法法.docx》由会员分享，可在线阅读，更多相关《潮流计算matlab牛顿拉夫逊法法.docx（36页珍藏版）》请在冰豆网上搜索。
潮流计算matlab牛顿拉夫逊法法
用matlab潮流计算（牛顿拉夫逊法）
%主程序：
[d]=uigetfile（'ieee14.m','SelectDataFile'）;
ifpathname==0
error（'youmustselectavaliddatafile'）
else
l（dfile）;
%stripoff.m
eval（d））;%打开数据文件
end
globaln;
globalm;
[nb,mb]=size（bus）;%节点重新编号
[nl,ml]=size（line）;
nSW=0;%平衡节点数目
nPV=0;%PV节点数目
nPQ=0;%PQ节点数目
fori=1:
nb,%nb为总节点数
type=bus（i,6）;
iftype==3,
nSW=nSW+1;%统计平衡节点数目
SW（nSW,:
）=bus（i,:
）;
elseiftype==2,
nPV=nPV+1;%统计PV节点数目
PV（nPV,:
）=bus（i,:
）;
else
nPQ=nPQ+1;%统计PQ节点数目
PQ（nPQ,:
）=bus（i,:
）;
end
end;
bus=[PQ;PV;SW];
newbus=[1:
nb]';
f=bus（:
1）;
nodenum=[newbusbus（:
1）];
bus（:
1）=newbus;
fori=1:
nl
forj=1:
2
fork=1:
nb
ifline（i,j）==nodenum（k,2）
line（i,j）=nodenum（k,1）;
break
end
end
end
end
Y=y（bus,line）;%形成节点导纳矩阵
迭代次数初值%K=0;
Kmax=10;%最大迭代次数
eps1=1.0e-10;
eps2=1.0e-10;
m=nPQ;n=nb;
Um=eye（m,m）;
myf=fopen（'output1.dat','w'）;
forK=1:
Kmax
fori=1:
m
forj=1:
m
ifi==j
Um（i,j）=bus（i,2）;
end
end
end
b=dPQ（Y,bus）;
C=jac（bus,Y）;
dX=C\b';
dx=dX';
[nx,mx]=size（dx）;
fori=1:
n-1%计算相角
bus（i,3）=bus（i,3）-dX（i,1）;
end
B=dx（nx,n:
mx）*Um;%计算电压差
bus（1:
m,2）=bus（1:
m,2）-B';%计算电压值
dx（nx,n:
mx）=B;
fprintf（myf,'--第%d次迭代时雅可比矩阵--',K）
fprintf（myf,'\n'）;
fori=1:
（n+m-1）
forj=1:
（n+m-1）
fprintf（myf,'%8.6f',C（i,j））;
fprintf（myf,''）;
end
fprintf（myf,'\n'）;
end
fprintf（myf,'--第%d次迭代时dPQ的误差--',K）
fprintf（myf,'\n'）;
fori=1:
（n+m-1）
fprintf（myf,'%8.6e',b（1,i））;
fprintf（myf,'\n'）;
end
fprintf（myf,'\n'）;
fprintf（myf,'--第%d次迭代时dx（误差）--',K）
fprintf（myf,'\n'）;
fori=1:
（n+m-1）
fprintf（myf,'%8.6e',dX（i,1））;
fprintf（myf,'\n'）;
end
fprintf（myf,'\n'）;
fprintf（myf,'第%d次迭代后节点电压（仅PQ节点）',K）
fprintf（myf,'\n'）;
fori=1:
m
fprintf（myf,'%8.6f',bus（i,2））;
fprintf（myf,'\n'）;
end
fprintf（myf,'\n'）;
fprintf（myf,'第%d次迭代后相角（角度）',K）
fprintf（myf,'\n'）;
fori=1:
n
fprintf（myf,'%8.6f',bus（i,3）*180/pi）;
fprintf（myf,'\n'）;
end
fprintf（myf,'\n'）;
if（max（abs（dx））break;
end
end
%计算功率
P1=0;T=0;%计算平衡节点的有功
forj=1:
n
T=bus（n,2）*bus（j,2）*（real（Y（n,j））*cos（bus（n,3）-bus（j,3））+imag（Y（n,j））*sin（bus（n,3）-bus（j,3）））;
P1=P1+T;
end
bus（n,4）=P1;
fork=m+1:
n%计算PV节点、平衡节点的无功
Q1=0;T=0;
forj=1:
n
T=bus（k,2）*bus（j,2）*（real（Y（k,j））*sin（bus（k,3）-bus（j,3））-imag（Y（k,j））*cos（bus（k,3）-bus（j,3）））;
Q1=Q1+T;
end
bus（k,5）=Q1;
end
bus（:
1）=f;%换回各节点、支路的初始编号
r=zeros（1,mb）;
fort=1:
n
forl=t+1:
n
ifbus（t,1）>bus（l,1）
r=bus（t,:
）;bus（t,:
）=bus（l,:
）;bus（l,:
）=r;
end
end
end
fori=1:
nl
forj=1:
2
fork=1:
nb
ifline（i,j）==nodenum（k,1）
line（i,j）=nodenum（k,2）;
break
end
end
end
end
fclose（myf）;
Pf=loss（bus,line）;%计算支路潮流及损耗
%将节点导纳矩阵、节点潮流计算结果写入文件output2
myf=fopen（'output2.dat','w'）;
fprintf（myf,'--节点导纳矩阵--\n'）;
fork=1:
n
forj=1:
n
fprintf（myf,'%8.6f',real（Y（k,j）））;
fprintf（myf,'+i*'）;
fprintf（myf,'%8.6f',imag（Y（k,j）））;
fprintf（myf,''）;
end
fprintf（myf,'\n'）;
end
fprintf（myf,'------------------牛顿－拉夫逊法潮流计算结果-------------\n'）;
fprintf（myf,'------------节点计算结果------------\n'）;
fprintf（myf,'--节点节点电压节点相角注入有功功率（P）注入无功功率（Q）类型--\n'）;
forl=1:
nb
forj=1:
mb
ifj==1|j==6
',bus（l,j））;fprintf（myf,'%8.1f
elseifj==3
',bus（l,j）*180/pi）;fprintf（myf,'%8.6f
else
',bus（l,j））;fprintf（myf,'%8.6f
end
end
'\n'）;fprintf（myf,
end
--\n'）;
支路计算结果fprintf（myf,'--
S（J,I）线路损耗线路功率S（I,J）线路功率节点（I）fprintf（myf,'--节点（J）
dS（I,J）--\n'）;
fork=1:
nl
forj=1:
5
ifj<=2
',Pf（k,j））;fprintf（myf,'%8.1f
'）;fprintf（myf,'
else
fprintf（myf,'%8.6f',real（Pf（k,j）））;
fprintf（myf,'+i*'）;
fprintf（myf,'%8.6f',imag（Pf（k,j）））;
'）;fprintf（myf,'
end
end
fprintf（myf,'\n'）;
end
fclose（myf）;
%根据支路参数建立节点导纳矩阵程序：
functionY=y（bus,line）
%目的：
根据支路参数建立节点导纳矩阵
%输入：
节点参数矩阵--bus；支路参数矩阵--line
%输出：
节点导纳矩阵--Y
[nb,mb]=size（bus）;
[nl,ml]=size（line）;
Y=zeros（nb,nb）;
fork=1:
nl
I=line（k,1）;
J=line（k,2）;
Zt=line（k,3）+j*line（k,4）;
ifZt==0
disp（'对地支路'）;
Yt=inf;
else
Yt=1/Zt;
end
Ym=line（k,5）+j*line（k,6）;
K=line（k,7）;
if（K==0）&（J~=0）%普通线路
Y（I,I）=Y（I,I）+Yt+Ym;
Y（J,J）=Y（J,J）+Yt+Ym;
Y（I,J）=Y（I,J）-Yt;
Y（J,I）=Y（I,J）;
end
if（K==0）&（J==0）%对地支路K=J=0，R=X=0
Y（I,I）=Y（I,I）+Ym;
end
ifK>0%K>0时变压器支路
Y（I,I）=Y（I,I）+Yt+Ym;
Y（J,J）=Y（J,J）+Yt/K^2;
Y（I,J）=Y（I,J）-Yt/K;
Y（J,I）=Y（I,J）;
end
ifK<0%K<0时变压器支路
Y（I,I）=Y（I,I）+Yt+Ym;
Y（J,J）=Y（J,J）+Yt*K^2;
Y（I,J）=Y（I,J）+Yt*K;
Y（J,I）=Y（I,J）;
end
end
%形成雅克矩阵程序：
functionJ=jac（bus,Y）
%形成雅可比矩阵
globaln;
globalm;
fori=1:
（n-1）%计算PQ、PV节点的有功
P1=0;T=0;
forj=1:
n
T=bus（i,2）*bus（j,2）*（real（Y（i,j））*cos（bus（i,3）-bus（j,3））+imag（Y（i,j））*sin（bus（i,3）-bus（j,3）））;
P1=P1+T;
end
bus（i,4）=P1;
end
fori=1:
n-1%计算PV、PQ节点的无功
Q1=0;T=0;
forj=1:
n
T=bus（i,2）*bus（j,2）*（real（Y（i,j））*sin（bus（i,3）-bus（j,3））-imag（Y（i,j））*cos（bus（i,3）-bus（j,3）））;
Q1=Q1+T;
end
bus（i,5）=Q1;
end
fori=1:
n-1%计算H
forj=1:
n-1
ifi~=j
H（i,j）=-bus（i,2）*bus（j,2）*（real（Y（i,j））*sin（bus（i,3）-bus（j,3））-imag（Y（i,j））*cos（bus（i,3）-bus（j,3）））;
N（i,j）=-bus（i,2）*bus（j,2）*（real（Y（i,j））*cos（bus（i,3）-bus（j,3））+imag（Y（i,j））*sin（bus（i,3）-bus（j,3）））;
K（i,j）=bus（i,2）*bus（j,2）*（real（Y（i,j））*cos（bus（i,3）-bus（j,3））+imag（Y（i,j））*sin（bus（i,3）-bus（j,3）））;
L（i,j）=H（i,j）;
end
end
end
fori=1:
n-1%计算H
forj=1:
n-1
ifj==i
H（i,i）=bus（i,5）+imag（Y（i,i））*bus（i,2）^2;
N（i,i）=-bus（i,4）-real（Y（i,i））*bus（i,2）^2;
K（i,i）=N（i,i）+2*real（Y（i,i））*bus（i,2）^2;
L（i,i）=-H（i,i）+2*imag（Y（i,i））*bus（i,2）^2;
end
end
end
N=N（1:
n-1,1:
m）;
K=K（1:
m,1:
n-1）;
L=L（1:
m,1:
m）;
J=[H,N;K,L];
%计算dPQ的程序：
functiondPQ=dPQ（Y,bus）
TlP--有功偏差量
TlQ--无功偏差量
%Y--节点导纳矩阵
%bus--节点参数（P,Q,U及相角）矩阵
globaln;
globalm;
delP=zeros（1,n-1）;
delQ=zeros（1,m）;
fori=1:
（n-1）%形成delP矩阵（PQ、PV节点）
P1=0;T=0;
forj=1:
n
T=bus（i,2）*bus（j,2）*（real（Y（i,j））*cos（bus（i,3）-bus（j,3））+imag（Y（i,j））*sin（bus（i,3）-bus（j,3）））;
P1=P1+T;
end
delP（1,i）=bus（i,4）-P1;
end
fori=1:
m%形成delQ矩阵（PQ节点）
Q1=0;T=0;
forj=1:
n
T=bus（i,2）*bus（j,2）*（real（Y（i,j））*sin（bus（i,3）-bus（j,3））-imag（Y（i,j））*cos（bus（i,3）-bus（j,3）））;
Q1=Q1+T;
end
delQ（1,i）=bus（i,5）-Q1;
end
dPQ=[delP,delQ];
%计算线路损耗、线路潮流程序：
functionPf=loss（bus,line）
%计算线路损耗、线路潮流
[nl,ml]=size（line）;
Pf=zeros（nl,5）;
fork=1:
nl
I=line（k,1）;
J=line（k,2）;
Zt=line（k,3）+i*line（k,4）;
ifZt==0
Yt=inf;
else
Yt=1/Zt;
end
Ym=line（k,5）+i*line（k,6）;
K=line（k,7）;
普通线路潮流%if（K==0）&（J~=0）
S（I,J）=bus（I,2）^2*（conj（Yt）+conj（Ym））-bus（I,2）*（cos（bus（I,3））+i*sin（bus（I,3）））*bus（J,2）*（cos（bus（J,3））-i*sin（bus（J,3）））*conj（Yt）;
S（J,I）=bus（J,2）^2*（conj（Yt）+conj（Ym））-bus（J,2）*（cos（bus（J,3））+i*sin（bus（J,3）））*bus（I,2）*（cos（bus（I,3））-i*sin（bus（I,3）））*conj（Yt）;
delS（I,J）=S（I,J）+S（J,I）;
end
if（K==0）&（J==0）%对地支路潮流
J=5;
S（I,5）=bus（I,2）^2*conj（Ym）;
end
ifK>0%变压器支路k>0时的潮流
S（I,J）=bus（I,2）^2*（conj（Ym+Yt*（1-1/K））+conj（Yt/K））-bus（I,2）*（cos（bus（I,3））+i*sin（bus（I,3）））*bus（J,2）*（cos（bus（J,3））-i*sin（bus（J,3）））*conj（Yt/K）;
S（J,I）=bus（J,2）^2*（conj（Yt））/K^2-bus（J,2）*（cos（bus（J,3））+i*sin（bus（J,3）））*bus（I,2）*（cos（bus（I,3））-i*sin（bus（I,3）））*conj（Yt/K）;
delS（I,J）=S（I,J）+S（J,I）;
end
ifK<0%变压器支路k<0时的潮流
S（I,J）=bus（I,2）^2*（conj（Ym+Yt））+bus（I,2）*（cos（bus（I,3））+i*sin（bus（I,3）））*bus（J,2）*（cos（bus（J,3））-i*sin（bus（J,3）））*conj（Yt*K）;
S（J,I）=bus（J,2）^2*（conj（Yt））*K^2+bus（J,2）*（cos（bus（J,3））+i*sin（bus（J,3）））*bus（I,2）*（cos（bus（I,3））-i*sin（bus（I,3）））*conj（Yt*K）;
delS（I,J）=S（I,J）+S（J,I）;
end
ifJ==5&Zt==0
Sp=[line（k,1）line（k,2）S（I,5）0S（I,5）];
else
Sp=[line（k,1）line（k,2）S（I,J）S（J,I）delS（I,J）];
end
Pf（k,:
）=Sp;
end
%输入的参数数据：
%datafortestcase
%各节点参数：
节点编号,注入有功,注入无功,（Sn=100MVA）电压幅值,电压相位,类型
%类型：
1=PQ节点,2=PV节点,3=平衡节点
%（bus#）（volt）（ang）（p）（q）（bustype）
bus=[
1,1.0,0.0,-0.478,0.039,1;
2,1.0,0.0,-0.076,-0.016,1;
3,1.0,0.0,0.0,0.0,1;
4,1.0,0.0,-0.295,-0.166,1;
5,1.0,0.0,-0.09,-0.058,1;
6,1.0,0.0,-0.035,-0.018,1;
7,1.0,0.0,-0.061,-0.016,1;
8,1.0,0.0,-0.135,-0.058,1;
9,1.0,0.0,-0.149,-0.05,1;
10,1.045,0.0,0.183,0.0,2;
11,1.010,0.0,-0.942,0.0,2;
12,1.70,0.0,-0.112,0.047,2;
13,1.90,0.0,0.0,0.174,2;
14,1.060,0.0,0.0,0.0,3];
%各支路参数：
起点编号,终点编号,电阻,电抗,电导,电纳
line=[
1,2,0.01335,0.04211,0.0,0.0,0;
1,3,0.0,0.20912,0.0,0.0,0;
1,4,0.0,0.55618,0.0,0.0,0;
1,10,0.05811,0.17632,0.0,0.0340,0;
1,11,0.06701,0.17103,0.0,0.0128,0;
2,10,0.05695,0.17388,0.0,0.0346,0;
2,12,0.0,0.25202,0.0,0.0,0;
2,14,0.05403,0.22304,0.0,0.0492,0;
3,4,0.0,0.11001,0.0,0.0,0;
3,13,0.0,0.17615,0.0,0.0,0;
4,5,0.03181,0.08450,0.0,0.0,0;
4,9,0.12711,0.27038,0.0,0.0,0;
5,6,0.08205,0.19207,0.0,0.0,0;
6,12,0.09498,0.19890,0.0,0.0,0;
7,8,0.22092,0.19988,0.0,0.0,0;
7,12,0.12291,0.25581,0.0,0.0,0;
8,9,0.17093,0.34802,0.0,0.0,0;
8,12,0.06615,0.13027,0.0,0.0,0;
10,11,0.04699,0.19797,0.0,0.0438,0;
10,14,0.01938,0.05917,0.0,0.0528,0];
输出结果数据1：
--第1次迭代时雅可比矩阵--
-38.62403321.5785544.7819431.797979-0.000000-0.000000-0.000000-0.000000-0.000000
6.840981-0.000000-0.000000-0.0000005.3460515.119505-0.000000-0.000000-10.417258
-0.000000-0.000000-0.000000-0.000000
-0.000000-0.000000-0.000000-0.000000-0.00000021.578554-38.240787-0.000000-0.000000
6.7454965.427654-0.000000-0.000000-0.000000-0.000000-0.0000006.840981-9.429913
-0.000000-0.000000-0.000000-0.000000
-0.000000-0.0000009.090083-24.658288-0.0000004.781943-0.000000-0.000000-0.000000
-0
升级会员 