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表示“系统数据”：552311branch.dat表示“支路数据”：表4-1支路数据支路状态节点号节点号RXB（K）1120.0250.080.141130.030.10.181230.020.060.1032400.19051.052233500.19051.0522generator.dat表示“发电机数据”：表4-2发电机数据140.500.01500.01500.01500.0150140.82151.000.01500.01500.01500.0150140.82loads.dat表示“负荷数据”：表4-3负荷数据节点状态节点号PQ11-0.8055-0.53212-0.1800-0.1201300pvnode.dat表示“节点数据”：41.00.54.600000.00001balance.dat表示“平衡节点数据”：51.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)clearclc%定义全局变量2010-1globalnode_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_countsdisp(sprintf('初始电压(%d)%d%d',ii,U(ii,1),U(ii,2)));end%%%形成delta_p,delta_q%%%形成循环迭代esp=0.0001;end_flag=0;KKK=0;ticwhileend_flag==0[end_flag,U]=flow_solution(esp);endt=toc%求平衡节点功率及各支路功率[SIJ]=branch_flow(U,Y,pi);%显示潮流结果end子程序1（形成节点导纳矩阵）：%节点导纳矩阵的形成function[jdY]=produce_Y(node_counts,branch_counts,branch)forx=1:node_countsfory=1:node_countsjdY(x,y)=0;%各节点导纳Y清零endendfort=1:branch_counts%接地支路ifbranch(t,1)==4zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗ZifzlZ~=0%三绕组变压器联结点jdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ;endend%变压器支路ifbranch(t,1)==2|branch(t,1)==3zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗ZjdY(branch(t,2),branch(t,3))=-1/zlZ/branch(t,6);%算节点I与J节点间互导纳YijjdY(branch(t,3),branch(t,2))=-1/zlZ/branch(t,6);%算节点J与I节点间互导纳YjijdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ/branch(t,6)^2;%算I节点自导纳YiijdY(branch(t,3),branch(t,3))=jdY(branch(t,3),branch(t,3))+1/zlZ;%算J节点自导纳Yjjend%线路支路ifbranch(t,1)==1zlZ=branch(t,4)+i*branch(t,5);%算各支路阻抗ZjdY(branch(t,2),branch(t,3))=-1/zlZ;%算节点I与J节点间互导纳YijjdY(branch(t,3),branch(t,2))=-1/zlZ;%算节点J与I节点间互导纳YjijdY(branch(t,2),branch(t,2))=jdY(branch(t,2),branch(t,2))+1/zlZ+i*branch(t,6)/2;%算I节点自导纳YiijdY(branch(t,3),branch(t,3))=jdY(branch(t,3),branch(t,3))+1/zlZ+i*branch(t,6)/2;;%算J节点自导纳Yjjendenddisp(sprintf('---------------节点导纳矩阵---------------'));%显示子程序2（形成）：function[jdB1]=produce_B1(node_counts,balance,Y)disp(sprintf('---------------B′---------------'));fora=1:node_countsifa==balance(1)Y(a,:)=[];y=Y;y(:,a)=[];y1=y;endendjdB1=imag(y1);子程序3（形成）：function[jdB2]=produce_B2(node_counts,pvnode,B1,pvnode_counts)disp(sprintf('---------------B″---------------'));fora=1:node_countsforb=1:pvnode_countsifa==pvnode(b,1)B1(a,:)=[];m=B1;m(:,a)=[];n=m;endendendjdB2=n;子程序4（电压赋初值）：%给节点电压赋初值function[U0]=first_SetVol(node_counts,bal_counts,pvnode_counts,balance,pvnode)disp(sprintf('---------------节点电压赋初值---------------'));fori=1:node_countsU0(i,1)=1.000;%PQ节点电压值相角U0(i,2)=0;endforj=1:bal_counts%重置平衡节点电压U0(balance(j,1),1)=balance(j,2);U0(balance(j,1),2)=0;endforj=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_countsp(ii)=0;q(ii)=0;endpi=3.1415126/180;forii=1:node_countsforjj=1:node_countsGIJ=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(j