用内点法收敛了。但是之前用序分量算潮流的代码被破坏了，准备恢复一下。

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
2014-10-11 10:57:55 +08:00 · 2014-10-11 10:57:55 +08:00 · cfa77f98e3
parent 4daca01ca0
commit cfa77f98e3
11 changed files with 204 additions and 148 deletions
--- 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
--- 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
--- 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
--- 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
--- 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
--- 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);
--- 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
--- 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
--- 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
--- 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
--- 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
+end
+%检查目标函数
+f=sum([real(measurement);imag(measurement)]-[-I1r;-I1i]);