九节点系统潮流计算编程牛N-R法.pdf
如图所示系统，试计算潮流分布，相关数据见 版潮流计算用户手册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 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,sumpi1,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=,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;static float ybi99=,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,0,0,0,-16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,;static float yd99=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;float ei09=,;float fi09=;float pi9=0,0,0,0;float qi9=0,0,0,0,0,0;h=0;km=15;eps=;do h+=1;printf(nNow The%dth timesn,h);for(i=0;i8;i+)printf(ei0%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;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=i;printf(detpi%d=%ft,i,detpii);printf(detqi%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*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=-ybiij*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 /*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=;else aij=;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=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=;for(i=0;iT;i+)xxx=;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;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,-ybi8i,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,ydij);sumqi1=mulIm(ei0i,fi0i,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/;
