facat
/
pes2014-basic
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

又调了一下，更理想了。 

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2015-02-02 21:28:56 +08:00
parent a2b2ba134f
commit beaa6855d4
1 changed files with 56 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
subOPF.m
View File

@ -90,7 +90,7 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);%
uVolt(noMeasurei)=0.7*mVolt(noMeasurei); uVolt(noMeasurei)=0.7*mVolt(noMeasurei);
%错误数据 %错误数据
%mVolt(2)=5; %mVolt(2)=5;
bigM=0.000003; % bigM=0.000003;
Vbi=sparse(0.5*ones(Busnum,1)); Vbi=sparse(0.5*ones(Busnum,1));
% Vbi(2)=1; % Vbi(2)=1;
PDbi=sparse(0.5*ones(length(Loadi),1)); PDbi=sparse(0.5*ones(length(Loadi),1));
@ -124,7 +124,7 @@ while(abs(Gap)>Precision*10)
%% 形成方程矩阵 %% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
bigM=0.7; bigM=1;
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H; Ly=Mat_H;
@ -146,69 +146,70 @@ end
% Gap=(Init_L*Init_Z'-Init_U*Init_W'); % Gap=(Init_L*Init_Z'-Init_U*Init_W');
Gap=1; Gap=1;
% KK=0; % KK=0;
eps=0.00001; % eps=0.00001;
fprintf('\n'); fprintf('\n');
fprintf('第2次迭代，算离散量。\n'); fprintf('第2次迭代，算离散量。\n');
while eps>0.000001 % while eps>0.000001
% Init_Z=sparse(ones(1,RestraintCount)); % Init_Z=sparse(ones(1,RestraintCount));
% Init_W=sparse(-1*ones(1,RestraintCount)); % Init_W=sparse(-1*ones(1,RestraintCount));
% Init_L=1*sparse(ones(1,RestraintCount)); % Init_L=1*sparse(ones(1,RestraintCount));
% Init_U=1*sparse(ones(1,RestraintCount)); % Init_U=1*sparse(ones(1,RestraintCount));
% Init_Y=sparse(1,2*Busnum);%与学姐一致 % Init_Y=sparse(1,2*Busnum);%与学姐一致
while(abs(Gap)>Precision*100) while(abs(Gap)>Precision*10)
if KK>kmax if KK>kmax
% break; % break;
end
plotGap(KK+1)=Gap;
Init_u=Gap/2/RestraintCount*CenterA;
AngleIJMat=0;
%% 开始计算OPF
%% 形成等式约束的雅克比
deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
%% 形成不等式约束的雅克比
deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
%%
L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U);
%% 形成海森阵
deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
%% 形成ddHy
ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
%% 开始构建ddg
ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
%% 开始构建deltF
deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
%% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
bigM=0.7;
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%% 开始解方程
fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
%%取各分量
% Vbi_=Vbi;
% PDbi_=PDbi;
% QDbi_=QDbi;
% Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1;
[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
% Vbi(Vbi>0.01)=1;
KK=KK+1;
end end
eps=eps*0.4; plotGap(KK+1)=Gap;
% eps=Gap Init_u=Gap/2/RestraintCount*CenterA;
Gap=100; AngleIJMat=0;
%% 开始计算OPF
%% 形成等式约束的雅克比
deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
%% 形成不等式约束的雅克比
deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
%%
L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U);
%% 形成海森阵
deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
%% 形成ddHy
ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
%% 开始构建ddg
ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
%% 开始构建deltF
deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
%% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
bigM=0.7;
eps=Gap*0.1;
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%% 开始解方程
fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
%%取各分量
% Vbi_=Vbi;
% PDbi_=PDbi;
% QDbi_=QDbi;
% Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1;
[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
% Vbi(Vbi>0.01)=1;
KK=KK+1;
end end
% eps=Gap
% Gap=100;
% end
%% 计算最大不平衡量 %% 计算最大不平衡量
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt'; dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';