基于MATLAB的电力系统潮流计算.docx
《基于MATLAB的电力系统潮流计算.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的电力系统潮流计算.docx（19页珍藏版）》请在冰豆网上搜索。
基于MATLAB的电力系统潮流计算
基于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; %置迭代次数为初始值
%创建节点导纳矩阵
fori=1:
n1
ifB1（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
Y
OrgS=zeros（2*n-2,1）;
DetaS=zeros（2*n-2,1）; %将OrgS、DetaS初始化
%创建OrgS，用于存储初始功率参数
h=0;
j=0;
fori=1:
n %对PQ节点的处理
ifi~=isb&B2（i,6）==2
h=h+1;
forj=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）））+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
fori=1:
n %对PV节点的处理，注意这时不可再将h初始化为0
ifi~=isb&B2（i,6）==3
h=h+1;
forj=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）））+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
OrgS
%创建PVU用于存储PV节点的初始电压
PVU=zeros（n-h-1,1）;
t=0;
fori=1:
n
ifB2（i,6）==3
t=t+1;
PVU（t,1）=B2（i,3）;
end
end
PVU
%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量
h=0;
fori=1:
n %对PQ节点的处理
ifi~=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;
fori=1:
n %对PV节点的处理，注意这时不可再将h初始化为0
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n %对PQ节点的处理
ifB2（i,6）==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==（n-1）%将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=0;
end
end
end
end
end
k=0;
fori=1:
n %对PV节点的处理
ifB2（i,6）==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==（n-1） %将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=0;
end
end
end
end
end
Jacbi
%求解修正方程，获取节点电压的不平衡量
DetaU=zeros（2*n-2,1）;
DetaU=inv（Jacbi）*DetaS;
DetaU
%修正节点电压
j=0;
fori=1:
n %对PQ节点处理
ifB2（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
fori=1:
n %对PV节点的处理
ifB2（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
%开始循环**********************************************************************
whileabs（max（DetaU））>pr
OrgS=zeros（2*n-2,1）; %!
!
!
初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵
h=0;
j=0;
fori=1:
n
ifi~=isb&B2（i,6）==2
h=h+1;
forj=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）））+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
fori=1:
n
ifi~=isb&B2（i,6）==3
h=h+1;
forj=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）））+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
OrgS
%创建DetaS
h=0;
fori=1:
n
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n
ifi~=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;
fori=1:
n
ifB2（i,6）==2
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==（n-1）
k=0;
end
end
end
end
end
k=0;
fori=1:
n
ifB2（i,6）==3
h=h+1;
forj=1:
n
ifj~=isb
k=k+1;
ifi==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
ifk==（n-1）
k=0;
end
end
end
end
end
Jacbi
DetaU=zeros（2*n-2,1）;
DetaU=inv（Jacbi）*DetaS;
DetaU
%修正节点电压
j=0;
fori=1:
n
ifB2（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
fori=1:
n
ifB2（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
Times=Times+1; %迭代次数加1
end
Times
一个原始数据的例子
节点数5
支路数5
平衡节点编号5
精度pr0.000001
B1（支路参数矩阵）
[120.04+0.25i0.5i10;130.1+0.35i010;230.08+0.30i0.5i10;420.015i01.051;530.03i01.051]
B2（节点参数矩阵）
[0-1.6-0.8i1002;0-2-1i1002;0-3.7-1.3i1002;05+0i1.051.0503;001.051.0501]
X（节点号和对地参数）
[10;20;30;40;50]
电力系统潮流计算
——9结点算例-PQ法
原始数据录入data.txt文档:
标号,起始结点,终止结点,支路电阻参数,支路电抗参数,支路对地导纳参数
1,2,5,0.0,0.063,0.0,
2,5,9,0.019,0.072,0.075,
3,6,9,0.012,0.101,0.105,
4,3,6,0.0,0.059,0.0,
5,6,8,0.039,0.17,0.179,
6,4,8,0.017,0.092,0.079,
7,5,7,0.032,0.161,0.153,
8,4,7,0.01,0.085,0.08
升级会员 