135 lines
4.3 KiB
Matlab
135 lines
4.3 KiB
Matlab
function [ V1r,V1i,I1r,I1i,isConverged,calcuTime,KK ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref,dI_F,flag,guessIf1,rIf1,noLoadi )
|
|
|
|
%把每个序的循环写在这个函数中。其实也就是内点法循环。
|
|
V1r=Vref*ones(busNum,1);
|
|
V1i=0*ones(busNum,1);
|
|
% I1r=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
|
|
% I1i=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
|
|
I1r=real(guessIf1);
|
|
I1i=imag(guessIf1);
|
|
KK=0;
|
|
plotGap=zeros(1,60);
|
|
|
|
% ContrlCount=busNum*2+length(Loadi)*2;
|
|
noBoundedLoadi=setdiff(Loadi,noLoadi);
|
|
% noBoundedLoadi=[];
|
|
boundedLoadi=setdiff(Loadi,noBoundedLoadi);
|
|
%初始化
|
|
%状态量为 SEPD SEQD SEVmf1 SEVaf1
|
|
% RestraintCount=length(Loadi)*2;
|
|
RestraintCount=length(boundedLoadi)*2;
|
|
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);%等式约束乘子
|
|
Gap=(Init_L'*Init_Z-Init_U'*Init_W);
|
|
kmax=15;
|
|
isSetBound=0;
|
|
largerBound=0;
|
|
realBound=1;
|
|
|
|
|
|
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
|
|
pIr=find(Ir>0);
|
|
nIr=find(Ir<0);
|
|
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
|
|
pIi=find(Ii>0);
|
|
nIi=find(Ii<0);
|
|
lower=ones(length(boundedLoadi)*2,1);
|
|
if realBound==1
|
|
isSetBound=1;
|
|
lower(pIr)=(0.7)*Ir(pIr);
|
|
lower(nIr)=(1.3)*Ir(nIr);
|
|
lower(pIi+length(Ir))=(0.7)*Ii(pIi);
|
|
lower(nIi+length(Ir))=(1.3)*Ii(nIi);
|
|
|
|
% lower(pIr)=-0.1;
|
|
% lower(nIr)=-0.1;
|
|
% lower(pIi+length(Ir))=-0.1;
|
|
% lower(nIi+length(Ir))=-0.1;
|
|
|
|
end
|
|
|
|
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
|
|
pIr=find(Ir>0);
|
|
nIr=find(Ir<0);
|
|
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
|
|
pIi=find(Ii>0);
|
|
nIi=find(Ii<0);
|
|
upper=ones(length(boundedLoadi)*2,1);
|
|
% if abs(Gap)<1 || isSetBound==1
|
|
if realBound==1
|
|
isSetBound=1;
|
|
upper(pIr)=(1.3)*Ir(pIr);
|
|
upper(nIr)=(0.7)*Ir(nIr);
|
|
upper(pIi+length(Ir))=(1.3)*Ii(pIi);
|
|
upper(nIi+length(Ir))=(0.7)*Ii(nIi);
|
|
% upper(pIr)=0.1;
|
|
% upper(nIr)=0.1;
|
|
% upper(pIi+length(Ir))=0.1;
|
|
% upper(nIi+length(Ir))=0.1;
|
|
end
|
|
tic
|
|
while(abs(Gap)>1e-5)
|
|
if KK>=kmax
|
|
break;
|
|
end
|
|
plotGap(KK+1)=Gap;
|
|
Init_u=Gap/2/RestraintCount*CenterA;
|
|
%% 开始计算OPF
|
|
%% 形成等式约束的雅克比
|
|
deltH=func_deltH(busNum,fsY1,Loadi,Balance);
|
|
%% 形成不等式约束的雅克比
|
|
% deltG=func_deltGstate1(busNum,Loadi,I1r,I1i);
|
|
deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
|
|
%%
|
|
% L_1Z=diag(Init_Z./Init_L);
|
|
% U_1W=diag(Init_W./Init_U);
|
|
%% 形成海森阵
|
|
deltdeltF=func_deltdeltF(busNum,fsY1,Loadi,wV1r,wV1i,wI1r,wI1i,V1measurement,V1r,V1i);
|
|
% deltdeltF=0;
|
|
%% 形成ddHy
|
|
% ddh=func_ddh(busNum,Loadi,Init_Z,Init_W);
|
|
ddh=0;
|
|
%% 开始构建ddg
|
|
ddg=0;
|
|
%% 开始构建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_G=FormGstate1(I1r,I1i,ismember(Loadi,boundedLoadi) );
|
|
Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
|
|
Ly=Mat_H;
|
|
|
|
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);
|
|
%% 开始解方程
|
|
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);
|
|
%%取各分量
|
|
[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));
|
|
% Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]);
|
|
fprintf('Gap %f :%d\n',full(Gap),KK+1);
|
|
KK=KK+1;
|
|
end
|
|
calcuTime=0;
|
|
% toc
|
|
calcuTime=toc;
|
|
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
|
|
if abs(Gap)<1e-5 && KK<kmax
|
|
isConverged=1;
|
|
else
|
|
isConverged=0;
|
|
end
|
|
end
|
|
|