Power0704—较理想(1).zip_潮流计算_环网 潮流_环网 潮流计算_环网潮流_环网潮流计算

clear; clc; tic; [ENodes] =xlsread('wokao.xlsx','C3:H7');%节点信息 [EBrchs] =xlsread('wokao.xlsx','I3:R7');%支路信息 ENodeTYPE=ENodes(:,2); %节点类型 U =ENodes(:,3); %电压幅值(为PQ赋电压初值，直接指定PV电压值) Ps =ENodes(:,4); %有功功率 Qs =ENodes(:,5); %无功功率 Qsh =ENodes(:,6); %并联电容器的电纳 ENodeNum =size(ENodeTYPE,1); %节点数目 alpha =zeros(ENodeNum,1); %电压相角 Linei =EBrchs(:,1); %支路首端节点 Linej =EBrchs(:,2); %支路末端节点 Rij =EBrchs(:,3); %支路电阻 Xij =EBrchs(:,4); %支路电抗 B0 =EBrchs(:,5); %线路电纳 Zij =Rij+1i*Xij; KT =EBrchs(:,6); %变压器非标准变比 EBrchNum =size(Linei,1); %支路数目 %% 求节点导纳矩阵Y m_lineY=1./(KT.*Zij); Y=sparse(Linei,Linej,-m_lineY,ENodeNum,ENodeNum)+ sparse(Linej,Linei,-m_lineY,ENodeNum,ENodeNum); %所有支路互导纳 Y=Y+sparse(Linei(KT==1),Linei(KT==1),m_lineY(KT==1)+1i*B0(KT==1)/2+1i*Qsh(KT==1),ENodeNum,ENodeNum)+... sparse(Linej(KT==1),Linej(KT==1),m_lineY(KT==1)+1i*B0(KT==1)/2+1i*Qsh(KT==1),ENodeNum,ENodeNum)+... sparse(Linei(KT~=1),Linei(KT~=1),m_lineY(KT~=1)+(1-KT(KT~=1))/(KT(KT~=1)^2)*(1/Zij(KT~=1))+1j*B0(KT~=1)/2+1j*Qsh(KT~=1),ENodeNum,ENodeNum)+... sparse(Linej(KT~=1),Linej(KT~=1),m_lineY(KT~=1)+(KT(KT~=1)-1)/KT(KT~=1)*(1/Zij(KT~=1))+1j*B0(KT~=1)/2+1j*Qsh(KT~=1), ENodeNum,ENodeNum); %% NR潮流计算 Det = ones(2*ENodeNum,1); %误差初值 KK=0; %迭代次数 PV = find(ENodeTYPE==2); %PV节点号 要求节点号必须顺序排列，这样位置即为节点号 swing = find(ENodeTYPE==3); %平衡节点号 while max(abs(Det))>1e-5 %正式进入迭代求解循环 %% 计算电不平衡量Fe(全部节点2*ENodeNum个方程，需修正swing与PV的不平衡量,使其为0） Ui= sparse(U.*(cos(alpha)+1i*sin(alpha))); %将矩阵转化为稀疏矩阵存贮 S = Ui.*(conj(Y)*conj(Ui)); %% ————————每一次迭代都需要检查节点类型———————————— %每一步迭代后，应在此判断PV点无功是否越限，PV的无功限值存在Qs中， %如越限，在此处更改越限节点类型1并取节点无功为限值。巧合的是无功恰好放在了Qs中，第二步可不做。 NodesQ=imag(S); delPVPos=abs(NodesQ(PV))>abs(Qs(PV)); %越限，返回越限节点在PV中的位置，逻辑数组 ENodeTYPE(PV(delPVPos))=1; %越限PV变为PQ,数组A(逻辑数组B)返回A中与B=1对应位置的数 PV =find(ENodeTYPE==2); %PV节点可能越限，每一循环需更新节点类型 %% ——————计算电量不平衡并修正swing的PQ、PV的Q———————— detS=(Ps+1i*Qs)-S; detP=real(detS); detQ=imag(detS); detP(swing)=0; detQ(swing)=0; detQ(PV)=0; Fe=[detP;detQ]; %% Jee(需修正平衡与PV,对角为1，非对角为0) Slm=diag(S); Sij=diag(Ui)*conj(Y)*diag(conj(Ui)); HJ=1i*(Sij-Slm); NL=-Sij-Slm; Jee_H=real(HJ); % H Jee_J=imag(HJ); % J Jee_N=real(NL); % N Jee_L=imag(NL); % L %————————————————————————————————————— Jee =sparse([Jee_H,Jee_N;Jee_J,Jee_L]); %————————————————————————————————————— %修正PV节点和平衡节点的Jaccobian Jee(swing,:) = 0; %swing的P Jee(:,swing) = 0; Jee(swing,swing) = 1; Jee(ENodeNum+swing,:) = 0; %swing的Q Jee(:,ENodeNum+swing) = 0; Jee(ENodeNum+swing,ENodeNum+swing) = 1; Jee(ENodeNum+PV,:) = 0; %PV的Q Jee(:,ENodeNum+PV) = 0; Jee=Jee+sparse(ENodeNum+PV,ENodeNum+PV,ones(length(PV),1),2*ENodeNum,2*ENodeNum); %% ————解修正方程式并修正电压/相角——————— Det=-Jee\Fe; alpha= alpha+Det(1:ENodeNum); U = U+(Det((ENodeNum+1):2*ENodeNum)).*U; %% KK=KK+1; end %% 善后工作 Ui= sparse(U.*(cos(alpha)+1i*sin(alpha))); PQ_loss=sum(Ui.*(conj(Y)*conj(Ui))); %线路损耗 S = Ui.*(conj(Y)*conj(Ui)); Sswing=S(swing); %平衡节点功率 %——————以下为线路功率损耗，没有写变压器的损耗，懒得写了wanglei0704———— Sij = sparse(Linei(KT==1),Linej(KT==1),U(Linei(KT==1)).*U(Linei(KT==1)).*conj(1i*B0(KT==1)/2+1i*Qsh(KT==1))+... U(Linei(KT==1)).*(conj(Linei(KT==1))-conj(Linej(KT==1))).*conj(m_lineY(KT==1)), ENodeNum,ENodeNum)+... sparse(Linej(KT==1),Linei(KT==1),U(Linej(KT==1)).*U(Linej(KT==1)).*conj(1i*B0(KT==1)/2+1i*Qsh(KT==1))+... U(Linej(KT==1)).*(conj(Linej(KT==1))-conj(Linei(KT==1))).*conj(m_lineY(KT==1)), ENodeNum,ENodeNum); detSij = Sij+Sij'; %% toc