用matlab电力系统潮流计算.docx
《用matlab电力系统潮流计算.docx》由会员分享，可在线阅读，更多相关《用matlab电力系统潮流计算.docx（47页珍藏版）》请在冰点文库上搜索。
用matlab电力系统潮流计算
题目：
潮流计算与matlab
教学单位电气信息学院
姓名
学号
年级
专业电气工程及其自动化
指导教师
职称副教授
摘要
电力系统稳态分析包括潮流计算和静态安全分析。
本文主要运用的事潮流计算，潮流计算是电力网络设计与运行中最基本的运算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。
本位就是运用潮流计算具体分析，并有MATLAB仿真。
关键词：
电力系统潮流计算MATLAB
Abstract
Electricpowersystemsteadyflowcalculationandanalysisofthestaticsafetyanalysis.Thispaper,bymeansofthecalculation,flowcalculationisthetrendofthepowernetworkdesignandoperationofthemostbasicoperationsofelectricpowernetwork,variousdesignschemeandtheoperationwaystotidecomputation,cangetallkindsofeachnodeofthepowergridvoltageandseekthetrendofthenetworkandthenetworkofthecomponentsofthepowerloss,andgettingelectricpower.Thestandardistousethepowerflowcalculationandanalysis,thespecifichaveMATLABsimulation.
Keywords:
Powersystem;Flowcalculation;MATLABsimulation
1任务提出与方案论证
潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法,是电力系统规划和运营中不可缺少的一个重要组成部分。
可以说,它是电力系统分析中最基本、最重要的计算,是系统安全、经济分析和实时控制与调度的基础。
常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。
潮流计算的结果是电力系统稳定计算和故障分析的基础。
在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。
同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。
因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。
在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。
是电力系统研究人员长期研究的一个课题。
它既是对电力系统规划设计和运行方式的合理性、可靠性及经济性进行定量分析的依据，又是电力系统静态和暂态稳定计算的基础。
潮流计算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都以计算机的应用为前提用计算机进行潮流计算主要步骤在于编制计算机程序，这是一项非常复杂的工作。
对系统进行潮流分析，本文利用MATLAB中的SimpowerSystems工具箱设计电力系统，在simulink环境下，不仅可以仿真系统的动态过程，还可以对系统进行稳态潮流分析。
2总体设计
SimpowerSystems使用Simulink环境，可以将该系统中的发电机、变压器，线路等模型联结起来，形成电力系统仿真模拟图。
在加人测量模块，并对各元件的参数进行设置后，用measurement和sink中的仪器可以观察各元件的电压、电流、功率的大小。
2.1潮流计算等值电路
2.2建立电力系统模型
在Simulink中按照电力系统原型选择元件进行建模。
所建立的模型和建立的方法在详细设计中详述。
在电力系统模型的建立工程中主要涉及到的是：
元器件的选择及其参数的设置；发电机选型；变压器选择；线路的选择；负荷模型的选择；母线选择。
2.3模型的调试与运行
建立系统模型，并设置好参数以后，就可以在Simulink环境下进行仿真运行。
运行的具体结果和分析也在详细设计中详述。
3详细设计
3.1计算前提
首先是发电机的参数计算，先对5个发电厂简化为5台发电机来计算。
发电机G1：
发电机G2：
发电机G3：
发电机G4：
发电机G5：
其次是变电站的参数计算，我们还是对7个变电站简化为7台变压器来计算。
变压器T1：
变压器T2：
（双并联）
变压器T3：
（四并联）
变压器T4：
（双并联）
变压器T5：
变压器T6：
（两个三绕组变压器并联）
变压器T7：
（双并联）
再次是传输线参数计算，5条传输线的具体计算如下。
根据教材查得
线路L1：
线路L2：
线路L3：
（双回路）
线路L4：
线路L5：
（双回路）
3.2手工计算
FLR1：
FLR2：
FLR3：
FLR4：
FLR5:
计算每一个FLR的功率分布和电压分布计算如下：
FLR1：
FLR2：
功率分布：
电压分布：
FLR3：
功率分布：
电压分布：
FLR4：
功率分布：
电压分布：
FLR5:
这里我们先将f点和发电机G5当做电源，经过
和
构成两端供电网络以g点作为运算负荷进行计算。
电压分布：
4设计图及源程序
4.1MATLAB仿真
相关的原始数据输入格式如下：
1、B1是支路参数矩阵，第一列和第二列是节点编号。
节点编号由小到大编写。
2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器。
3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
4、X为节点号和对地参数矩阵。
其中第一列为节点编号，第二列为节点对地参数。
4.2潮流计算源程序
%本程序的功能是用牛顿——拉夫逊法进行11节点潮流计算
clear;
n=11;%input（'请输入节点数：
n='）;
nl=11;%input（'请输入支路数：
nl='）;
isb=1;%input（'请输入平衡母线节点号：
isb='）;
pr=0.00001;%input（'请输入误差精度：
pr='）;
B1=[120.03512+0.08306i0.13455i10;
230.0068+0.18375i01.023811;
140.05620+0.13289i0.05382i10;
450.00811+0.24549i01.023811;
160.05620+0.13289i0.05382i10;
460.04215+0.09967i0.04037i10;
670.0068+0.18375i01.023811;
680.02810+0.06645i0.10764i10;
140.05620+0.13289i0.05382i10;
8100.00811+0.24549i011;
890.03512+0.08306i0.13455i10]
B2=[001.11.101;
001002;
00.343+0.21256i1002;
001002;
00.204+0.12638i1002;
001002;
00.306+0.18962i1002;
001002;
0.501.11.103;
00.343+0.21256i1002;
001002]
;
%B1矩阵：
1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳
%5、支路的变比；6、支路首端处于K侧为1，1侧为0
%B2矩阵：
1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值
%4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量
%6、节点分类标号：
1为平衡节点（应为1号节点）；2为PQ节点；
%3为PV节点；
%input（'请输入各节点参数形成的矩阵：
B2='）;
Y=zeros（n）;e=zeros（1,n）;f=zeros（1,n）;V=zeros（1,n）;sida=zeros（1,n）;S1=zeros（nl）;
%%%---------------------------------------------------
fori=1:
nl%支路数
ifB1（i,6）==0%左节点处于1侧
p=B1（i,1）;q=B1（i,2）;
else%左节点处于K侧
p=B1（i,2）;q=B1（i,1）;
end
Y（p,q）=Y（p,q）-1./（B1（i,3）*B1（i,5））;%非对角元
Y（q,p）=Y（p,q）;%非对角元
Y（q,q）=Y（q,q）+1./（B1（i,3）*B1（i,5）^2）+B1（i,4）./2;%对角元K侧
Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）./2;%对角元1侧
end
%求导纳矩阵
disp（'导纳矩阵Y='）;
disp（Y）
%----------------------------------------------------------
G=real（Y）;B=imag（Y）;%分解出导纳阵的实部和虚部
fori=1:
n%给定各节点初始电压的实部和虚部
e（i）=real（B2（i,3））;
f（i）=imag（B2（i,3））;
V（i）=B2（i,4）;%PV节点电压给定模值
end
fori=1:
n%给定各节点注入功率
S（i）=B2（i,1）-B2（i,2）;%i节点注入功率SG-SL
B（i,i）=B（i,i）+B2（i,5）;%i节点无功补偿量
end
%===================================================================
P=real（S）;Q=imag（S）;%分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2
whileIT2~=0%N0=2*n雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
fori=1:
n
ifi~=isb%非平衡节点
C（i）=0;D（i）=0;
forj1=1:
n
C（i）=C（i）+G（i,j1）*e（j1）-B（i,j1）*f（j1）;%Σ（Gij*ej-Bij*fj）
D（i）=D（i）+G（i,j1）*f（j1）+B（i,j1）*e（j1）;%Σ（Gij*fj+Bij*ej）
end
P1=C（i）*e（i）+f（i）*D（i）;%节点功率P计算eiΣ（Gij*ej-Bij*fj）+fiΣ（Gij*fj+Bij*ej）
Q1=C（i）*f（i）-e（i）*D（i）;%节点功率Q计算fiΣ（Gij*ej-Bij*fj）-eiΣ（Gij*fj+Bij*ej）
%求i节点有功和无功功率P',Q'的计算值
V2=e（i）^2+f（i）^2;%电压模平方
%=========以下针对非PV节点来求取功率差及Jacobi矩阵元素=========
ifB2（i,6）~=3%非PV节点
DP=P（i）-P1;%节点有功功率差
DQ=Q（i）-Q1;%节点无功功率差
%===============以上为除平衡节点外其它节点的功率计算=================
%=================求取Jacobi矩阵===================
forj1=1:
n
ifj1~=isb&j1~=i%非平衡节点&非对角元
X1=-G（i,j1）*e（i）-B（i,j1）*f（i）;%dP/de=-dQ/df
X2=B（i,j1）*e（i）-G（i,j1）*f（i）;%dP/df=dQ/de
X3=X2;%X2=dp/dfX3=dQ/de
X4=-X1;%X1=dP/deX4=dQ/df
p=2*i-1;q=2*j1-1;
J（p,q）=X3;J（p,N）=DQ;m=p+1;%X3=dQ/deJ（p,N）=DQ节点无功功率差
J（m,q）=X1;J（m,N）=DP;q=q+1;%X1=dP/deJ（m,N）=DP节点有功功率差
J（p,q）=X4;J（m,q）=X2;%X4=dQ/dfX2=dp/df
elseifj1==i&j1~=isb%非平衡节点&对角元
X1=-C（i）-G（i,i）*e（i）-B（i,i）*f（i）;%dP/de
X2=-D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dP/df
X3=D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dQ/de
X4=-C（i）+G（i,i）*e（i）+B（i,i）*f（i）;%dQ/df
p=2*i-1;q=2*j1-1;J（p,q）=X3;J（p,N）=DQ;%扩展列△Q
m=p+1;
J（m,q）=X1;q=q+1;J（p,q）=X4;J（m,N）=DP;%扩展列△P
J（m,q）=X2;
end
end
else
%===============下面是针对PV节点来求取Jacobi矩阵的元素===========
DP=P（i）-P1;%PV节点有功误差
DV=V（i）^2-V2;%PV节点电压误差
forj1=1:
n
ifj1~=isb&j1~=i%非平衡节点&非对角元
X1=-G（i,j1）*e（i）-B（i,j1）*f（i）;%dP/de
X2=B（i,j1）*e（i）-G（i,j1）*f（i）;%dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J（p,q）=X5;J（p,N）=DV;%PV节点电压误差
m=p+1;
J（m,q）=X1;J（m,N）=DP;q=q+1;J（p,q）=X6;%PV节点有功误差
J（m,q）=X2;
elseifj1==i&j1~=isb%非平衡节点&对角元
X1=-C（i）-G（i,i）*e（i）-B（i,i）*f（i）;%dP/de
X2=-D（i）+B（i,i）*e（i）-G（i,i）*f（i）;%dP/df
X5=-2*e（i）;
X6=-2*f（i）;
p=2*i-1;q=2*j1-1;J（p,q）=X5;J（p,N）=DV;%PV节点电压误差
m=p+1;
J（m,q）=X1;J（m,N）=DP;q=q+1;J（p,q）=X6;%PV节点有功误差
J（m,q）=X2;
end
end
end
end
end
%=========以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差=====================
fork=3:
N0%N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;%N=N0+1即N=2*n+1扩展列△P、△Q或△U
fork2=k1:
N1%从k+1列的Jacobi元素到扩展列的△P、△Q或△U
J（k,k2）=J（k,k2）./（J（k,k）+eps）;%用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J（k,k）=1;%对角元规格化K行K列对角元素赋1
%====================回代运算=======================================
ifk~=3%不是第三行k>3
k4=k-1;
fork3=3:
k4%用k3行从第三行开始到当前行的前一行k4行消去
fork2=k1:
N1%k3行后各行上三角元素
J（k3,k2）=J（k3,k2）-J（k3,k）*J（k,k2）;%消去运算（当前行k列元素消为0）
end%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J（k3,k）=0;%当前行第k列元素已消为0
end
ifk==N0%若已到最后一行
break;
end
%==================前代运算==================================
fork3=k1:
N0%从k+1行到2*n最后一行
fork2=k1:
N1%从k+1列到扩展列消去k+1行后各行下三角元素
J（k3,k2）=J（k3,k2）-J（k3,k）*J（k,k2）;%消去运算
end%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J（k3,k）=0;%当前行第k列元素已消为0
end
else%是第三行k=3
%======================第三行k=3的前代运算========================
fork3=k1:
N0%从第四行到2n行（最后一行）
fork2=k1:
N1%从第四列到2n+1列（即扩展列）
J（k3,k2）=J（k3,k2）-J（k3,k）*J（k,k2）;%消去运算（当前行3列元素消为0）
end%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J（k3,k）=0;%当前行第3列元素已消为0
end
end
end
%====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
fork=3:
2:
N0-1
L=（k+1）./2;
e（L）=e（L）-J（k,N）;%修改节点电压实部
k1=k+1;
f（L）=f（L）-J（k1,N）;%修改节点电压虚部
end
%------修改节点电压-----------
fork=3:
N0
DET=abs（J（k,N））;
ifDET>=pr%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2（a）=IT2;%不满足要求的节点数
ICT1=ICT1+1;%迭代次数
end
%用高斯消去法解"w=-J*V"
disp（'迭代次数：
'）;
disp（ICT1）;
disp（'没有达到精度要求的个数：
'）;
disp（ICT2）;
fork=1:
n
V（k）=sqrt（e（k）^2+f（k）^2）;%计算各节点电压的模值
sida（k）=atan（f（k）./e（k））*180./pi;%计算各节点电压的角度
E（k）=e（k）+f（k）*j;%将各节点电压用复数表示
end
%===============计算各输出量===========================
disp（'各节点的实际电压标幺值E为（节点号从小到大排列）：
'）;
disp（E）;%显示各节点的实际电压标幺值E用复数表示
disp（'-----------------------------------------------------'）;
disp（'各节点的电压大小V为（节点号从小到大排列）：
'）;
disp（V）;%显示各节点的电压大小V的模值
disp（'-----------------------------------------------------'）;
disp（'各节点的电压相角sida为（节点号从小到大排列）：
'）;
disp（sida）;%显示各节点的电压相角
forp=1:
n
C（p）=0;
forq=1:
n
C（p）=C（p）+conj（Y（p,q））*conj（E（q））;%计算各节点的注入电流的共轭值
end
S（p）=E（p）*C（p）;%计算各节点的功率S=电压X注入电流的共轭值
end
disp（'各节点的功率S为（节点号从小到大排列）：
'）;
disp（S）;%显示各节点的注入功率
disp（'-----------------------------------------------------'）;
disp（'各条支路的首端功率Si为（顺序同您输入B1时一致）：
'）;
fori=1:
nl
p=B1（i,1）;q=B1（i,2）;
ifB1（i,6）==0
Si（p,q）=E（p）*（conj（E（p））*conj（B1（i,4）./2）+（conj（E（p）*B1（i,5））...
-conj（E（q）））*conj（1./（B1（i,3）*B1（i,5））））;
Siz（i）=Si（p,q）;
else
Si（p,q）=E（p）*（conj（E（p））*conj（B1（i,4）./2）+（conj（E（p）./B1（i,5））...
-conj（E（q）））*conj（1./（B1（i,3）*B1（i,5））））;
Siz（i）=Si（p,q）;
end
disp（Si（p,q））;
SSi（p,q）=Si（p,q）;
ZF=['S（',num2str（p）,',',num2str（q）,'）=',num2str（SSi（p,q））];
disp（ZF）;
disp（'-----------------------------------------------------'）;
end
disp（'各条支路的末端功率Sj为（顺序同您输入B1时一致）：
'）;
fori=1:
nl
p=B1（i,1）;q=B1（i,2）;
ifB1（i,6）=
升级会员 