潮流计算程序至雅克比矩阵

%本程序的功能是用牛顿——拉夫逊法进行潮流计算%X1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳%5、支路的变比；6、支路首端处于K侧为1，1侧为0%X2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值%4、PV节点电压V的给定值；5、初相角给定值6、节点所接的无功补偿设备的容量%7、节点分类标号：1为PQ节点；2为PV节点；3为平衡节点；clear;n=6;%input('请输入节点数：n=')nl=7;%input('请输入支路数：nl=')n2=5;%input('请输入已知条件中含有有功功率的节点数：n2=')n3=4;%input('请输入已知条件中含有无功功率的节点数：n3=')isb=6;%input('请输入平衡母线节点号：isb=')pr=0.00001;%input('请输入误差精度：pr=')X1=[120.000+0.300i01.0251;140.097+0.407i010;160.123+0.518i010;250.282+0.640i010;350.723+1.050i010;340.000+0.133i01.1000;460.080+0.370i010]%input('请输入由支路参数形成的矩阵：X1=')X2=[050.0+5.0i1.0000001;030.0+18.0i1.0000001;055.0+13.0i1.0000001;001.0000001;50.101.0001.100002;001.0001.050003]%input('请输入各节点参数形成的矩阵：X2=')X3=[10;20;30;40;50;60]%input('请输入由节点号及其对地阻抗形成的矩阵：X3=')Y=zeros(n);u=zeros(1,n);v=zeros(1,n);delt=zeros(1,n);JJ=zeros(n2+n3);%设置各个参数矩阵的形式fori=1:nifX3(i,2)~=0;a=X3(i,1);Y(a,a)=1./X3(i,2);endendfori=1:nlifX1(i,6)==0a=X1(i,1);b=X1(i,2);elsea=X1(i,2);b=X1(i,1);endY(a,b)=Y(a,b)-1./(X1(i,3)*X1(i,5));Y(b,a)=Y(a,b);Y(b,b)=Y(b,b)+1./(X1(i,3)*X1(i,5)^2)+X1(i,4)./2;Y(a,a)=Y(a,a)+1./X1(i,3)+X1(i,4)./2;end%求导纳矩阵disp('导纳矩阵Y=')disp(Y)G=real(Y);B=imag(Y);%将导纳矩阵的实部和虚部整理出来fori=1:nu(i)=X2(i,3);v(i)=X2(i,4);delt(i)=X2(i,5);end%根据X2矩阵得出各个节点电压的初始给定值和相角fori=1:nS(i)=X2(i,1)-X2(i,2);B(i,i)=B(i,i)+X2(i,6);end%计算各个节点的发电机功率与负荷功率的差值，以及检查各个节点是否有无功补偿装置P=real(S);Q=imag(S);%将各个节点的发电机功率与负荷功率的差值的实部和虚部整理出来disp(P);disp(Q);k=0;%设置迭代次数初始值pr=1;%给定一个误差精度whilepr0.00001fora=1:n2forb=1:npt(a)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));endpp(a)=P(a)-sum(pt);endfora=1:n3forb=1:nqt(a)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));endqq(a)=Q(a)-sum(qt);end%求得各个节点的有功功率和无功功率的误差disp(qq);disp(pp);fora=1:n2forb=1:nhh(b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));endw1(a)=sum(hh);fora=1:n3forb=1:nnn(b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));kk(b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));ll(b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));endw2(a)=sum(nn);w3(a)=sum(kk);w4(a)=sum(ll);end%求式（3-33）中的包括j=i在内所计算的结果fora=1:n2forb=1:n2ifa==bH(a,b)=w1(a)-u(a)^2*(G(a,a)*sin(delt(a)-delt(a))-B(a,a)*cos(delt(a)-delt(a)));elseH(a,b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));endendend%计算H单元的元素fora=1:n2forb=1:n3ifa==bN(a,b)=w2(a)+2*u(a)^2*G(a,a)+u(a)^2*(G(a,a)*cos(delt(a)-delt(a))+B(a,a)*sin(delt(a)-delt(a)));elseN(a,b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));endendend%计算N单元的元素fora=1:n3forb=1:n2ifa==bK(a,b)=w3(a)+u(a)^2*(G(a,a)*cos(delt(a)-delt(a))-B(a,a)*sin(delt(a)-delt(a)));elseK(a,b)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));endendend%计算K单元的元素fora=1:n3forb=1:n3ifa==bL(a,b)=w4(a)+2*u(a)^2*G(a,a)+u(a)^2*(G(a,a)*sin(delt(a)-delt(a))+B(a,a)*cos(delt(a)-delt(a)));elseL(a,b)=-u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));endendend%计算L单元的元素JJ=[H(1,1)H(1,2)H(1,3)H(1,4)H(1,5)N(1,1)N(1,2)N(1,3)N(1,4);H(2,1)H(2,2)H(2,3)H(2,4)H(2,5)N(2,1)N(2,2)N(2,3)N(2,4);H(3,1)H(3,2)H(3,3)H(3,4)H(3,5)N(3,1)N(3,2)N(3,3)N(3,4);H(4,1)H(4,2)H(4,3)H(4,4)H(4,5)N(4,1)N(4,2)N(4,3)N(4,4);H(5,1)H(5,2)H(5,3)H(5,4)H(5,5)N(5,1)N(5,2)N(5,3)N(5,4);K(1,1)K(1,2)K(1,3)K(1,4)K(1,5)L(1,1)L(1,2)L(1,3)L(1,4);K(2,1)K(2,2)K(2,3)K(2,4)K(2,5)L(2,1)L(2,2)L(2,3)L(2,4);K(3,1)K(3,2)K(3,3)K(3,4)K(3,5)L(3,1)L(3,2)L(3,3)L(3,4);K(4,1)K(4,2)K(4,3)K(4,4)K(4,5)L(4,1)L(4,2)L(4,3)L(4,4)]%排列雅可比矩阵的位置disp(JJ)fora=1:n2DP(a,1)=pp(a);endfora=1:4DP(a+n2,1)=qq(a);end%放置有功功率和无功功率误差的位置uu=-inv(JJ)*DP;%inv表示JJ的逆矩阵pr=max(abs(uu));%检验是否满足要求fora=1:n2delt(a)=delt(a)+uu(a);endfora=1:n3u(a)=u(a)+uu(a+n2);end%对方程式进行修正k=k+1;k-1,delt',u'end