各种潮流计算方法的matlab代码

%% 牛顿拉夫逊法 %% function [count,B,Q]=newtpf(n1,n,isb,m,pr) %相关原始数据输入格式如下： %B1是支路参数矩阵：列数1、2为节点编号，编号从小到大编写，从低压侧到高压侧编写；3为支路的串联阻抗参数； %4为支路的对地导纳参数；5为含变压器支路的变压器变比；6为支路首端置于K侧为1,1侧为0。 %B2是节点参数矩阵：1.节点注入发电机功率参数；2.节点负荷功率参数；3.节点电压初始值； %4.PV节点电压给定值；5.节点相连的无功补偿设备容量；6.节点参数类型。 % 求节点导纳矩阵Y n1=input('请输入支路数n1='); % n1为支路数=5 n=input('请输入节点数n='); % n为节点数=5 isb=input('请输入平衡节点号，isb='); %平衡节点isb=5 m=input('请输入PQ节点数,m='); %PQ节点数m=3 pr=input('请输入误差精度，pr='); Y=zeros(n1); %节点导纳矩阵 B1=[1 2 0.04+0.25i 0.25i 1 0; 1 3 0.10+0.35i 0 1 0; 2 3 0.08+0.30i 0.25i 1 0; 2 4 0+0.015i 0 1.05 0; 3 5 0+0.03i 0 1.05 1]; %% 求节点导纳矩阵Y for i=1:n1 %支路数 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); if B1(i,5)==1 Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4); Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); else Y(p,p)= Y(p,p)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4); %对角元K侧 Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4); %对角元1侧 end end Y % 设置各节点电压初值 B2=[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; 5 5+0i 1.05 1.05 0 3; 0 0 1.05 1.05 0 1]; % 建立雅克比矩阵 G=real(Y);B=imag(Y); U=zeros(n,1);thita=zeros(n,1);deltaThita=zeros(n,n); E=zeros(n,n);F=zeros(n,n); J=zeros(7,7); JH=zeros(n-1,n-1);JN=zeros(n-1,m);JK=zeros(m,n-1);JL=zeros(m,m); %初始化雅克比矩阵中的四个部分 Ui=zeros(n,1);Uj=zeros(n,1); % 牛顿迭代 Pi=real(B2(:,2));Qi=imag(B2(:,2)); P1=zeros(n,1);Q1=zeros(n,1); DeltaU=ones(m,1);DeltaThita=ones(n-1,1); count=0; %迭代次数 deltaP=zeros(n-1,1); deltaQ=zeros(m,1); c=zeros(n,1);d=zeros(n,1); Resu=zeros(7,1); TXdeltaU=zeros(1,10); while max(abs([DeltaThita;DeltaU]))>pr %判断收敛性 count=count+1; U=B2(:,3); for i=1:n for j=1:n E(i,j)=G(i,j)*sin(thita(i)-thita(j))-B(i,j)*cos(thita(i)-thita(j)); %(G*sin-B*cos) F(i,j)=G(i,j)*cos(thita(i)-thita(j))+B(i,j)*sin(thita(i)-thita(j)); %(G*cos+B*sin) end end c(:,1)=0; d(:,1)=0; for i=1:n for j=1:n c(i,1)=c(i,1)+F(i,j)*B2(j,3); %求和∑F d(i,1)=d(i,1)+E(i,j)*B2(j,3); %求和∑E end end for i=1:n-1 deltaP(i,1)=Pi(i,1)-c(i,1)*B2(i,3); %计算节点功率△P if i<=m deltaQ(i,1)=Qi(i,1)-d(i,1)*B2(i,3); %计算节点功率△Q end end %迭代雅克比矩阵 Ui=B2(:,3); Uj=B2(:,3); for i=1:n for j=1:n if i==j %对角元 JH(i,j)=Ui(i,1)*(d(i,1)+Uj(j,1)*B(i,j)); JN(i,j)=-(c(i,1)-Uj(j,1)*G(i,j))-2*Ui(i,1)*G(i,i); JK(i,j)=-Ui(i,1)*(c(i,1)-Uj(j,1)*G(i,j)); JL(i,j)=-(d(i,1)+Uj(j,1)*B(i,j))+2*Ui(i,1)*B(i,i); else %非对角元 JH(i,j)=-Ui(i,1)*E(i,j)*Uj(j,1); JN(i,j)=-F(i,j)*Ui(i,1); JK(i,j)=Ui(i,1)*F(i,j)*Uj(j,1); JL(i,j)=-E(i,j)*Ui(i,1); end end end J=[JH JN;JK JL]; J([5 9 10],:)=[]; %去除不需要的行 J(:,[5 9 10])=[]; %去除不需要的列 %% 求解修正方程式，得到修正量 Resu=(-J)\[deltaP;deltaQ]; DeltaThita=Resu(1:4,1); DeltaU=Resu(5:7,1); TXdeltaU(1,count)=-DeltaU(1,1); %% 修正节点电压 for i=1:n if B2(i,6)~=1 %不是平衡节点 thita(i,1)=thita(i,1)+DeltaThita(i); if B2(i,6)~=3 B2(i,3)=B2(i,3)+DeltaU(i); end end end if count==0 J; break; end end J count DeltaUU(1,:)=DeltaU(:,1); DeltaUU DeltaThita deltaP deltaQ %% 收敛后，计算各个计算量 P=zeros(n,1);Q=zeros(n,1); sum=0; for i=1:n if B2(i,6)==1 %平衡节点的有功和无功公式 for k=1:n sum=sum+conj(Y(i,k))*conj(U(k)); end P(i,1)=real(B2(i,3)*sum); Q(i,1)=imag(B2(i,3)*sum); end if B2(i,6)==2 %PQ节点 P(i,1)=real(B2(i,2)); Q(i,1)=imag(B2(i,2)); end if B2(i,6)==3 %PV节点 P(i,1)=5; for j=1:n Q(i,1)=B2(i,3)*c(i,1); end end end P Q thita=thita*180/pi; UU=B2(:,3)'; UU Thita=thita'; Thita TXdeltaU subplot(2,1,1) plot(TXdeltaU(1,:)),xlabel('迭代次数'),ylabel('电压不平衡量最大值'),title('收敛特性曲线'); axis([0,6,0,0.1]); set(gca,'XTick',0:1:6); set(gca,'YScale','linear'); hold on;