From cfa77f98e3953af5e95057c90b90f9961083de7f Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Sat, 11 Oct 2014 10:57:55 +0800 Subject: [PATCH] =?UTF-8?q?=E7=94=A8=E5=86=85=E7=82=B9=E6=B3=95=E6=94=B6?= =?UTF-8?q?=E6=95=9B=E4=BA=86=E3=80=82=E4=BD=86=E6=98=AF=E4=B9=8B=E5=89=8D?= =?UTF-8?q?=E7=94=A8=E5=BA=8F=E5=88=86=E9=87=8F=E7=AE=97=E6=BD=AE=E6=B5=81?= =?UTF-8?q?=E7=9A=84=E4=BB=A3=E7=A0=81=E8=A2=AB=E7=A0=B4=E5=9D=8F=E4=BA=86?= =?UTF-8?q?=EF=BC=8C=E5=87=86=E5=A4=87=E6=81=A2=E5=A4=8D=E4=B8=80=E4=B8=8B?= =?UTF-8?q?=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- FormG.m | 7 +-- FormH.m | 22 ++++++-- FormLw.m | 20 ++++--- FormLz.m | 22 ++++---- Modification.m | 24 ++++---- SolveIt.m | 16 +++--- func_deltF.m | 17 ++++-- func_deltG.m | 33 ++++++----- func_deltH.m | 23 ++++++-- func_deltdeltF.m | 24 ++++++-- run.m | 144 +++++++++++++++++++++++++---------------------- 11 files changed, 204 insertions(+), 148 deletions(-) diff --git a/FormG.m b/FormG.m index b264438..1c98270 100644 --- a/FormG.m +++ b/FormG.m @@ -1,7 +1,6 @@ -function Mat_G=FormG(Volt,PD,QD,Loadi) +function Mat_G=FormG(I1r,I1i) Mat_G=[ - PD(Loadi); - QD(Loadi); - Volt; + I1r; + I1i; ]; end \ No newline at end of file diff --git a/FormH.m b/FormH.m index 7ce3fba..9aff295 100644 --- a/FormH.m +++ b/FormH.m @@ -1,14 +1,24 @@ -function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle) +function Mat_H=FormH(fsY11,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance) %% %QDcos=textread('D:\Project\最小化潮流\最小潮流算例\仙海919PDQDglys.txt'); %QD(QD~=0)=PD(QD~=0)./tan(QDcos); %QD(QD_NON_ZERO_IND)=QD_NON_ZERO; %% %PD(Loadi)=QD(Loadi)./tan(acos(0.98)); -AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); -dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt; -dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; - -Mat_H=[dP;dQ;]; +% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); +% dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt; +% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt; +% +% Mat_H=[dP;dQ;]; +H=[ + real(fsY11),-imag(fsY11); + imag(fsY11),real(fsY11); +]; +busNum=length(fsY11); +% C=sparse(1:length(Loadi),Loadi,1,length(Loadi),busNum*2); +Mat_H=H*[V1r;V1i]; +Mat_H=Mat_H+sparse(Loadi,1,-I1r,busNum*2,1)+sparse(Loadi+busNum,1,-I1i,busNum*2,1); +Mat_H(Balance)=Mat_H(Balance)-BalI1r; +Mat_H(Balance+busNum)=Mat_H(Balance+busNum)-BalI1i; end \ No newline at end of file diff --git a/FormLw.m b/FormLw.m index 7c40469..aa7e137 100644 --- a/FormLw.m +++ b/FormLw.m @@ -1,10 +1,12 @@ -function Lw=FormLw(Mat_G,Init_U,Busnum,Loadi) -VoltU=sparse(1.2*ones(Busnum,1)); -PDU=sparse(2*ones(length(Loadi),1)); -QDU=sparse(2*ones(length(Loadi),1)); -t1=[PDU; - QDU; - VoltU]; -t2=Mat_G+Init_U-t1; -Lw=t2; +function Lw=FormLw(Loadi) +% VoltU=sparse(1.2*ones(Busnum,1)); +% PDU=sparse(2*ones(length(Loadi),1)); +% QDU=sparse(2*ones(length(Loadi),1)); +% t1=[PDU; +% QDU; +% VoltU]; +% t2=Mat_G+Init_U-t1; +% Lw=t2; + +Lw=0.1*sparse(ones(length(Loadi)*2,1)); end \ No newline at end of file diff --git a/FormLz.m b/FormLz.m index ccffe93..fa0ba8e 100644 --- a/FormLz.m +++ b/FormLz.m @@ -1,12 +1,12 @@ -function Lz=FormLz(Mat_G,Init_L,Busnum,Loadi) -VoltL=sparse(0.8*ones(Busnum,1)); -PDL=sparse(-2*ones(length(Loadi),1)); -QDL=sparse(-2*ones(length(Loadi),1)); -t1=[PDL; - QDL; - VoltL - ]; -t2=Mat_G-Init_L-t1; -Lz=t2; - +function Lz=FormLz(Loadi) +% VoltL=sparse(0.8*ones(Busnum,1)); +% PDL=sparse(-2*ones(length(Loadi),1)); +% QDL=sparse(-2*ones(length(Loadi),1)); +% t1=[PDL; +% QDL; +% VoltL +% ]; +% t2=Mat_G-Init_L-t1; +% Lz=t2; +Lz=-0.1*sparse(ones(length(Loadi)*2,1)); end \ No newline at end of file diff --git a/Modification.m b/Modification.m index d55962c..e3d4c40 100644 --- a/Modification.m +++ b/Modification.m @@ -1,4 +1,5 @@ -function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi) +function [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY, ... + V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi) AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); Init_Z=Init_Z+AlphaD*deltZ; @@ -6,13 +7,16 @@ Init_L=Init_L+AlphaP*deltL; Init_W=Init_W+AlphaD*deltW; Init_U=Init_U+AlphaP*deltU; Init_Y=Init_Y+AlphaD*deltY; -t=deltX(1:size(Loadi,1)*2); -PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi)); -QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2); -t=deltX(size(Loadi,1)*2+1:ContrlCount); -t(Busnum+Balance)=0; -balVolt=Volt(Balance); -Volt=Volt+AlphaP*t(1:Busnum); -Volt(Balance)=balVolt; -UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); +t=deltX(1:busNum*2); +V1r=V1r+AlphaP*t(1:busNum); +V1r(Balance)=1; +V1i=V1i+AlphaP*t(busNum+1:end); +V1i(Balance)=0; +t=deltX(busNum*2+1:ContrlCount); +I1r=I1r+AlphaP*t(1:length(Loadi)); +I1i=I1i+AlphaP*t(length(Loadi)+1:end); +%balVolt=1;%Volt(Balance); +%Volt=Volt+AlphaP*t(1:Busnum); +%Volt(Balance)=balVolt; +%UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum); end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m index c41122f..4e22f8d 100644 --- a/SolveIt.m +++ b/SolveIt.m @@ -8,14 +8,14 @@ aa=[ deltH',zeros(length(Init_Y)); ]; yy=[LxComa;-Ly]; -%% 平衡节点电压不变 -t=size(Loadi,1)*2; -aa(t+Balance,:)=0; -aa(:,t+Balance)=0; -aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); -deltG(t+Balance,:)=0; -%% -t=size(Loadi,1)*2+Busnum*1; +%% 平衡节点电压实部 + +aa(Balance,:)=0; +aa(:,Balance)=0; +aa=aa+sparse(Balance,Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); +deltG(Balance,:)=0; +%% 平衡节点电压虚部 +t=Busnum*1; aa(t+Balance,:)=0; aa(:,t+Balance)=0; aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); diff --git a/func_deltF.m b/func_deltF.m index 7bfb3c6..980e4c7 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -1,8 +1,15 @@ -function deltF=func_deltF(SEPD,ContrlCount) -t=sparse(ones(length(SEPD),1)); +function deltF=func_deltF(measurement,busNum,fsY11,Loadi,V1r,V1i,I1r,I1i) +% H=[ +% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); +% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); +% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); +% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); +% ]; +% X=[V1r;V1i;I1r;I1i]; +% deltF=-H'*(measurement-H*X); deltF=[ - t; - sparse(ContrlCount-length(t),1); + zeros(busNum*2,1); + -2*( [real(measurement);imag(measurement)]-[I1r;I1i]); ]; - +% deltF=0; end \ No newline at end of file diff --git a/func_deltG.m b/func_deltG.m index 0522de4..2657442 100644 --- a/func_deltG.m +++ b/func_deltG.m @@ -1,22 +1,21 @@ function deltG=func_deltG(busNum,Loadi) -dgV_dPD=sparse(length(Loadi),busNum); -dgV_dQD=sparse(length(Loadi),busNum); -dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum); - sparse(busNum,busNum); - ]; -% -dgPD_dPD=sparse(eye(length(Loadi))); -dgPD_dQD=sparse(length(Loadi),length(Loadi)); -dgPD_dx=sparse(busNum*2,length(Loadi)); -% -dgQD_dPD=sparse(length(Loadi),length(Loadi)); -dgQD_dQD=sparse(eye(length(Loadi))); -dgQD_dx=sparse(busNum*2,length(Loadi)); -%% +% dgV_dPD=sparse(length(Loadi),busNum); +% dgV_dQD=sparse(length(Loadi),busNum); +% dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum); +% sparse(busNum,busNum); +% ]; +% % +% dgPD_dPD=sparse(eye(length(Loadi))); +% dgPD_dQD=sparse(length(Loadi),length(Loadi)); +% dgPD_dx=sparse(busNum*2,length(Loadi)); +% % +% dgQD_dPD=sparse(length(Loadi),length(Loadi)); +% dgQD_dQD=sparse(eye(length(Loadi))); +% dgQD_dx=sparse(busNum*2,length(Loadi)); +% %% deltG=[ - dgPD_dPD,dgQD_dPD,dgV_dPD; - dgPD_dQD,dgQD_dQD,dgV_dQD; - dgPD_dx,dgQD_dx,dgV_dx; + zeros(busNum*2,length(Loadi)*2); + eye(length(Loadi)*2); ]; end \ No newline at end of file diff --git a/func_deltH.m b/func_deltH.m index fa65c0a..c69b752 100644 --- a/func_deltH.m +++ b/func_deltH.m @@ -1,7 +1,18 @@ -function deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Angle,Loadi) -dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)]; -dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)]; -dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 -deltH=[dH_dPD;dH_dQD;dH_dx']; - +function deltH=func_deltH(busNum,fsY11,Loadi,Balance) +% dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)]; +% dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)]; +% dH_dx = jacobian(Busnum,Volt,Y,Angle,UAngel,r,c); %形成雅克比矩阵 +% deltH=[dH_dPD;dH_dQD;dH_dx']; +H=[ + real(fsY11),-imag(fsY11); + imag(fsY11),real(fsY11); + +]; +deltH=[ + H'; + sparse(1:length(Loadi),Loadi,-1,length(Loadi),busNum*2); + sparse(1:length(Loadi),Loadi+busNum,-1,length(Loadi),busNum*2); + %sparse(1:length(Balance),Balance,-1,1,busNum*2); + %sparse(1:length(Balance),Balance+busNum,-1,1,busNum*2); +]; end \ No newline at end of file diff --git a/func_deltdeltF.m b/func_deltdeltF.m index 309e1e9..2003458 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -1,6 +1,22 @@ -function deltdeltF=func_deltdeltF(SEPD,ContrlCount) +function deltdeltF=func_deltdeltF(busNum,fsY11,Loadi) +% H=[ +% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); +% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); +% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); +% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); +% ]; + +% H=[ +% real(fsY11),-imag(fsY11),zeros(busNum,length(Loadi)*2); +% imag(fsY11),real(fsY11),zeros(busNum,length(Loadi)*2); +% zeros(length(Loadi),busNum*2), eye(length(Loadi)),zeros(length(Loadi)); +% zeros(length(Loadi),busNum*2), zeros(length(Loadi)),eye(length(Loadi)); +% ]; + +% deltdeltF=H'*H; deltdeltF=[ - 0*zeros(ContrlCount); - ]; -deltdeltF=0; + zeros(busNum*2,busNum*2+length(Loadi)*2); + zeros(length(Loadi)*2,busNum*2),2*eye(length(Loadi)*2); +]; +% deltdeltF=0; end \ No newline at end of file diff --git a/run.m b/run.m index c2bf845..4eba75a 100644 --- a/run.m +++ b/run.m @@ -57,6 +57,8 @@ tic VoltpA=sparse(ones(busNum,1)); VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi); + + while(k<=kmax && maxD> EPS) k=k+1; %把补偿电容看作负荷 @@ -81,51 +83,17 @@ while(k<=kmax && maxD> EPS) 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; +% fsY11*V1; %形成负荷序电流的测量值 mIf0=If0; mIf1=If1; mIf1(3)=-mIf1(2); -% 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，否则是奇异的 +% fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,bus +% Num);%这里要置0，置1，否则是奇异的 %%做最小二乘法 - %先做正序的 - Z=[ - -real(mIf1(Loadi));%这里要加-号，因为用的Z是注入电流 - -imag(mIf1(Loadi));%这里要加-号，因为用的Z是注入电流 - %加电压量测 - [1 1 1]'; - ]; - H=fsY11(Loadi,:); - H=[ - -real(H),imag(H),; - -imag(H),-real(H),; - eye(3),eye(3); - ]; - % 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;])); @@ -224,28 +192,65 @@ fprintf(' % VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end)); %% 用牛顿法求解end %% 开始进入状态估计 +%准备量测量 +iterPhaseASpotLoadP=phaseASpotLoadP; +iterPhaseBSpotLoadP=phaseBSpotLoadP; +iterPhaseCSpotLoadP=phaseCSpotLoadP; +iterPhaseASpotLoadQ=phaseASpotLoadQ; +iterPhaseBSpotLoadQ=phaseBSpotLoadQ; +iterPhaseCSpotLoadQ=phaseCSpotLoadQ; +%全部转换为负荷电流 +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,:)'); +%试着算一下正序电流 +% fsY11*V1; +%形成负荷序电流的测量值 +mIf0=If0; +mIf1=-If1; +% mIf1(3)=-mIf1(2); +mIf2=If2; +% fsY11(:,Balance)=0; +% fsY11(Balance,:)=0; +% fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); +% mIf1(3)=1; +%平衡节点电流 +BalI1r=real(-sum(mIf1)); +BalI1i=imag(-sum(mIf1)); +inv(fsY11)*(mIf1); +V1r=1*ones(busNum,1); +V1i=1*ones(busNum,1); +I1r=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 +I1i=0.01*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反 +measurement=-mIf1(2); clear PD QD PG QG; %状态量 -SEVoltpA=sparse(ones(busNum,1)); -SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); -SEVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi); -SEphaseASpotLoadP=zeros(length(phaseASpotLoadP),1); -SEphaseBSpotLoadP=zeros(length(phaseBSpotLoadP),1); -SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1); -SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1); -SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1); -SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1); +% SEVoltpA=sparse(ones(busNum,1)); +% SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); +% SEVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi); +% SEphaseASpotLoadP=zeros(length(phaseASpotLoadP),1); +% SEphaseBSpotLoadP=zeros(length(phaseBSpotLoadP),1); +% SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1); +% SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1); +% SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1); +% SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1); % -SEVmf1=sparse(ones(busNum,1)); -SEVaf1=sparse(zeros(busNum,1)); -SEPD=sparse(zeros(busNum,1)); -SEQD=sparse(zeros(busNum,1)); +% SEVmf1=sparse(ones(busNum,1)); +% SEVaf1=sparse(zeros(busNum,1)); +% SEPD=sparse(zeros(busNum,1)); +% SEQD=sparse(zeros(busNum,1)); KK=0; plotGap=zeros(1,60); %初始化 %状态量为 SEPD SEQD SEVmf1 SEVaf1 -RestraintCount=length(SEVmf1)+length(Loadi)*2; -ContrlCount=length(SEVmf1)*2+length(Loadi)*2; +RestraintCount=length(Loadi)*2; +ContrlCount=busNum*2+length(Loadi)*2; CenterA=0.1; Init_Z=sparse(ones(RestraintCount,1)); Init_W=sparse(-1*ones(RestraintCount,1)); @@ -253,30 +258,31 @@ Init_L=1*sparse(ones(RestraintCount,1)); Init_U=1*sparse(ones(RestraintCount,1)); Init_Y=sparse(2*busNum,1);%等式约束乘子 Gap=(Init_L'*Init_Z-Init_U'*Init_W); -PG=sparse(busNum,1); -PG(Balance)=0.1105; -QG=sparse(busNum,1); -QG(Balance)=0.0984; -SEPD(2)=0.1105; -SEQD(2)=0.0984; +% PG=sparse(busNum,1); +% PG(Balance)=0.1105; +% QG=sparse(busNum,1); +% QG(Balance)=0.0984; +% SEPD(2)=0.1105; +% SEQD(2)=0.0984; while Gap>1e-5 && KK<20 KK=KK+1; Init_u=Gap/2/RestraintCount*CenterA; - deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi); + deltH=func_deltH(busNum,fsY1,Loadi,Balance); deltG=func_deltG(busNum,Loadi); L_1Z=diag(Init_Z./Init_L); U_1W=diag(Init_W./Init_U); - deltdeltF=func_deltdeltF(SEPD,ContrlCount); - ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); + deltdeltF=func_deltdeltF(busNum,fsY1,Loadi); +% ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); + ddh=0; ddg=func_ddg(); - deltF=func_deltF(SEPD,ContrlCount); + deltF=func_deltF(measurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i); Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); - Mat_G=FormG(SEVmf1,SEPD,SEQD,Loadi); - Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang); + Mat_G=FormG(I1r,I1i); + Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance); Ly=Mat_H; - Lz=FormLz(Mat_G,Init_L,busNum,Loadi); - Lw=FormLw(Mat_G,Init_U,busNum,Loadi); + Lz=FormLz(Loadi); + Lw=FormLw(Loadi); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 @@ -284,6 +290,8 @@ while Gap>1e-5 && KK<20 fprintf('迭代次数 %d Gap %f\n',KK,plotGap(KK)); XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi); [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum); - [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,SEVmf1,SEVaf1,SEPD,SEQD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,SEPD,SEQD,Loadi); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi); Gap=(Init_L'*Init_Z-Init_U'*Init_W); -end \ No newline at end of file +end +%检查目标函数 +f=sum([real(measurement);imag(measurement)]-[-I1r;-I1i]); \ No newline at end of file