distributionnetwork-power2current/IPMLoop.m

function [ V1r,V1i,I1r,I1i,isConverged,calcuTime ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref,dI_F,flag,guessIf1,rIf1,noLoadi )
%<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>
V1r=Vref*ones(busNum,1);
V1i=0*ones(busNum,1);
% I1r=0.1*ones(length(Loadi),1);%ע<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>෴
% I1i=0.1*ones(length(Loadi),1);%ע<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>෴
I1r=real(guessIf1);
I1i=imag(guessIf1);
KK=0;
plotGap=zeros(1,60);
%<EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD>
%״̬<EFBFBD><EFBFBD>Ϊ SEPD SEQD SEVmf1 SEVaf1
RestraintCount=length(Loadi)*2;
% ContrlCount=busNum*2+length(Loadi)*2;
noBoundedLoadi=setdiff(Loadi,noLoadi);
noBoundedLoadi=[];
ContrlCount=busNum*2+length(Loadi)*2;
CenterA=0.1;
Init_Z=1*sparse(ones(RestraintCount,1));
Init_W=sparse(-1*ones(RestraintCount,1));
Init_L=1*sparse(ones(RestraintCount,1));
Init_U=1*sparse(ones(RestraintCount,1));
Init_Y=sparse(2*busNum,1);%<EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Gap=(Init_L'*Init_Z-Init_U'*Init_W);
kmax=100;
isSetBound=0;
largerBound=1;
realBound=0;
tic
while(abs(Gap)>1e-4)
    if KK>=kmax
        break;
    end
    plotGap(KK+1)=Gap;
    Init_u=Gap/2/RestraintCount*CenterA;
    %% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>OPF
    %% <EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
    deltH=func_deltH(busNum,fsY1,Loadi,Balance);
    %% <EFBFBD>γɲ<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
     deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
    %%
%     L_1Z=diag(Init_Z./Init_L);
%     U_1W=diag(Init_W./Init_U);
    %% <EFBFBD>γɺ<EFBFBD>ɭ<EFBFBD><EFBFBD>
    deltdeltF=func_deltdeltF(busNum,fsY1,Loadi,wV1r,wV1i,wI1r,wI1i,V1measurement,V1r,V1i);
%     deltdeltF=0;
    %% <EFBFBD>γ<EFBFBD>ddHy
%     ddh=func_ddh(busNum,Loadi,Init_Z,Init_W);
    ddh=0;
    %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ddg
    ddg=0;
    %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>deltF
    deltF=func_deltF(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
%     deltF=0;
    %%
    Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
    Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
    Mat_G=FormGstate1(I1r,I1i);
    Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
    Ly=Mat_H;
    
%     if isSetBound==0
    
        Ir=real(guessIf1);
        pIr=find(Ir>0);
        nIr=find(Ir<0);
        Ii=imag(guessIf1);
        pIi=find(Ii>0);
        nIi=find(Ii<0);
        lower=ones(length(Loadi)*2,1);
%     if abs(Gap)<1 || isSetBound==1

    if abs(Gap)<1000
        largerBound=0;
        realBound=1;
    end

     if realBound==1
%         lower(pIr)=(0.9-3/(KK+1))*Ir(pIr);
%         lower(nIr)=(1.100+3/(KK+1))*Ir(nIr);
%         lower(pIi+length(Ir))=(0.9-3/(KK+1))*Ii(pIi);
%         lower(nIi+length(Ir))=(1.10+3/(KK+1))*Ii(nIi);
        isSetBound=1;
        lower(pIr)=(0.7)*Ir(pIr);
        lower(nIr)=(1.300)*Ir(nIr);
        lower(pIi+length(Ir))=(0.7)*Ii(pIi);
        lower(nIi+length(Ir))=(1.30)*Ii(nIi);
        lower( ismember(Loadi,noBoundedLoadi))=-0.9;
        lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-0.9;
     end
        
    if largerBound==1
        lower(pIr)=-2.8;
        lower(nIr)=-2.8;
        lower(pIi+length(Ir))=-2.8;
        lower(nIi+length(Ir))=-2.8;
        
%         lower(pIr)=(0.7)*Ir(pIr)-0.1;
%         lower(nIr)=(1.300)*Ir(nIr)-0.1;
%         lower(pIi+length(Ir))=(0.7)*Ii(pIi)-0.1;
%         lower(nIi+length(Ir))=(1.30)*Ii(nIi)-0.1;
     end

     
     
%     lowerR=lower(1:length(lower)/2);
%     lowerR(I1r./lowerR>0.98)=lowerR(I1r./lowerR>0.98)-0.001;
%     lowerI=lower(length(lower)/2+1:end);
%     lowerI(I1i./lowerI>0.98)=lowerI(I1i./lowerI>0.98)-0.001;    
%     lower=[lowerR;lowerI];

%         I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998);
% I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998);
    
        Ir=real(guessIf1);
        pIr=find(Ir>0);
        nIr=find(Ir<0);
        Ii=imag(guessIf1);
        pIi=find(Ii>0);
        nIi=find(Ii<0);
        upper=ones(length(Loadi)*2,1);
%     if abs(Gap)<1 || isSetBound==1
   if realBound==1
%         upper(pIr)=(1.10+3/(KK+1))*Ir(pIr);
%         upper(nIr)=(0.9-3/(KK+1))*Ir(nIr);
%         upper(pIi+length(Ir))=(1.10+3/(KK+1))*Ii(pIi);
%         upper(nIi+length(Ir))=(0.9-3/(KK+1))*Ii(nIi);
        isSetBound=1;
        upper(pIr)=(1.30)*Ir(pIr);
        upper(nIr)=(0.7)*Ir(nIr);
        upper(pIi+length(Ir))=(1.30)*Ii(pIi);
        upper(nIi+length(Ir))=(0.7)*Ii(nIi);
        upper(ismember(Loadi,noBoundedLoadi))=0.9;
        upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=0.9;
   end

    if largerBound==1

        upper(pIr)=2.8;
        upper(nIr)=2.8;
        upper(pIi+length(Ir))=2.8;
        upper(nIi+length(Ir))=2.8;

%         upper(pIr)=(1.30)*Ir(pIr)+0.1;
%         upper(nIr)=(0.7)*Ir(nIr)+0.1;
%         upper(pIi+length(Ir))=(1.30)*Ii(pIi)+0.1;
%         upper(nIi+length(Ir))=(0.7)*Ii(nIi)+0.1;
    end
        
        
%     end
    
    Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag,guessIf1,rIf1,lower);
    Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag,guessIf1,rIf1,upper);
    Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
    YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
    %% <EFBFBD><EFBFBD>ʼ<EFBFBD>ⷽ<EFBFBD><EFBFBD>
    XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,busNum);
    %%ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum,Balance);
    [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,Vref);
    Gap=(Init_L'*Init_Z-Init_U'*Init_W)+max(abs(deltX));
    fprintf('Gap %f :%d\n',full(Gap),KK+1);
    KK=KK+1;
end
toc
calcuTime=toc/KK;
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
if abs(Gap)<0.00001 && KK<kmax
    isConverged=1;
else
    isConverged=0;
end
end