基于matlab的PQ分解法程序

%clear;%參數(shù)設(shè)置區(qū)域：%線路參數(shù)%tr=【入端，出端，電阻，電抗，電導(dǎo)，電納,有無調(diào)壓變壓器，變壓器變比】；tr=[1,2,0.01938,0.05917,0,0.02640*2,0,1;2,3,0.04699,0.19797,0,0.02190*2,0,1;2,4,0.05811,0.17632,0,0.01870*2,0,1;1,5,0.05403,0.22304,0,0.02460*2,0,1;2,5,0.05695,0.17388,0,0.01700*2,0,1;3,4,0.06701,0.17103,0,0.01730*2,0,1;4,5,0.01335,0.04211,0,0.00640*2,0,1;7,8,0,0.17615,0,0,0,1;7,9,0,0.11001,0,0,0,1;9,10,0.03181,0.08450,0,0,0,1;6,11,0.09498,0.19890,0,0,0,1;6,12,0.12291,0.15581,0,0,0,1;6,13,0.06615,0.13027,0,0,0,1;9,14,0.12711,0.27038,0,0,0,1;10,11,0.08205,0.19207,0,0,0,1;12,13,0.22092,0.19988,0,0,0,1;13,14,0.17093,0.34802,0,0,0,1;4,7,0,0.20912,0,0,1,1.0225;4,9,0,0.55618,0,0,1,1.0320;5,6,0,0.25202,0,0,1,1.0730;];%1是平衡節(jié)點，%2是PV節(jié)點，%3是PQ節(jié)點；%節(jié)點編號類型發(fā)電機發(fā)出有功負荷吸收復(fù)功率電壓幅值電壓相位node=[1,1,0,0,1.0600,0;2,2,0.40,0.2170+0.1270i,1.045,0;3,2,0,0.9420+0.1900i,1.0100,0;4,3,0,0.4780-0.0390i,1,0;5,3,0,0.0760+0.0160i,1,0;6,2,0,0.1120+0.0750i,1.0700,0;7,3,0,0,1,0;8,2,0,0,1.0900,0;9,3,0,0.295+0.166i,1,0;10,3,0,0.0900+0.0580i,1,0;11,3,0,0.0350+0.0180i,1,0;12,3,0,0.0610+0.0160i,1,0;13,3,0,0.1350+0.0580i,1,0;14,3,0,0.1490+0.0500i,1,0;];%對地阻抗矩陣，節(jié)點編號、對地電阻、對地電納；GND=[9,0,0.19;];ngnd=1;%對地支路數(shù)nbb=14;%節(jié)點數(shù)ntr=20;%線路數(shù)nitmax=1000;%最大迭代次數(shù)error=1e-12;%允許的誤差%調(diào)整輸入矩陣元素順序pvjd=0;fori=1:nbb%找出PV節(jié)點個數(shù)if(node(i,2)==2)pvjd=pvjd+1;endendfori=1:nbb%找出平衡節(jié)點序號if(node(i,2)==1)nban=i;endendnode1=zeros(1,6);%將平衡節(jié)點移至最后一行，以便于架構(gòu)Jacobiannode1(1,:)=node(nban,:);node(nban,:)=node(nbb,:);node(nbb,:)=node1(1,:);fork=1:ntr%將線路矩陣中相應(yīng)序號對換i=1;fori=1:2if(tr(k,i)==nban)tr(k,i)=nbb;elseif(tr(k,i)==nbb)tr(k,i)=nban;endendend%tr%nodeYR=zeros(nbb,nbb);%初始化各導(dǎo)納矩陣YI=zeros(nbb,nbb);Y=zeros(nbb,nbb);%調(diào)壓變壓器i=1;fori=1:ntr%含有雙繞組調(diào)壓變壓器支路的導(dǎo)納矩陣if(tr(i,7)==1)mo=tr(i,3)^2+tr(i,4)^2;YR(tr(i,1),tr(i,1))=(tr(i,8)^2)*tr(i,3)/mo+YR(tr(i,1),tr(i,1));YI(tr(i,1),tr(i,1))=-(tr(i,8)^2)*tr(i,4)/mo+YI(tr(i,1),tr(i,1));YR(tr(i,2),tr(i,2))=tr(i,3)/mo+YR(tr(i,2),tr(i,2));YI(tr(i,2),tr(i,2))=-tr(i,4)/mo+YI(tr(i,2),tr(i,2));YR(tr(i,1),tr(i,2))=-tr(i,8)*tr(i,3)/mo+YR(tr(i,1),tr(i,2));YI(tr(i,1),tr(i,2))=-tr(i,8)*-tr(i,4)/mo+YI(tr(i,1),tr(i,2));YR(tr(i,2),tr(i,1))=-tr(i,8)*tr(i,3)/mo+YR(tr(i,2),tr(i,1));YI(tr(i,2),tr(i,1))=-tr(i,8)*-tr(i,4)/mo+YI(tr(i,2),tr(i,1));endendY=YR+YI*1i;forcc=1:ntr%線路的導(dǎo)納矩陣i=tr(cc,1);j=tr(cc,2);mo=tr(cc,3)^2+tr(cc,4)^2;if(tr(cc,7)~=1)YR(i,i)=YR(i,i)+tr(cc,3)/mo+0.5*tr(cc,5);YI(i,i)=YI(i,i)-tr(cc,4)/mo+0.5*tr(cc,6);YR(i,j)=YR(i,j)-tr(cc,3)/mo;YI(i,j)=YI(i,j)+tr(cc,4)/mo;YR(j,i)=YR(j,i)-tr(cc,3)/mo;YI(j,i)=YI(j,i)+tr(cc,4)/mo;YR(j,j)=YR(j,j)+tr(cc,3)/mo+0.5*tr(cc,5);YI(j,j)=YI(j,j)-tr(cc,4)/mo+0.5*tr(cc,6);endendY=YR+YI*1i;forcc=1:ngnd%對地并聯(lián)支路修正；i=GND(cc,1);%mo=GND(cc,2)^2+GND(cc,3)^2;YR(i,i)=YR(i,i)+GND(cc,2);YI(i,i)=YI(i,i)-GND(cc,3);endY=YR+YI*1i%輸出最終的導(dǎo)納矩陣；X=zeros(1,2*nbb-2-pvjd);DPQ=zeros(1,2*nbb-2-pvjd);f=0;nit=1;snet=zeros(1,nbb);%每個節(jié)點處發(fā)電機與負載的凈注入功率Snet：fori=1:nbbsnet(i)=node(i,3)-node(i,4);end%snetwhile(f==0&&nit<=nitmax)str=zeros(1,nbb);%每個節(jié)點各條線路的注入凈功率Str（或稱為計算功率）：%strfori=1:nbbpsum=0;qsum=0;forj=1:nbbpsum=psum+node(i,5)*node(j,5)*YR(i,j)*cos(node(i,6)-node(j,6))+node(i,5)*node(j,5)*YI(i,j)*sin(node(i,6)-node(j,6));%計算每兩個節(jié)點間傳輸有功qsum=qsum+node(i,5)*node(j,5)*YR(i,j)*sin(node(i,6)-node(j,6))-node(i,5)*node(j,5)*YI(i,j)*cos(node(i,6)-node(j,6));%計算每兩個節(jié)點間傳輸無功endstr(i)=psum+qsum*1i;%將各個節(jié)點的計算功率存入數(shù)組中。end%strP=zeros(1,nbb);Q=zeros(1,nbb);fori=1:nbbP(i)=real(str(i));%提取各個節(jié)點的有功功率并存入數(shù)組Q(i)=imag(str(i));%提取各個節(jié)點的無功功率并存入數(shù)組end%P%QDS=snet-str;%計算功率偏差值：DP=zeros(1,nbb);DQ=zeros(1,nbb);%DSfori=1:nbb%根據(jù)節(jié)點類型校正功率偏差if(node(i,2)==1)%平衡節(jié)點不參加迭代，功率偏差設(shè)為0DS(i)=0;elseif(node(i,2)==2)%若是PV節(jié)點,則只有有功偏差方程，無功偏差設(shè)為零y=imag(DS(i));DS(i)=DS(i)-y*1i;endend%DSfori=1:nbb%提取有功、無功偏差絕對值的最大值用于進行誤差比較DP(i)=abs(real(DS(i)))+DP(i);DQ(i)=abs(imag(DS(i)))+DQ(i);end%DP%DQA=max(DP);B=max(DQ);MAX=max(A,B);if(MAX<error)f=1;%若誤差滿足要求，則標志位置1，表示停止迭代endDP=zeros(1,nbb);DQ=zeros(1,nbb);fori=1:nbb%提取實際的無功、有功偏差。DP(i)=real(DS(i))+DP(i);DQ(i)=imag(DS(i))+DQ(i);endDVA=zeros(nbb,nbb);DVM=zeros(nbb,nbb);if(nit==1)B1=zeros(nbb,nbb);B2=zeros(nbb,nbb);B1=-YI;B2=-YI;forii=1:nbbif(node(ii,2)==1)forjj=1:nbbifii==jjB1(ii,ii)=1;B2(ii,ii)=1;elseB1(ii,jj)=0;B1(jj,ii)=0;B2(ii,jj)=0;B2(jj,ii)=0;endendendif(node(ii,2)==2)forjj=1:nbbifii==jjB2(ii,ii)=1;elseB2(ii,jj)=0;B2(jj,ii)=0;endendendendend%DP%DQDVA=B1\DP';DVM=B2\DQ';kk=1;forii=1:nbbD(1,kk)=DVA(ii);D(1,kk+1)=DVM(ii);kk=kk+2;endiii=1;forii=1:nbbnode(ii,6)=node(ii,6)+D(1,iii)/node(ii,5);node(ii,5)=node(ii,5)+D(1,iii+1)/node(ii,5);iii=iii+2;endnit=nit+1;enddisp('允許最大迭代次數(shù)')nitmaxdisp('本次計算迭代次數(shù)')nit=nit-1%計算功率潮流Strcal=zeros(ntr,5);V=zeros(nbb,1);%計算每個節(jié)點的電壓相量fori=1:nbbV(i)=node(i,5)*cos(node(i,6))+node(i,5)*sin(node(i,6))*1i;enddisp('每個節(jié)點的電壓相量')Vfork=1:ntr%計算每一條輸電線路首段輸入功率i=tr(k,1);j=tr(k,2);if(tr(k,7)==1)T=tr(k,8);elseT=1;endS=node(i,5)^2*(Y(i,j)'*T+tr(k,6)'/2*1i)-V(i)*V(j)'*Y(i,j)';Strcal(k,3)=-S;Strcal(k,1)=tr(k,1);Strcal(k,2)=tr(k,2);endfork=1:ntr%計算每一條輸電線路末端輸入功率j=tr(k,1);i=tr(k,2);if(tr(k,7)==1)T=tr(k,8);elseT=1;endS=node(i,5)^2*(Y(i,j)'/T+tr(k,6)'/2*1i)-V(i)*V(j)'*Y(i,j)';Strcal(k,4)=-S;endfork=1:ntr%計算每一條線路上的功率損耗Strcal(k,5)=Strcal(k,3)+Strcal(k,4);endfork=1:ntr%換回原來的序號i=1;fori=1:2if(Strcal(k,i)==nban)Strcal(k,i)=nbb;elseif(Strcal(k,i)==nbb)Strcal(k,i)=nban;endendenddisp('線路、支路的潮流計算結(jié)果')disp('入端序號出端序號入端注入功率出端注入功率功率損耗')Strcalsum=0;%計算平衡節(jié)點的發(fā)電機輸出功率fori=1:nbbif(node(i,2)==1)forj=1:nbbsum=Y(i,j)'*V(j)'+sum;endSnod