基于matlab实现流计算万能程序，能够计算潮流，网络潮流程序编写的不错.rar

clear; clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %B1来表示支路参数，以下是关于B1矩阵的说明 %B1的第一列表示支路的首字母 %B1的第二列表示支路的末字母 %B1的第三列表示支路的电阻参数 %B1的第四列表示支路的电抗参数 %B1的第五列表示支路的电纳参数;注明：表示电纳参数的时候只是用了其数值 %B1的第六列表示支路的变压器变比;注明：变压器阻抗总是折算到变压器的低压侧，支路的首端与变压器的低压侧阻抗相连； %B1的第七列表示支路是否含有变压器，如果有则为1，否则为0 B1=[1,2,1.099,56.911,0,1.05,1; 2,3,4,32.5,0.0000312,1,0; 2,7,5,40.625,0.000039,1,0; 3,7,3,24.375,0.0000234,1,0; 6,3,1.22,76.963,0,0.956,1; 6,4,0.749,-5.273,0,1.05,1; 6,5,0.791,50.845,0,1,1]; %B2来表示节点参数，以下是关于B2矩阵的说明 %第一列表示该节点电源注入的功率 %第二列表示该节点负荷所需要的功率 %第三列表示该节点电压的初始值 %第四列表示PV节点电压的规定值 %第五列表示该节点的无功补偿容量 %第六列表示该节点类型，1表示平衡节点，2表示PQ节点，3表示PV节点； B2=[1.5+0.9i,0.0019+0.00736i,1,0,0,2; 0,1+0.936i,1,0,0,2; 0,-0.0499i,1,0,0,2; 0,1+0.8i,1,0,0,2; 0,0.2+0.15i,1,0,0,2; 0,0.00075+0.00079i,1,0,0,2; 0,0,1.05,0,0,1]; %n为网络的节点数 n=7; %b为网络支路的条数 b=7; %per为误差的精度 per=0.00001; %ns为不满足误差精度节点的个数，迭代次数用k表示 %PE,PF,QE,QF是PQ节点在雅格比矩阵中的元素 %PE,PF,VE,VF是PV节点在雅格比矩阵中的元素 PE=zeros(n-1,n-1);PF=zeros(n-1,n-1);QE=zeros(n-1,n-1);QF=zeros(n-1,n-1); VE=zeros(n-1,n-1);VF=zeros(n-1,n-1); %定义雅格比矩阵 J=zeros(2*(n-1),2*(n-1)); %雅格比矩阵形成的增广矩阵 JW=zeros(2*(n-1),2*(n-1)+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第一步形成节点导纳矩阵Y Y=zeros(n,n); for i=1:n if B1(i,n)==0 %the road don't have transformer p=B1(i,1); q=B1(i,2); Y(p,q)=-(Y(p,q)+1/(B1(i,3)+sqrt(-1)*B1(i,4))); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+B1(i,5)*sqrt(-1)+1/(B1(i,3)+sqrt(-1)*B1(i,4)); Y(q,q)=Y(q,q)+B1(i,5)*sqrt(-1)+1/(B1(i,3)+sqrt(-1)*B1(i,4)); else %the road have transformer p=B1(i,1); q=B1(i,2); Y(p,q)=-(Y(p,q)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2))); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+B1(i,5)*sqrt(-1)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2)); Y(q,q)=Y(q,q)+B1(i,5)*sqrt(-1)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2)); end end disp('得出网络的节点导纳矩阵Y'); disp(Y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第二步做一些准备工作，为求出功率的偏移量做准备,这些做准备的偏移量都是以多行一列的形式表达 %1）分解出节点导纳矩阵的实部G和虚部B G=real(Y); B=imag(Y); %求出各节点输入输出功率的初值；注明：流入网络的功率为正 %各节点注入的有功为P,注入的无功为Q;注明此时得出的无功功率是指数值大小，不含有i; P=zeros(n-1,1); Q=zeros(n-1,1); for i=1:n-1 P(i,1)=real(B2(i,1))-real(B2(i,2)); Q(i,1)=imag(B2(i,1))-imag(B2(i,2))+B(i,5); end disp(P); disp(Q); %2）先定义原始的PV电压值，再分解出电压值的实部e和虚部f,这里要特别清楚复数实部和虚部的定义 for i=1:n V(i,1)=B2(i,n-3); %PV节点电压的初值 end disp(V); E=zeros(n,1); F=zeros(n,1); for i=1:n if B2(i,n-1)==3 e(i,1)=real(B2(i,n-3)); f(i,1)=imag(B2(i,n-3)); else e(i,1)=real(B2(i,n-4)); f(i,1)=imag(B2(i,n-4)); end end disp(e); disp(f); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第三步奏开始判断且求出第一次迭代量 %定义迭代次数k k=0; %定义含有不符合误差精度的节点个数ns,先假设ns=1; ns=1; %由于下面的程序必须涉及到功率偏移量或者电压的偏移量 %定义有功偏移量XP,无功偏移量XQ,电压的偏移量XV,为了万能程序的形成，此处先不定义明确 XP=zeros(n-1,1); XQ=zeros(n-1,1); XV=zeros(n-1,1); %%%%%%%%%%%%%!!!!说明由于要反复迭代计算，里面的元素不能够赋值为0！！！！！！！！ while ns~=0; ns=0; %先赋值为0，为后面的判断做铺垫 k=0; for i=1:n-1 for j=1:n XP(i,1)=XP(i,1)+e(i,1)*(G(i,j)*e(j,1)-B(i,j)*f(j,1))+f(i,1)*(G(i,j)*f(j,1)+B(i,j)*e(j,1)); XQ(i,1)=XQ(i,1)+f(i,1)*(G(i,j)*e(j,1)-B(i,j)*f(j,1))-e(i,1)*(G(i,j)*f(j,1)+B(i,j)*e(j,1)); XV(i,1)=XV(i,1)+e(i,1)^2+f(i,1)^2; end if B(i,n-1)==3 XV(i,1)=V(i,1)^2-XV(i,1); else XP(i,1)=P(i,1)- XP(i,1); XQ(i,1)=Q(i,1)- XQ(i,1); end end disp('功率的偏移量') disp(XP); disp(XQ); disp(XV); %XV的结果都市7； %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %做几点说明：本来在此处判断三个偏移量是否满足收敛要求，但是为了程序的连续性，先求雅格比矩阵求出电压的偏移量再做判断，并且程序在此处一般不能满足收敛要求%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第四部求雅格比矩阵 for i=1:(n-1) for j=1:(n-1) if i~=j PE(i,j)=-(G(i,j)*e(i,1)+B(i,j)*f(i,1)); QF(i,j)=-PE(i,j); PF(i,j)=B(i,j)*e(i,1)-G(i,j)*f(i,1); PF(i,j)=QE(i,j); VE(i,j)=0; VF(i,j)=0; else for x=1:n PE(i,j)=PE(i,j)-(G(i,x)*e(x,1)-B(i,x)*f(x,1)); PF(i,j)=PF(i,j)-(G(i,x)*f(x,1)+B(i,x)*e(x,1)); QE(i,j)=QE(i,j)+(G(i,x)*f(x,1)+B(i,x)*e(x,1)); QF(i,j)=QF(i,j)-(G(i,x)*e(x,1)-B(i,x)*f(x,1)); VE(i,i)=-2*e(i,1); VF(i,i)=-2*f(i,1); end PE(i,j)=PE(i,j)-G(i,i)*e(i,1)-B(i,i)*f(i,1); PF(i,j)=PF(i,j)+B(i,i)*e(i,1)-G(i,i)*f(i,1); QE(i,j)=QE(i,j)+B(i,i)*e(i,1)-G(i,i)*f(i,1); QF(i,j)=QF(i,j)+G(i,i)*e(i,1)+B(i,i)*f(i,1); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:2*(n-1) for j=1:2*(n-1) if rem(i,2)~=0 %奇数行 if rem(j,2)~=0 %为奇数行奇数列 a=(i+1)/2; b=(j+1)/2; J(i,j)=PE(a,b); else %为奇数行偶数列 a=(i+1)/2; b=j/2; J(i,j)=PF(a,b); end else %偶数行 a=i/2; if rem(j,2)~=0 %偶数行奇数列 b=(j+1)/2; if B2(a,n-1)==3 %如果是PV节点的 J(i,j)=VE(a,b); else %如果是PQ节点 J(i,j)=QE(a,b); end else %偶数行偶数列 b=j/2; if B2(a,n-1)==3 %如果是PV节点的 J(i,j)=VF(a,b); else %如果是PQ节点 J(i,j)=QF(a,b); end end end end end disp('输出雅格比矩阵') disp(J); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第五步雅格比矩阵和功率的偏移量形成增广矩阵JW %说明这里的增广矩阵是定义为JW,在形成时乘以*（-1） J=(-1)*J; for i=1:12 for j=1:13 if j<=12 JW(i,j)=J(i,j); else if rem(i,2)~=0 %不论是PQ,PV节点，奇数行都是一样的 a=(i+1)/2; JW(i,j)=XP(a,1); else a=i/2; if B2(a,n-1)==3 %偶数行为PV节点 JW(i,j)=XV(a,1); else %偶数行为PQ节点 JW(i,j)=XQ(a,1); end end end end end disp('输出增广矩阵JW') disp(JW); %%%%%%%%%%%%%%%%%%%%%%%%%