distributionnetwork-power2current/run.m

function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,isConverged]=run()
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȰѸ<EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><EFBFBD>֪<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%
close all
clc
clear
lineZ=readLineZ('feeder13\lineParameter.txt');
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
    cap]=dataRead(lineZ,'feeder13\data1.txt');
a=exp(1j*2*pi/3);
Tp2f=1/3*[1 1 1;
    1 a a^2;
    1 a^2 a];
Tp2f=sparse(Tp2f);
Tf2p=inv(Tp2f);
fsY1amp=abs(fsY1);
[r,c,fsY1ang]=find(fsY1);
fsY1ang=angle(fsY1ang);
Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP;
Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ;
busNum=length(phaseASpotLoadP);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Vmf1=sparse(ones(busNum,1));
Vaf1=sparse(zeros(busNum,1));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
PQi=nodeNum;
PG=sparse(busNum,1);
QG=sparse(busNum,1);
QGi=[Balance];
PD=Pabc/3;
QD=Qabc/3;
Loadi=find(PD~=0);
maxD=100000;% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
EPS=1e-5;
k=0;
kmax=20;
fsY11=fsY1;
fsY00=fsY0;
fsY22=fsY2;
Vf2=sparse(busNum,1);
If2=sparse(busNum,1);
Vf0=sparse(busNum,1);
If0=sparse(busNum,1);
%׼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
fsY2(Balance,:)=0;
fsY2(:,Balance)=0;
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
fsY0(Balance,:)=0;
fsY0(:,Balance)=0;
fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%%LU<EFBFBD>ֽ<EFBFBD>
[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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;
    %<EFBFBD>Ѳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    SA=VoltpA.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,1),busNum,1));
    SB=VoltpB.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
    SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
    %<EFBFBD>Ȳ<EFBFBD>Ҫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    SA=0;
    SB=0;
    SC=0;
    iterPD=PD+real(SA+SB+SC)/3;
    iterQD=QD+imag(SA+SB+SC)/3;
    iterPhaseASpotLoadP=phaseASpotLoadP+real(SA);
    iterPhaseBSpotLoadP=phaseBSpotLoadP+real(SB);
    iterPhaseCSpotLoadP=phaseCSpotLoadP+real(SC);
    iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA);
    iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB);
    iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC);
    %%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD>˷<EFBFBD>
    [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);%<EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    maxD=max(abs([dP;dQ;]));
    jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%<EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
    [dV, dVangle]=Solv(busNum,jaco,dP,dQ);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [Vmf1, Vaf1]=Modify(Vmf1,Vaf1,dV,dVangle,1);
    fprintf('<EFBFBD><EFBFBD> %d <EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ %f\n',k,full(maxD));
    %ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
    VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%<EFBFBD><EFBFBD>Tp2f\ <EFBFBD><EFBFBD><EFBFBD><EFBFBD>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;
    Vf0=fsY0Q*(fsY0U\(fsY0L\(fsY0P*(fsY0R\If0))));
    %Vf2=fsY2\If2;
    Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
    fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>%f\n',toc);
end
FortiscueToc=toc;
fprintf('Fortiscue<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> %f\n',FortiscueToc);
Vf1=Vmf1.*exp(1j*Vaf1);
%%
(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݵĹ<EFBFBD><EFBFBD><EFBFBD>
conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
IpABC=Tf2p*conj([If0';If1';If2']);
%ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
VoltpABC=Tf2p*conj([ Vf0'; Vf1'; Vf2']);
disp(['                  A                        B                          C'])
full(abs(VoltpABC'))
fprintf('<EFBFBD>ڵ<EFBFBD><EFBFBD>Ŷ<EFBFBD>Ӧ\n');
disp([setIJ,nodeNum ])
%%<EFBFBD><EFBFBD><EFBFBD>鷴<EFBFBD>ƻ<EFBFBD>ȥ<EFBFBD>Ĺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
    phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ%f\n\n',full(max(abs(ub))))
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>״̬<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% clear PD QD PG QG;
%׼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rIf0=If0;
rIf1=If1;
rIf2=If2;
sigma=0.01;
iterPhaseASpotLoadP=phaseASpotLoadP;
iterPhaseBSpotLoadP=phaseBSpotLoadP;
iterPhaseCSpotLoadP=phaseCSpotLoadP;
iterPhaseASpotLoadQ=phaseASpotLoadQ;
iterPhaseBSpotLoadQ=phaseBSpotLoadQ;
iterPhaseCSpotLoadQ=phaseCSpotLoadQ;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
mphaseASpotLoadP=phaseASpotLoadP.*(1+normrnd(0,sigma,length(phaseASpotLoadP),1));
mphaseBSpotLoadP=phaseBSpotLoadP.*(1+normrnd(0,sigma,length(phaseBSpotLoadP),1));
mphaseCSpotLoadP=phaseCSpotLoadP.*(1+normrnd(0,sigma,length(phaseCSpotLoadP),1));
mphaseASpotLoadQ=phaseASpotLoadQ.*(1+normrnd(0,sigma,length(phaseASpotLoadQ),1));
mphaseBSpotLoadQ=phaseBSpotLoadQ.*(1+normrnd(0,sigma,length(phaseBSpotLoadQ),1));
mphaseCSpotLoadQ=phaseCSpotLoadQ.*(1+normrnd(0,sigma,length(phaseCSpotLoadQ),1));
%<EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1));
mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1));
mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1));
%
mVoltpA=abs(VoltpA).*(1+normrnd(0,sigma,length(VoltpA),1));
mVoltpB=abs(VoltpB).*(1+normrnd(0,sigma,length(VoltpB),1)).*exp(1j*-120/180*pi);
mVoltpC=abs(VoltpC).*(1+normrnd(0,sigma,length(VoltpC),1)).*exp(1j*+120/180*pi);
%
% mVoltpA=sparse(ones(busNum,1));
% mVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
% mVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
%ת<EFBFBD><EFBFBD><EFBFBD>õ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
fV012=Tp2f*conj([mVoltpA';mVoltpB';mVoltpC']);
mfV0=conj(fV012(1,:)');
mfV1=conj(fV012(2,:)');
mfV2=conj(fV012(3,:)');
% 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);
%ȫ<EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD>ɵ<EFBFBD><EFBFBD><EFBFBD>
mCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA);
mCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB);
mCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC);
%ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
mf012=Tp2f*conj([mCurpA';mCurpB';mCurpC']);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
mIf0=conj(mf012(1,:)');
mIf1=conj(mf012(2,:)');%<EFBFBD><EFBFBD><EFBFBD>ɵ<EFBFBD><EFBFBD><EFBFBD>
mIf2=conj(mf012(3,:)');
mIf0=-mIf0;
mIf1=-mIf1;%mIf1<EFBFBD><EFBFBD>ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD>൱<EFBFBD>ڷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
mIf2=-mIf2;
%% <EFBFBD>趨<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
BalI1r=real(-sum(rIf1));
BalI1i=imag(-sum(rIf1));
%<EFBFBD><EFBFBD>ѹ
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% mfsY11=fsY11;
% mfsY11(:,Balance)=0;
% mfsY11(Balance,:)=0;
% mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV1=inv(mfsY11)*(mIf1)+1;
% sigmaV1=normrnd(0,sigma,length(Loadi),1);
% V1measurement=rV1(Loadi).*(1+sigmaV1);
V1measurement=mfV1(Loadi);
% wV1r=abs(real( rV1(Loadi).*sigmaV1 ));
% wV1i=abs(imag( rV1(Loadi).*sigmaV1 ));
% sigmaI1=normrnd(0,sigma,length(Loadi),1);
% I1measurement=mIf1(Loadi).*(1+sigmaI1);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD>ǵ<EFBFBD>Ч<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
I1measurement=mIf1(Loadi);
% wI1r=abs( real(mIf1(Loadi).*sigmaI1) );
% wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵڶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ

wI1r=abs(real(I1measurement)).*sigma;
wI1i=abs(imag(I1measurement)).*sigma;
wV1r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV1i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV1r=wV1r(Loadi);
wV1i=wV1i(Loadi);
% [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
% fprintf('Ŀ<EFBFBD><EFBFBD>ֵ %f\n',full(f));
%% <EFBFBD>㸺<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
BalI2r=real(-sum(rIf2));
BalI2i=imag(-sum(rIf2));
%<EFBFBD><EFBFBD>ѹ
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% mfsY22=fsY22;
% mfsY22(:,Balance)=0;
% mfsY22(Balance,:)=0;
% mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV2=inv(mfsY22)*(mIf2);
% sigmaV2=normrnd(0,sigma,length(Loadi),1);
% V2measurement=rV2(Loadi).*(1+sigmaV2);
V2measurement=mfV2(Loadi);
% wV2r=abs(real( rV2(Loadi).*sigmaV2 ));
% wV2i=abs(imag( rV2(Loadi).*sigmaV2 ));
% sigmaI2=normrnd(0,sigma,length(Loadi),1);
% I2measurement=mIf2(Loadi).*(1+sigmaI2);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD>ǵ<EFBFBD>Ч<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
I2measurement=mIf2(Loadi);
% wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
% wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵڶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ
wI2r=abs(real(I2measurement)).*sigma;
wI2i=abs(imag(I2measurement)).*sigma;
wV2r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV2i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV2r=wV2r(Loadi);
wV2i=wV2i(Loadi);

% [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
% f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
% fprintf('Ŀ<EFBFBD><EFBFBD>ֵ %f\n',full(f));
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
BalI0r=real(-sum(rIf0));
BalI0i=imag(-sum(rIf0));
%<EFBFBD><EFBFBD>ѹ
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% mfsY00=fsY00;
% mfsY00(:,Balance)=0;
% mfsY00(Balance,:)=0;
% mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV0=inv(mfsY00)*(mIf0);
% sigmaV0=normrnd(0,sigma,length(Loadi),1);
% V0measurement=rV0(Loadi).*(1+sigmaV0);
V0measurement=mfV0(Loadi);
% wV0r=abs(real( rV0(Loadi).*sigmaV0 ));
% wV0i=abs(imag( rV0(Loadi).*sigmaV0 ));
% sigmaI0=normrnd(0,sigma,length(Loadi),1);
% I0measurement=mIf0(Loadi).*(1+sigmaI0);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD>ǵ<EFBFBD>Ч<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
I0measurement=mIf0(Loadi);
% wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
% wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵڶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ

wI0r=abs(real(I0measurement)).*sigma;
wI0i=abs(imag(I0measurement)).*sigma;
wV0r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV0i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV0r=wV0r(Loadi);
wV0i=wV0i(Loadi);
% matlabpool local 3

%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD>µ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% %A<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% dIrdP_Ar=real(mVoltpA)./(real(mVoltpA).^2+imag(mVoltpA).^2);
% dIrdQ_Ar=imag(mVoltpA)./(real(mVoltpA).^2+imag(mVoltpA).^2);
% % dIrdP_Ar=dIrdP_Ar(Loadi);
% % dIrdQ_Ar=dIrdQ_Ar(Loadi);
% dI_pAr=dIrdP_Ar.*mphaseASpotLoadP*0.2+mphaseASpotLoadQ.*dIrdQ_Ar*0.2;
% dIidP_Ai=-imag(mVoltpA)./(real(mVoltpA).^2+imag(mVoltpA).^2);
% dIidQ_Ai=-real(mVoltpA)./(real(mVoltpA).^2+imag(mVoltpA).^2);
% % dIidP_Ai=dIidP_Ai(Loadi);
% % dIidQ_Ai=dIidQ_Ai(Loadi);
% dI_pAi=dIidP_Ai.*mphaseASpotLoadP*0.2+dIidQ_Ai.*mphaseASpotLoadQ*0.2;
% 
% %B<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% dIrdP_Br=real(mVoltpB)./(real(mVoltpB).^2+imag(mVoltpB).^2);
% dIrdQ_Br=imag(mVoltpB)./(real(mVoltpB).^2+imag(mVoltpB).^2);
% % dIrdP_Br=dIrdP_Br(Loadi);
% % dIrdQ_Br=dIrdQ_Br(Loadi);
% dI_pBr=dIrdP_Br.*mphaseBSpotLoadP*0.2+dIrdQ_Br.*mphaseBSpotLoadQ*0.2;
% dIidP_Bi=-imag(mVoltpB)./(real(mVoltpB).^2+imag(mVoltpB).^2);
% dIidQ_Bi=-real(mVoltpB)./(real(mVoltpB).^2+imag(mVoltpB).^2);
% % dIidP_Bi=dIidP_Bi(Loadi);
% % dIidQ_Bi=dIidQ_Bi(Loadi);
% dI_pBi=dIidP_Bi.*mphaseBSpotLoadP*0.2+dIidQ_Bi.*mphaseBSpotLoadQ*0.2;
% 
% %C<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% dIrdP_Cr=real(mVoltpC)./(real(mVoltpC).^2+imag(mVoltpC).^2);
% dIrdQ_Cr=imag(mVoltpC)./(real(mVoltpC).^2+imag(mVoltpC).^2);
% % dIrdP_Cr=dIrdP_Cr(Loadi);
% % dIrdQ_Cr=dIrdQ_Cr(Loadi);
% dI_pCr=dIrdP_Cr.*mphaseCSpotLoadP*0.2+dIrdQ_Cr.*mphaseBSpotLoadQ*0.2;
% dIidP_Ci=-imag(mVoltpC)./(real(mVoltpC).^2+imag(mVoltpC).^2);
% dIidQ_Ci=-real(mVoltpC)./(real(mVoltpC).^2+imag(mVoltpC).^2);
% % dIidP_Ci=dIidP_Ci(Loadi);
% % dIidQ_Ci=dIidQ_Ci(Loadi);
% dI_pCi=dIidP_Ci.*mphaseCSpotLoadP*0.2+dIidQ_Ci.*mphaseCSpotLoadQ*0.2;
% % <EFBFBD>ϳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% dI_F=Tp2f*[dI_pAr'+1j*dI_pAi';
%     dI_pBr'+1j*dI_pBi';
%     dI_pCr'+1j*dI_pCi';
%     ];
% dI_F=conj(dI_F)';
% dI_F=dI_F(Loadi,:);
% <EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
deltmCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA)*0.2;
deltmCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB)*0.2;
deltmCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC)*0.2;
%ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
deltF012=Tp2f*conj([deltmCurpA';deltmCurpB';deltmCurpC']);
dI_F=conj(deltF012)';
dI_F=dI_F(Loadi,:);
tic
for II=1:3
    
    if II==1
        fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
        tic
        [ V1r,V1i,I1r,I1i,isConverged1 ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1,dI_F,II );
        toc
    end
    if II==2
        fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
        tic
        [ V2r,V2i,I2r,I2i,isConverged2 ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0,dI_F,II );
        toc
    end
    if II==3
        fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
        tic
        [ V0r,V0i,I0r,I0i,isConverged0 ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0,dI_F,II );
        toc
    end
    
end
toc
% matlabpool close

% f=sum(([real(I0measurement);imag(I0measurement)]-[I0r;I0i]).^2)+sum((real(rV0)-V0r).^2)+sum((imag(rV0)-V0i).^2);
% fprintf('Ŀ<EFBFBD><EFBFBD>ֵ %f\n',full(f));
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
SEVf0=V0r+1j*V0i;
SEVf1=V1r+1j*V1i;
SEVf2=V2r+1j*V2i;
SEVoltpABC=Tf2p*conj([ SEVf0'; SEVf1'; SEVf2']);
% SEVoltpABC2=Tf2p*conj([ rV0'; rV1'; rV2']);
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
SEIf0=I0r+1j*I0i;
SEIf1=I1r+1j*I1i;
SEIf2=I2r+1j*I2i;
SEIpABC=full(Tf2p*conj([SEIf0';SEIf1(1:3)';SEIf2']));
%<EFBFBD><EFBFBD>һ<EFBFBD>²<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
VError=(full(abs(VoltpABC))-abs(SEVoltpABC))./abs(VoltpABC)*100;
VError=reshape(VError,size(VError,1)*size(VError,2),1);
% barPlot( VError, 10,['<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%'],['<EFBFBD>ֲ<EFBFBD><EFBFBD>ܶ<EFBFBD>'],['<EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD>ֵ']);
% <EFBFBD><EFBFBD><EFBFBD>ฺ<EFBFBD><EFBFBD>
rThreeLoad=[ phaseASpotLoadP'+1j*phaseASpotLoadQ';
    phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
    phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
    ];
% rThreeLoad=rThreeLoad(:,setxor(1:size(SEVoltpABC,2),Balance));
SEThreeLoad=SEVoltpABC(:,setxor(1:size(SEVoltpABC,2),Balance)).*conj(-SEIpABC);
% phaseLoadPError=real(rThreeLoad-SEThreeLoad)./real(rThreeLoad)*100;
% phaseLoadQError=imag(rThreeLoad-SEThreeLoad)./imag(rThreeLoad)*100;
% phaseLoadPError=reshape(phaseLoadPError,size(phaseLoadPError,1)*size(phaseLoadPError,2),1);
% phaseLoadQError=reshape(phaseLoadQError,size(phaseLoadQError,1)*size(phaseLoadQError,2),1);
% figure()
% barPlot( phaseLoadPError, 10,['<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%'],['<EFBFBD>ֲ<EFBFBD><EFBFBD>ܶ<EFBFBD>'],['<EFBFBD>й<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>']);
% figure()
% barPlot( phaseLoadQError, 10,'<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%','<EFBFBD>ֲ<EFBFBD><EFBFBD>ܶ<EFBFBD>','<EFBFBD>޹<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>');
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
rVoltABCV=conj([VoltpA';VoltpC';VoltpC']);
mVoltABCV=conj([mVoltpA';mVoltpB';mVoltpC']);
rPD3P=real([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
    phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
    phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
    ]);
rPD3P=rPD3P(:,Loadi);
mPD3P=real([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
    mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
    mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
    ]);
mPD3P=mPD3P(:,Loadi);
PD3P=real(SEThreeLoad);
rQD3P=imag([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
    phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
    phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
    ]);
rQD3P=rQD3P(:,Loadi);
mQD3P=imag([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
    mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
    mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
    ]);
mQD3P=mQD3P(:,Loadi);
QD3P=imag(SEThreeLoad);
JMeasurement=sum(sum((( abs(mVoltABCV)-abs(SEVoltpABC) )./ abs(mVoltABCV)./sigma).^2))+sum(sum(((mPD3P-PD3P)./mPD3P./sigma).^2))+sum(sum(((mQD3P-QD3P)./mQD3P./sigma).^2));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͳ<EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Busnum=busNum;
mCount=Busnum*3+length(Loadi)*3*2;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
AME_Volt=sum(sum(abs( abs(mVoltABCV)-abs(SEVoltpABC))));
AME_VAngle=sum(sum(abs( angle(mVoltABCV)-angle(SEVoltpABC))));
AME_PD=sum(sum(abs(real(SEThreeLoad-rThreeLoad(:,Loadi)))));
AME_QD=sum(sum(abs(imag(SEThreeLoad-rThreeLoad(:,Loadi)))));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ
isConverged=isConverged1*isConverged2*isConverged0;
end