Q分解法潮流计算matlab小程序.docx
- 文档编号:18413982
- 上传时间:2023-08-16
- 格式:DOCX
- 页数:21
- 大小:21.53KB
Q分解法潮流计算matlab小程序.docx
《Q分解法潮流计算matlab小程序.docx》由会员分享,可在线阅读,更多相关《Q分解法潮流计算matlab小程序.docx(21页珍藏版)》请在冰点文库上搜索。
Q分解法潮流计算matlab小程序
-Q分解法潮流计算matlab小程序1
(2010-06-2314:
43:
08)
转载▼
标签:
杂谈
题目是东南大学陈怡编的电力系统分析里,五节点的算例,本程序的通用性有待提高。
各元件的有关数据如下:
(1)发电机G1的端电压为1p.u,发出的有功、无功可调节;
(2)发电机G2的端电压为1p.u,按指定的有功P=0.5发电,取ε=0.0001;
(3)其中,G1为平衡节点,编号为5;G2为PV节点,编号为4;其余为PQ节点,编号为1、2、3。
其它数据见上图所示。
(4)本程序均采用标幺值进行计算,角度为度( ),但程序运行中均是用弧度进行计算。
在编程开始时,先创建了MATLAB数据,如下所示,仅在MATLAB命令窗口中有效:
systemdata.dat表示“系统数据”:
5 5 2 3 1 1
branch.dat表示“支路数据”:
表4-1支路数据
支路状态
节点号
节点号
R
X
B(K)
1
1
2
0.025
0.08
0.14
1
1
3
0.03
0.1
0.18
1
2
3
0.02
0.06
0.10
3
2
4
0
0.1905
1.0522
3
3
5
0
0.1905
1.0522
generator.dat表示“发电机数据”:
表4-2发电机数据
1
4
0.5
0
0.0150
0.0150
0.0150
0.0150
140.82
1
5
1.0
0
0.0150
0.0150
0.0150
0.0150
140.82
loads.dat表示“负荷数据”:
表4-3负荷数据
节点状态
节点号
P
Q
1
1
-0.8055
-0.532
1
2
-0.1800
-0.120
1
3
0
0
pvnode.dat表示“节点数据”:
4 1.0 0.5 4.60000 0.0000 1
balance.dat表示“平衡节点数据”:
5 1.0000
下面对本程序作简单的说明:
其中,node_counts表示“节点数”,branch_counts表示“支路数(含并联支路、接地支路)”,generator_counts表示“发电机数”,load_counts表示“负荷数”,pvnode_counts表示“PV节点数”,bal_counts表示“平衡节点数”,branch表示“支路数据”,generator表示“发电机数据”,loads表示“负荷数据”,pvnode表示“节点数据”,balance表示“平衡节点数据”,Y表示“节点导纳矩阵”,B1表示“B′”,B2表示“B″”,U表示“节点电压”,dPQ表示“有功无功的不平衡量”,K表示“迭代次数”,eps表示“收敛精度”,SIJ、SJI表示“支路潮流”,p表示“注入有功”,q表示“注入无功”。
主程序:
%===================算例为5节点5支路网(东南大学陈怡编)
%*********************************************************************
%* P-Q法潮流计算主程序 *
%*********************************************************************
%M函数在Workspace看不到变量的值
function[UPQ,SIJ,SJI,BalanceNodePQ,LossIJPQ,LossP,LossQ]=flowmain(systemdata,branch,generator,loads,pvnode,balance)
clear
clc
%定义全局变量2010-1
globalnode_counts;%节点数
globalbranch_counts;%支路数(含并联支路、接地支路)
globalgenerator_counts;%发电机数
globalload_counts;%负荷数
globalpvnode_counts;%PV节点数
globalbal_counts;%平衡节点数
globalbranch;%支路数据
globalgenerator;%发电机数据
globalloads;%负荷数据
globalpvnode;%PV节点数据
globalbalance;%平衡节点数据
globalY;%节点导纳矩阵
globalB1;%
globalB2;
globalU;%节点电压
globaldPQ;
globalKKK;%迭代次数
globaleps;%收敛精度
globalSIJ;%支路潮流
globalp;%注入有功
globalq;%注入无功
%下述方法仅在MATLAB命令窗口中有效
loadsystemdata.dat;%系统数据
loadbranch.dat;%支路数据
loadgenerator.dat;%发电机数据
loadloads.dat;%负荷数据
loadpvnode.dat;%PV节点数据
loadbalance.dat;%平衡节点数据
%赋初值
node_counts=systemdata
(1);
branch_counts=systemdata
(2);
generator_counts=systemdata(3);
load_counts=systemdata(4);
pvnode_counts=systemdata(5);
bal_counts=systemdata(6);
%%%节点导纳矩阵Y的形成
Y=produce_Y(node_counts,branch_counts,branch)
%%%形成B撇和B撇撇
B1=produce_B1(node_counts,balance,Y)
B2=produce_B2(node_counts,pvnode,B1,pvnode_counts)
%%给节点电压赋初值
[U]=first_SetVol(node_counts,bal_counts,pvnode_counts,balance,pvnode);
forii=1:
node_counts
disp(sprintf('初始电压(%d) %d%d',ii,U(ii,1),U(ii,2)));
end
%%%形成delta_p,delta_q
%%%形成循环迭代
esp=0.0001;
end_flag=0;
KKK=0;
tic
whileend_flag==0
[end_flag,U]=flow_solution(esp);
end
t=toc
%求平衡节点功率及各支路功率
[SIJ]=branch_flow(U,Y,pi);%显示潮流结果
end
子程序1(形成节点导纳矩阵):
%节点导纳矩阵的形成
function [jdY]=produce_Y(node_counts,branch_counts,branch)
forx=1:
node_counts
fory=1:
node_counts
jdY(x,y)=0;%各节点导纳Y清零
end
end
fort=1:
branch_counts
%接地支路
if branch(t,1)==4
zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗Z
ifzlZ~=0%三绕组变压器联结点
jdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ;
end
end
%变压器支路
if branch(t,1)==2|branch(t,1)==3
zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗Z
jdY(branch(t,2),branch(t,3))=-1/zlZ/branch(t,6);%算节点I与J节点间互导纳Yij
jdY(branch(t,3),branch(t,2))=-1/zlZ/branch(t,6);%算节点J与I节点间互导纳Yji
jdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ/branch(t,6)^2;%算I节点自导纳Yii
jdY(branch(t,3),branch(t,3))=jdY(branch(t,3),branch(t,3))+1/zlZ;%算J节点自导纳Yjj
end
%线路支路
if branch(t,1)==1
zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗Z
jdY(branch(t,2),branch(t,3))=-1/zlZ;%算节点I与J节点间互导纳Yij
jdY(branch(t,3),branch(t,2))=-1/zlZ;%算节点J与I节点间互导纳Yji
jdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ+i*branch(t,6)/2;%算I节点自导纳Yii
jdY(branch(t,3),branch(t,3))=jdY(branch(t,3),branch(t,3))+1/zlZ+i*branch(t,6)/2;;%算J节点自导纳Yjj
end
end
disp(sprintf('---------------节点导纳矩阵---------------'));%显示
子程序2(形成):
function[jdB1]=produce_B1(node_counts,balance,Y)
disp(sprintf('---------------B′---------------'));
fora=1:
node_counts
ifa==balance
(1)
Y(a,:
)=[];
y=Y;
y(:
a)=[];
y1=y;
end
end
jdB1=imag(y1);
子程序3(形成):
function[jdB2]=produce_B2(node_counts,pvnode,B1,pvnode_counts)
disp(sprintf('---------------B″---------------'));
fora=1:
node_counts
forb=1:
pvnode_counts
ifa==pvnode(b,1)
B1(a,:
)=[];
m=B1;
m(:
a)=[];
n=m;
end
end
end
jdB2=n;
子程序4(电压赋初值):
%给节点电压赋初值
function[U0]=first_SetVol(node_counts,bal_counts,pvnode_counts,balance,pvnode)
disp(sprintf('---------------节点电压赋初值---------------'));
fori=1:
node_counts
U0(i,1)=1.000;%PQ节点 电压值 相角
U0(i,2)=0;
end
forj=1:
bal_counts%重置平衡节点电压
U0(balance(j,1),1)=balance(j,2);
U0(balance(j,1),2)=0;
end
forj=1:
pvnode_counts%重置PV节点电压
U0(pvnode(j,1),1)=pvnode(j,2);
U0(pvnode(j,1),2)=0;
end
子程序5(潮流计算的迭代):
function[end_flag,U]=flow_solution(esp)
globalnode_counts;
globalbranch_counts;
globalgenerator_counts;
globalload_counts;
globalpvnode_counts;
globalbal_counts;
globalgenerator;
globalloads;
globalpvnode;
globalbalance;
globalp;%注入有功
globalq;%注入无功
globalesp;
globalB1;
globalB2;
globalY;
globalU;
globaldelta_p_q;
globalKKK;%迭代次数
globalpi;
%初始化
forii=1:
node_counts
p(ii)=0;
q(ii)=0;
end
pi=3.1415126/180;
forii=1:
node_counts
forjj=1:
node_counts
GIJ=real(Y(ii,jj));
BIJ=imag(Y(ii,jj));
p(ii)=p(ii)+U(ii,1)*U(jj,1)*(GIJ*cos((U(ii,2)-U(jj,2))*pi)+BIJ*sin((U(ii,2)-U(jj,2))*pi));%注入有功求和
q(ii)=q(ii)+U(ii,1)*U(jj,1)*(GIJ*sin((U(ii,2)-U(jj,2))*pi)-BIJ*cos((U(ii,2)-U(jj,2))*pi));%注入无功求和
end
end
NN_counts=1;
forjj=1:
load_counts%是发电机PQ节点
ifloads(jj,1)==1%是节PQ点
dPQ(NN_counts,1)=loads(jj,2);
dPQ(NN_counts,2)=loads(jj,3)-p(dPQ(NN_counts,1));
NN_counts=NN_counts+1;
end
end
forjj=1:
generator_counts%是发电机PQ节点
ifgenerator(jj,1)==1&generator(jj,10)==2%是发电机PQ节点
dPQ(NN_counts,1)=generator(jj,2);
dPQ(NN_counts,2)=generator(jj,3)-p(dPQ(NN_counts,1));
NN_counts=NN_counts+1;
end
end
MM_counts=1;
forjj=1:
load_counts%是发电机PQ节点
ifloads(jj,1)==1%是节PQ点
dPQ(NN_counts,1)=loads(jj,2);%存PQ节点号
dPQ(NN_counts,2)=loads(jj,4)-q(dPQ(NN_counts,1));%存无功失配功率dq
NN_counts=NN_counts+1;MM_counts=MM_counts+1;
end
end
NN_counts=NN_counts-1;
MM_counts=MM_counts-1;
dPQ=(dPQ)
%%%求解delta_p_div_v
delta_p_div_v=dPQ([1234],2)./U([1234],1)
delta_q_div_v=dPQ([567],2)./U([123],1)
%%%判断是否收敛
if(abs(max(delta_p_div_v)) end_flag=1; disp(sprintf('计算收敛.')); return; elseend_flag=0; end %%%解修正方程 v_multip_delta_theta=-inv(B1)*delta_p_div_v; forii=1: (node_counts-bal_counts) delta_theta(ii,1)=(v_multip_delta_theta(ii,1)/U(ii,1))/pi; end delta_v=-inv(B2)*delta_q_div_v; KKK=KKK+1; disp(sprintf('---------------迭代次数---K=%d----------------',KKK)); U([1234],2)=U([1234],2)+delta_theta; U([123],1)=U([123],1)+delta_v; forii=1: node_counts disp(sprintf('迭代过程电压(%d) %d %d',ii,U(ii,1),U(ii,2))); end return; 子程序6(功率分布的计算): %求解潮流 function[SIJ]=branch_flow(U,Y,pi,node_counts,branch,branch_counts) globalnode_counts; globalU; globalY; globalpi; globalbranch; globalbranch_counts; form=1: node_counts Uz(m)=U(m,1)*(cos((U(m,2))*pi)+i*sin((U(m,2))*pi)); end forp=1: node_counts C(p)=0; forq=1: node_counts C(p)=C(p)+conj(Y(p,q))*conj(Uz(q)); end S(p,1)=p; S(p,2)=Uz(p)*C(p); end disp('---------------各节点的功率S为----------------'); disp(S); disp('---------------各支路的功率SIJ为----------------'); form=1: branch_counts ifbranch(m,1)==1 SIJ(branch(m,2),branch(m,3))=Uz(branch(m,2))*(conj(Uz(branch(m,2)))*conj(i*branch(m,6)/2)+(conj(Uz(branch(m,2)))-conj(Uz(branch(m,3))))*conj(Y(branch(m,2),branch(m,3)))); else SIJ(branch(m,2),branch(m,3))=Uz(branch(m,2))*(conj(Uz(branch(m,2)))*conj((1-branch(m,6))/(branch(m,6)^2)*Y(branch(m,2),branch(m,3)))+(conj(Uz(branch(m,2)))-conj(Uz(branch(m,3))))*conj(Y(branch(m,2),branch(m,3)))); end end SIJ disp('---------------各支路的功率SJI为---------------'); form=1: branch_counts ifbranch(m,1)==1 SJI(branch(m,3),branch(m,2))=Uz(branch(m,3))*(conj(Uz(branch(m,3)))*conj(i*branch(m,6)/2)+(conj(Uz(branch(m,3)))-conj(Uz(branch(m,2))))*conj(Y(branch(m,3),branch(m,2)))); else SJI(branch(m,3),branch(m,2))=Uz(branch(m,3))*(conj(Uz(branch(m,3)))*conj((branch(m,6)-1)/(branch(m,6))*Y(branch(m,3),branch(m,2)))+(conj(Uz(branch(m,3)))-conj(Uz(branch(m,2))))*conj(Y(branch(m,3),branch(m,2)))); end end SJI disp('---------------各支路功率损耗为--------------
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 解法 潮流 计算 matlab 程序