基于MATLAB的电力系统潮流计算.doc

基于MATLAB的电力系统潮流计算%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号，将变压器的串联阻抗置于低压侧处理。%第三列为支路的串列阻抗参数。%第四列为支路的对地导纳参数。%第五烈为含变压器支路的变压器的变比%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地%参数。n=input(请输入节点数:n=);n1=input(请输入支路数:n1=);isb=input(请输入平衡节点号:isb=);pr=input(请输入误差精度:pr=);B1=input(请输入支路参数:B1=);B2=input(请输入节点参数:B2=);X=input(节点号和对地参数:X=);Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵for i=1:n1 if B1(i,6)=0 %不含变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3); Y(q,q)=Y(q,q)+1/(B1(i,5)2*B1(i,3); endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendOrgS%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endendPVU%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); endendt=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理 if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendk=0;for i=1:n %对PV节点的处理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理 if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n %对PV节点的处理 if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2%开始循环*while abs(max(DetaU)prOrgS=zeros(2*n-2,1); %!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendOrgS%创建DetaSh=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); endendt=0;for i=1:n if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS%创建Ii=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) k=0; end end end endendk=0;for i=1:n if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) k=0; end end end endendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度pr 0.000001B1(支路参数矩阵)1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1B2(节点参数矩阵)0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1X(节点号和对地参数)1 0;2 0;3 0;4 0;5 0电力系统潮流计算9结点算例-PQ法原始数据录入data.txt文档:标号,起始结点,终止结点,支路电阻参数,支路电抗参数,支路对地导纳参数1,2,5,0.0,0.063,0.0,2,5,9,0.019,0.072,0.075,3,6,9,0.012,0.101,0.105, 4,3,6,0.0,0.059,0.0,5,6,8,0.039,0.17,0.179,6,4,8,0.017,0.092,0.079, 7,5,7,0.032,0.161,0.153,8,4,7,0.01,0.085,0.088,9,1,4,0.0,0.058,0.0,潮流程序chaoliu.txt文档:#include#include#define N 9 /*总结点数*/#define M 6 /*PQ结点数*/#define K 9 /*线路数*/#define eps 1e-4void guass(int n,int m,float c,float bN,float x) /*高斯函数*/float aNN,yN;int i,j,k;for(i=0;im;i+) yi=ci; for(j=0;jm;j+) aij=bi+nj+n;for(k=0;km-1;k+) for(i=k+1;im;i+) for(j=k+1;j=0;i-) for(j=i+1;jm;j+) yi-=aij*xj; xi=yi/aii;struct lineint Lindex; int Headnode; int Endnode; float R; float X; float b;lineK;struct line *t;main() float rN=0.0; float uN=1.04,1.025,1.025,1.0,1.0,1.0,1.0,1.0,1.0; float pN=1,1.63,0.85,0.0,0.0,0.0,-1.25,-0.9,-1.0; float qN=1,1,1,0.0,0.0,0.0,-0.5,-0.3,-0.35; float gNN=0.0,bNN=0.0; float hNN=0.0; float BNN=0.0; float temp; float HK6; float lr,lx,lb1,lg,lb; int i,j; int ku=0,kr=0,kp=1,kq=1; void val(float uN,float gNN,float bNN,float rN,int ku, int kr,float hNN); /*函数申明*/ FILE *fp; fp=fopen(data.txt,r); for(i=0;iK;i+) for(j=0;j6;j+) fscanf(fp,%f,&temp); Hij=temp; fclose(fp); for(i=0;iK;i+) linei.Lindex=(int)Hi0;linei.Headnode=(int)Hi1;linei.Endnode=(int)Hi2; linei.R=Hi3;linei.X=Hi4;linei.b=Hi5; for(t=line;tHeadnode-1; j=t-Endnode-1; lr=t-R; lx=t-X; lb1=t-b; lg=lr/(lr*lr+lx*lx); lb=-lx/(lr*lr+lx*lx); gii+=lg; gjj+=lg; bii+=lb+lb1; bjj+=lb+lb1; hij=hji=-lb1; gij-=lg; gji-=lg; bij-=lb; bji-=lb; getch(); printf(n=jie dian dao na ju zhen=n); for(i=0;iN;i+) for(j=0;jN;j+) Bij=bij; printf(%8f,Bij); printf(n); printf(n); getch(); printf(n=gei ding chu zhi=n); for(i=0;iN;i+) printf(u%d=%8f p%d=%8f q%d=%8fn,i+1,ui,i+1,pi,i+1,qi); printf(n=die dai qiu jie=n); while(kp=1) float ip,iq,max; float dpN,dqN,drN; float dpuN,dquN; float y1N-1,y2 ; float x1N-1,x2 ; for(i=1;iN;i+) /*算dp对应B,N-1维,除去平衡结点*/ ip=0; for(j=0;jN;j+) ip=ip+uj*(gij*cos(ri-rj)+bij*sin(ri-rj); dpi=pi-ui*ip; dpui=dpi/ui; printf(dp%d=%8fn,i+1,dpi); printf(n); getch(); max=fabs(dpu1); for(i=1;imax) max=fabs(dpui); if (max=eps) for(i=0;iN-1;i+) y1i=-dpui+1; /*起值不同,为了对应,故加一*/ guass(1,N-1,y1,B,x1); for (i=1;iN;i+) dri=x1i-1/ui; ri=ri+dri; printf(n= di %d ci die dai hou dian ya xiang jiao chu zhi =n,kr+1); for(i=1;iN;i+) printf(r%d=%8fn,i+1,ri); getch(); printf(nn); kr=kr+1;kq=1; top: for(i=N-M;iN;i+) /*算dq对应B,仅M维,除去平衡结点和PV结点*/ iq=0; for(j=0;jN;j+) iq=iq+uj*(gij*sin(ri-rj)-bij*cos(ri-rj); dqi=qi-ui*iq; dqui=dqi/ui; printf(dq%d=%8fn,i+1,dqi); printf(n); max=fabs(dquN-M); for (i=N-M;imax) max=fabs(dqui); if(max=eps) for(i=0;iM;i+) y2i=-dqui+N-M; /*同上,对应加N-M*/ guass(N-M,M,y2,B,x2); for(i=N-M;iN;i+) ui=ui+x2i-(N-M); printf(n=di %d ci die dai dian ya chu zhi=n,ku+1); for(i=N-M;iN;i+) printf(u%d=%8fn,i+1,ui); printf(nn); ku=ku+1;kp=1; else kq=0; if(kp=0) val(u,g,b,r,ku,kr,h); else kp=0; if(kq=0) val(u,g,b,r,ku,kr,h); else goto top; void val(float uN,float gNN,float bNN,float rN,int ku, int kr,float hNN) float ps=0,pv1=0,pv2=0; float qs=0,qv1=0,qv2=0; float pNN=0; float qNN=0; float sNN; float dpNN=0; float dqNN=0; float dsNN; float dSp=0,dSq=0; int i,j; FILE *fp1; printf(n=ping heng jie dian gong lv =n); getch(); for(i=0;iN;i+) ps=ps+u0*ui*(g0i*cos(r0-ri)+b0i*sin(r0-ri); qs=qs+u0*ui*(g0i*sin(r0-ri)-b0i*cos(r0-ri); printf(sp=%8f+j(%8f)n,ps,qs); printf(n=PV jie dian gong lv=n); getch(); for(i=0;iN;i+) pv1=pv1+u1*ui*(g1i*cos(r1-ri)+b1i*sin(r1-ri); qv1=qv1+u1*ui*(g1i*sin(r1-ri)-b1i*cos(r1-ri); printf(sv1=%8f+j(%8f)n,pv1,qv1); for(i=0;iN;i+) pv2=pv2+u2*ui*(g2i*cos(r2-ri)+b2i*sin(r2-ri); qv2=qv2+u2*ui*(g2i*sin(r2-ri)-b2i*cos(r2-ri); printf(sv2=%8f+j(%8f)n,pv2,qv2); for(i=0;iN;i+) for(j=0;jN;j+) pij=ui*ui*(-gij)+ui*uj*(gij*cos(ri-rj)+bij*sin(ri-rj); qij=ui*ui*(-hij+bij)+ui*uj*(gij*sin(ri-rj)-bij*cos(ri-rj); dpij=ui*ui*(-gij)+ui*uj*(gij*cos(ri-rj)+bij*sin(ri-rj) +uj*uj*(-gji)+uj*ui*(gji*cos(rj-ri)+bji*sin(rj-ri); dqij=ui