diff --git a/OPF.m b/OPF.m index c60a8b9..df310c9 100644 --- a/OPF.m +++ b/OPF.m @@ -28,7 +28,7 @@ mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; % mVolt=load('mVolt'); % mVolt=mVolt.mVolt; -mVolt(2)=rVolt(2)*(1-sigma*4); +mVolt(2)=rVolt(2)*(1+0.2); %% 画Case A的图 % figure('Color',[1 1 1]); [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有 diff --git a/subOPF.m b/subOPF.m index 76510d4..d0007b2 100644 --- a/subOPF.m +++ b/subOPF.m @@ -90,14 +90,15 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);% uVolt(noMeasurei)=0.7*mVolt(noMeasurei); %错误数据 %mVolt(2)=5; -bigM=2; -Vbi=sparse(0.9*ones(Busnum,1)); -PDbi=sparse(0.9*ones(length(Loadi),1)); -QDbi=sparse(0.9*ones(length(Loadi),1)); -eps=10; +bigM=0.3; +Vbi=sparse(0.0*ones(Busnum,1)); +% Vbi(2)=1; +PDbi=sparse(0.0*ones(length(Loadi),1)); +QDbi=sparse(0.0*ones(length(Loadi),1)); +eps=1; % 第一遍,算连续的值 fprintf('第1次迭代,算连续量。\n'); -while(abs(Gap)>Precision*1) +while(abs(Gap)>Precision*10) if KK>kmax break; end @@ -123,6 +124,7 @@ while(abs(Gap)>Precision*1) %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); + 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_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Ly=Mat_H; @@ -142,13 +144,13 @@ end % 第二遍,算离散 % [~,~,Init_Z,Init_W,Init_L,Init_U,~,~,~,RestraintCount,~,~,~,~,~,~,~,~,~]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD); % Gap=(Init_L*Init_Z'-Init_U*Init_W'); -Gap=1000; +Gap=1; % KK=0; eps=1; fprintf('\n'); fprintf('第2次迭代,算离散量。\n'); -while eps>0.001 - while(abs(Gap)>Precision*10) +% while eps>0.0001 + while(abs(Gap)>Precision*10 || eps>0.0001) if KK>kmax break; end @@ -174,9 +176,15 @@ while eps>0.001 %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); + bigM=0.2; 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; + if eps>0.01 + eps=Gap*2 + else + eps=eps*0.2 + end % if KK>20 % eps=eps*0.1; % eps=Gap; @@ -184,9 +192,9 @@ while eps>0.001 % eps=1e-5; % end % eps; - % % if any(Vbi>0.002) - % % Vbi(Vbi>0.002)=1; - % % end +% if any(Vbi>0.5) +% Vbi(Vbi>0.5)=1; +% end % end 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); @@ -204,9 +212,10 @@ while eps>0.001 Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=KK+1; end - eps=eps*0.9 - Gap=100; -end +% eps=eps*0.95; +% eps=Gap +% Gap=100; +% end %% 计算最大不平衡量 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';