From 1ffa3c0bf3873e8f6710c7646b384f9ca2d17620 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Wed, 8 Oct 2014 16:21:02 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=9D=E6=AD=A5=E5=AE=8C=E6=88=90=E4=BA=86?= =?UTF-8?q?=E5=8F=AA=E7=94=A8=E6=AD=A3=E5=BA=8F=E7=9A=84?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- run.m | 69 +++++++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 58 insertions(+), 11 deletions(-) diff --git a/run.m b/run.m index e88d539..6feda3e 100644 --- a/run.m +++ b/run.m @@ -1,5 +1,5 @@ %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。 -% +% clc clear lineZ=readLineZ('feeder13\lineParameter.txt'); @@ -41,9 +41,11 @@ If2=sparse(busNum,1); Vf0=sparse(busNum,1); If0=sparse(busNum,1); %准备序矩阵 +%平衡节点置0置1 fsY2(Balance,:)=0; fsY2(:,Balance)=0; fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); +%平衡节点置0置1 fsY0(Balance,:)=0; fsY0(:,Balance)=0; fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); @@ -69,6 +71,57 @@ while(k<=kmax && maxD> EPS) iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA); iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB); iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC); + %全部转换为负荷电流 + CurpA=conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA); + CurpB=conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB); + CurpC=conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC); + %转换为序电流 + f012=Tp2f*conj([CurpA';CurpB';CurpC']); + %把三序电流分离出来 + If0=conj(f012(1,:)'); + If1=conj(f012(2,:)'); + If2=conj(f012(3,:)'); + %得到三序电压 + V012=Tp2f*conj([VoltpA';VoltpB';VoltpC']); + %分离出三序电压 + V0=conj(V012(1,:)'); + V1=conj(V012(2,:)'); + V2=conj(V012(3,:)'); + %试着算一下正序电流 + fsY11*V1; + %形成负荷序电流的测量值 + mIf0=If0; + mIf1=If1; + mIf1(1)=0.02; + mIf1(3)=0.03; + Loadi=[ 1 2 3]; + mIf2=If2; + %计算 + fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);%这里要置0，置1，否则是奇异的 + %%做最小二乘法 + %先做正序的 + Z=[%这里要加-号，因为用的Z是注入电流 + -real(mIf1(Loadi)); + -imag(mIf1(Loadi)); + ]; + H=fsY11(Loadi,:); + H=[ + -real(H),imag(H); + -imag(H),-real(H); + ]; + % J=fsY11(Loadi,:); + % J=[ + % conj(-real(J))',-conj(imag(J))'; + % conj(imag(J))',-conj(real(J))'; + % ]; + J=conj(H'); + %J*Z+J*H*X=0 ->J*H*X=-J*Z + X=-inv(J*H)*J*Z; + % X=inv(conj(H)'*H)*conj(H)'*Z; + Xr=X(1:length(X)/2); + Xi=X(length(X)/2+1:end); + Xri=Xr+1j*Xi; + fsY11*Xri; [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ... PQi,PG,QG,QGi,iterPD,iterQD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%不平衡量 maxD=max(abs([dP;dQ;])); @@ -79,15 +132,9 @@ while(k<=kmax && maxD> EPS) %转换为三相电压 VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%用Tp2f\ 代替Tf2p* VoltpA=conj(VoltpABC(1,:)'); - CurpA=-conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA); - VoltpB=conj(VoltpABC(2,:)'); - CurpB=-conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB); - VoltpC=conj(VoltpABC(3,:)'); - CurpC=-conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC); - f012=Tp2f*conj([CurpA';CurpB';CurpC']); - If0=conj(f012(1,:)'); - If1=conj(f012(2,:)'); - If2=conj(f012(3,:)'); + + + If0(Balance)=0; If2(Balance)=0; %Vf0=fsY0\If0; @@ -116,8 +163,8 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ... phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ... phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ ); fprintf('最大不平衡量为%f\n\n',full(max(abs(ub)))) -fprintf('开始牛顿法迭代\n'); %% 用牛顿法求解begin +% fprintf('开始牛顿法迭代\n'); % [r,c,GB]=find(phaseABCY); % Y=abs(phaseABCY); % Yangle=angle(GB);