九节点系统潮流计算编程牛N_R法.doc

《九节点系统潮流计算编程牛N_R法.doc》由会员分享，可在线阅读，更多相关《九节点系统潮流计算编程牛N_R法.doc（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、如图所示系统，试计算潮流分布，相关数据见 PSASP7.0版潮流计算用户手册P121。#include #include float divRe(float b1,float b2,float b3,float b4) float a1r; a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4); return(a1r); float divIm(float b1,float b2,float b3,float b4) float a1i; a1i=(b2*b3-b1*b4)/(b3*b3+b4*b4); return(a1i); float mulRe(float b1,float

2、b2,float b3,float b4) float a2r; a2r=b1*b3-b2*b4; return(a2r); float mulIm(float b1,float b2,float b3,float b4) float a2i; a2i=b2*b3+b1*b4; return(a2i); float Max(float a,int n) int i;float max;max= fabs(a0);for(i=1;imax)max=fabs(ai); return(max); void main() int i,j,k,h,km; int T=16; float eps,sump

3、i1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,I1r,I1i,t,xx,xxx; float pi08,qi08,detpi8,detqi8,Iir08,Iii08,J01616, detsi16,detui16, delta_p99,delta_q99, a1632,ni1616,H88,N88, J88,L88,ei19,fi19,sp99,sq99; static float ybr99= 3.3074,-1.3652,0,0,0,-1.9422,0,0,0, -1.3652,2.5528,-1.1876,0,0,0,0,0,0, 0,-1.1876,2

4、.8047,-1.6171,0,0,0,0,0,0,0,-1.6171,2.7722,-1.1551,0,0,0,0, 0,0,0,-1.1551,2.437,-1.282,0,0,0, -1.9422,0,0,0,-1.282,3.2242,0,0,0, 0, 0, 0; static float ybi99= -39.3089,11.6041,0,0,0,10.5107,0,0,17.3611, 11.6041,-17.3382,5.9751,0,0,0,0,0,0, 0,5.9751,-35.4456,13.6980,0,0,16,0,0, 0,0,13.6980,-23.3033,9.

5、7843,0,0,0,0, 0,0,0,9.782,-32.1538,5.5882,0,17.0648,0, 10.5107,0,0,0,5.5882,-15.841,0,0,0, 0,0,16,0,0,0,-16,0,0, 0,0,0,0,17.0648,0,0,-17.0648,0, 17.3611,0,0,0,0,0,0,0,-17.3611;static float yd99=0,0.088,0,0,0,0.079,0,0,0,0.088,0,0.153,0,0,0,0,0,0,0,0.153,0,0.0745,0,0,0,0,0,0,0,0.0745,0,0.1045,0,0,0,0

6、,0,0,0,0.1045,0,0.179,0,0,0,0.079,0,0,0,0.179,0,0,0,0,0,0,0 ; float ei09=1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.04; float fi09=0.0; float pi9=0,-1.25,0,-1.0,0,-0.9,1.63,0.85,0; float qi9=0,-0.5,0,-0.35,0,-0.3,0,0,0; h=0; km=15; eps=0.00001; do h+=1; printf(nNow The %dth timesn,h); for(i=0;i8;i+) printf(e

7、i0%d=%ft,i,ei0i); printf(fi0%d=%ft,i,fi0i); for(i=0;i8;i+) printf(pi%d=%ft,i,pii); printf(qi%d=%ft,i,qii); sumpi2=0; sumqi2=0; for(i=0;i8;i+) for(j=0;j9;j+) sumpi1=(ei0i*(ybrij*ei0j-ybiij*fi0j)+fi0i*(ybrij*fi0j+ybiij*ei0j); sumpi2+=sumpi1; pi0i=sumpi2; printf(pi0%d=%ft,i,pi0i); sumpi2=0; for(i=0;i8;

8、i+) for(j=0;j9;j+) sumqi1=(fi0i*(ybrij*ei0j-ybiij*fi0j)-ei0i*(ybrij*fi0j+ybiij*ei0j); sumqi2+=sumqi1; qi0i=sumqi2; printf(qi0%d=%ft,i,qi0i); sumqi2=0; for(i=0;i8;i+) detpii=pii-pi0i; detqii=qii-qi0i; if(i=6|i=7) qi0i=ei0i*ei0i+fi0i*fi0i; detqii=1.051-qi0i; printf(detpi%d=%ft,i,detpii); printf(detqi%

9、d=%ft,i,detqii); /* /节点的注入电流表达式 for(i=0;i8;i+) Iii0i=0; Iir0i=0; for(i=0;i8;i+) for(j=0;j9;j+) Iir0i+=ybrij*ei0j-ybiij*fi0j; Iii0i+=ybrij*fi0j+ybiij*ei0j; /*/求解NHJL矩阵for(i=0;i8;i+) for(j=0;j8;j+) if(i=j) if(i=6|i=7) Hij=-ybiij*ei0j+ybrij*fi0j+Iii0i; Nij=ybrij*ei0j+ybiij*fi0j+Iir0i; Jij=2*fi0i; Lij=2

10、*ei0i; else Hij=-ybiij*ei0j+ybrij*fi0j+Iii0i; Nij=ybrij*ei0j+ybiij*fi0j+Iir0i; Jij=-ybrij*ei0j-ybiii*fi0j+Iir0i; Lij=-ybiij*ei0j+ybrij*fi0j-Iii0i; else if(i=6|i=7) Hij=ybrij*fi0j-ybiij*ei0j; Nij=ybrij*ei0j+ybiij*fi0j; Jij=0; Lij=0; else Hij=ybrij*fi0j-ybiij*ei0j; Nij=ybrij*ei0j+ybiij*fi0j; Jij=-ybii

11、j*fi0j-ybrij*ei0j; Lij=ybrij*fi0j-ybiij*ei0j; /形成jacobian矩阵 for(i=0;i16;i+) for(j=0;j16;j+)if(i%2=0&j%2=0) J0ij=Hi/2j/2;else if(i%2=0&j%2!=0) J0ij=Ni/2(j-1)/2;else if(i%2!=0&j%2=0) J0ij=J(i-1)/2j/2;else J0ij=L(i-1)/2(j-1)/2; /for(i=0;i16;i+) /for(j=0;j16;j+) / printf(J0%d%d=%ft,i,j,J0ij); /* /求detui

12、 /* for(i=0;i16;i+)if(i%2=0) detsii=detpii/2; else detsii=detqi(i-1)/2;/将detp和detq用一个数组表示for(i=0;iT;i+) for(j=0;j(2*T);j+) if (j16) aij=J0ij; else if (j=T+i) aij=1.0; else aij=0.0; for(i=0;iT;i+) for(k=0;kT;k+) if(k!=i) t=aki/aii; for(j=0;j(2*T);j+) xx=aij*t; akj=akj-xx; for(i=0;iT;i+) t=aii; for(j=

13、0;j(2*T);j+) aij=aij/t; for(i=0;iT;i+) for(j=0;jT;j+) niij=aij+T; /*printf(逆矩阵为：n); for (i=0;iT;i+) for (j=0;jT;j+) printf(%10.3f,niij); printf(n); */ xxx=0.0;for(i=0;iT;i+)xxx=0.0;for(j=0;jT;j+)xxx=xxx+niij*detsij;detuii=xxx;/检测detui满足要求与否max=Max(detui,16); printf(max=%fn,max); /*/算第n次迭代后的u for(i=0

14、;iT;i+) if(i%2=0) fi1i/2=fi0i/2+detuii; elseei1(i-1)/2=ei0(i-1)/2+detuii; /*for(i=0;i8;i+)/下一次迭代赋初值 ei0i=ei1i; fi0i=fi1i; for(i=0;i8;i+) printf(ei0%d=%ft,i,ei0i); printf(fi0%d=%ft,i,fi0i); for(i=0;ieps&hkm); printf(All do %d timesn,h); sumir=0; sumii=0; /*for(i=0;i9;i+)/平衡节点功率计算 I1r=mulRe(ybr8i,-ybi

15、8i,ei0i,-fi0i); I1i=mulIm(ybr8i,-ybi8i,ei0i,-fi0i); sumir+=I1r; sumii+=I1i; pi8=mulRe(ei08,fi08,sumir,sumii); qi8=mulIm(ei08,fi08,sumir,sumii); printf(S9=%f+j%fn,pi8,qi8); sumpi1=0; sumpi2=0; sumqi1=0; sumqi2=0; for(i=0;i9;i+) for(j=0;j9;j+) if(i!=j&ybiij!=0) sumpi1=mulRe(ei0i,fi0i,0.0,ydij); sumqi1

16、=mulIm(ei0i,fi0i,0.0,ydij); sumpi2=mulRe(ei0i-ei0j,fi0i-fi0j,-ybrij,-ybiij); sumqi2=mulIm(ei0i-ei0j,fi0i-fi0j,-ybrij,-ybiij); sumpi1+=sumpi2; sumqi1+=sumqi2; spij=mulRe(ei0i,fi0i,sumpi1,-sumqi1); sqij=mulIm(ei0i,fi0i,sumpi1,-sumqi1); printf(S%d=%f+j%fn,(i+1)*10+j+1,spij,sqij); for(i=0;i9;i+) for(j=0;j9;j+) if(j !=i&ybiij!=0) delta_pij=spij+spji; delta_qij=sqij+sqji; if( (delta_pij !=0) | (delta_qij !=0) printf(dS%d=%f+j%fn,(i+1)*10+j+1,delta_pij,delta_qij); ei18=ei08; fi18=fi08; for(i=0;i9;i+) printf(u%d=%f%fn,i+1,sqrt(ei1i*ei1i+fi1i*fi1i),atan(fi1i/ei1i)*180/3.14159);