neidianfa.rar_5节点潮流计算_5节点系统潮流_内点法_内点法matlab_直角潮流计算

clear; %clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %数据加载 n=input('请输入要计算的节点系统(5):') load Node5.txt;%节点数据 load Branch5.txt;%支路数据 load Generator5.txt;%发电机数据 Node=Node5; Branch=Branch5; Generator=Generator5; %节点数据处理 N=Node(:,1);%节点号 Type=Node(:,2);%节点类型 Uamp=Node(:,3);%节点电压幅值 Dlta=Node(:,4);%节点电压相角 Pd=Node(:,5);%节点负荷有功 Qd=Node(:,6);%节点负荷无功 Pg=Node(:,7);%节点出力有功 Qg=Node(:,8);%节点出力无功 Umax=Node(:,9);%节点电压幅值上限 Umin=Node(:,10);%节点电压幅值下限 Bc=Node(:,11);%节点补偿电容电纳值 %支路数据处理 Nbr=Branch(:,1);%支路号 Nl=Branch(:,2);%支路首节点 Nr=Branch(:,3);%支路末节点 R=Branch(:,4);%支路电阻 X=Branch(:,5);%支路电抗 Z=R+1i*X;%支路阻抗=支路电阻+支路电抗 Bn=Branch(:,6);%支路对地电纳 K=Branch(:,7);%支路变压器变比，0表示无变压器 Ptmax=Branch(:,8);%线路传输功率上限 %发电机数据处理 Ng=Generator(:,1);%发电机序号 Nbus=Generator(:,2);%所在母线号 Pumax=Generator(:,3);%发电机有功出力上界 Qumax=Generator(:,4);%发电机无功出力上界 Pumin=Generator(:,5);%发电机有功出力下界 Qumin=Generator(:,6);%发电机无功出力下界 a2=Generator(:,7);%燃料耗费曲线二次系数 a1=Generator(:,8);%燃料耗费曲线一次系数 a0=Generator(:,9);%燃料耗费曲线常数项 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n=length(N);%节点个数 ng=length(Ng);%发电机台数 nbr=length(Nbr);%支路个数 x=zeros(2*(ng+n),1);%控制变量+状态变量 x(1:ng)=Pg(Nbus); x(ng+1:2*ng)=Qg(Nbus); x((2*ng+2):2:2*(ng+n))=Uamp; x((2*ng+1):2:2*(ng+n)-1)=Dlta; l=0.8*ones(2*ng+n+nbr,1);%松弛变量 u=1.1*ones(2*ng+n+nbr,1);%松弛变量 w=-1.5*ones(2*ng+n+nbr,1);%拉格朗日乘子 z=ones(2*ng+n+nbr,1);%拉格朗日乘子 y=zeros(2*n,1);%拉格朗日乘子 y(1:2:2*n-1)=1e-3; y(2:2:2*n)=-1e-3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算不等式约束的上下限%%%%%%%%%%%%%%%%%%%%%%%% %gmin gmin=zeros(2*ng+n+nbr,1); gmin(1:ng)=Pumin; gmin(ng+1:2*ng)=Qumin; gmin(2*ng+1:2*ng+n)=Umin; gmin(2*ng+n+1:2*ng+n+nbr)=-Ptmax; %gmax gmax=zeros(2*ng+n+nbr,1); gmax(1:ng)=Pumax; gmax(ng+1:2*ng)=Qumax; gmax(2*ng+1:2*ng+n)=Umax; gmax(2*ng+n+1:2*ng+n+nbr)=Ptmax; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Y=zeros(n,n); %%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%% for ii=1:nbr if K(ii)==0%非变压器支路 Y(Nl(ii),Nr(ii))=-1/Z(ii); Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii)); else%变压器支路 Y(Nl(ii),Nr(ii))=-1/Z(ii)/K(ii); Y(Nr(ii),Nl(ii))= Y(Nl(ii),Nr(ii)); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%计算对角元素%%%%%%%%%%%%%%%%%%%%%% for ii=1:n%将支路导纳加入到对角元素中 for jj=1:nbr if K(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非变压器支路 Y(ii,ii)=Y(ii,ii)+1/Z(jj); else if K(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%变压器支路 Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj); end end end end for ii=1:nbr%将对地电纳加入到对角元素中 if K(ii)==0%非变压器支路 Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+1i*Bn(ii); Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+1i*Bn(ii); else%变压器支路 Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+(K(ii)-1)/K(ii)/Z(ii); Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+(1-K(ii))/K(ii)/K(ii)/Z(ii); end end for ii=1:n Y(ii,ii)=Y(ii,ii)+i*Bc(ii); end G=real(Y);%电导 B=imag(Y);%电纳 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k=0;%迭代次数 Kmax=150;%最大迭代次数 iteration=1e-4;%误差精度 delta=0.08; Gap=(l'*z-u'*w)*ones(Kmax,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while k<50 %计算互补间隙Gap Gap(k+1)=l'*z-u'*w; if Gap>iteration miu=delta*Gap(k+1)/(2*(2*ng+n+nbr)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%相角差计算%%%%%%%%%%%%%%%%%%%%%% theta=zeros(n,n); for ii=1:n for jj=1:n theta(ii,jj)=Dlta(ii)-Dlta(jj); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%1、等式约束雅克比矩阵%%%%%%%%%%%%%%%% pxh=zeros(2*(ng+n),2*n); %%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%% for ii=1:ng pxh(Ng(ii),2*Nbus(ii)-1)=1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%% for ii=1:ng pxh(Ng(ii)+ng,2*Nbus(ii))=1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%ah/ax%%%%%%%%%%%%%%%%%%%%%%% HH=zeros(n,n); JJ=zeros(n,n); NN=zeros(n,n); LL=zeros(n,n); for ii=1:n for jj=1:n if ii~=jj%i!=j时的情况 %非对角元素 HH(ii,jj)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); JJ(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); NN(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); LL(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); %对角元素 HH(ii,ii)=HH(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); JJ(ii,ii)=JJ(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); NN(ii,ii)=NN(ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj))); LL(ii,ii)=LL(ii,ii)-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); end end NN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii); LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii); end pxh(1+2*ng:2:2*(n+ng)-1,1:2:2*n-1)=HH'; pxh(1+2*ng:2:2*(n+ng)-1,2:2:2*n)=JJ'; pxh(2+2*ng:2:2*(n+ng),1:2:2*n-1)=NN'; pxh(2+2*ng:2:2*(n+ng),2:2:2*n)=LL'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%2、不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%% %g1:电源有功出力上下限约束 ag1aP=eye(ng,ng); ag1aQ=zeros(ng,ng); ag1ax=zeros(2*n,ng); %g2:电源无功出力上下限约束 ag2aP=zeros(ng,ng); ag2aQ=eye(ng,ng); ag2ax=zeros(2*n,ng); %g3:节点电压幅值上下限约束 ag3aP=zeros(ng,n); ag3aQ=zeros(ng,n); ag3ax=zeros(2*n,n); for ii=1:n ag3ax(2*ii,ii)=1; end %g4:线路潮流上下限约束 ag4aP=zeros(ng,nbr); ag4aQ=zeros(ng,nbr); ag4ax=zeros(2*n,nbr); for ii=1:n for jj=1:nbr if Nl(jj)==ii ag4ax(2*ii-1,jj)=-Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))); ag4ax(2*ii,jj)=Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))))-2*Uamp(Nl(jj))*G(Nl(jj),Nr(jj)); end if Nr(jj)==ii ag4ax(2*ii-1,jj)=Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))); ag4ax(2*ii,jj)=Uamp(Nl(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj)))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj)))); end end end pxg=[ag1aP ag2aP ag3aP ag4aP; ag1aQ ag2aQ ag3aQ ag4aQ; ag1ax ag2ax ag3ax ag4ax];%此即为不等式约束的雅克比矩