diff --git a/FormLxComa.m b/FormLxComa.m index 9d25a12..918ef87 100644 --- a/FormLxComa.m +++ b/FormLxComa.m @@ -4,8 +4,9 @@ t2=Lul+diag(Init_Z)*Lz; t3=inv(diag(Init_L)); t4=t3*t2;% t5=Luu-diag(Init_W)*Lw; -t6=inv(diag(Init_U)); -t7=t6*t5;% +% t6=inv(diag(Init_U)); +% t7=t6*t5;% +t7=diag(Init_U)\t5; t8=deltG*(t4+t7);%% LxComa=Lx+t8; end \ No newline at end of file diff --git a/OPF.m b/OPF.m index d38f1f3..c60a8b9 100644 --- a/OPF.m +++ b/OPF.m @@ -5,7 +5,7 @@ arrayA=zeros(21,10); sumCaseA_SE=0; sumCaseB_SE=0; sumCaseC_SE=0; -for I=1:20 +for I=1:1 close [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ... pf('E:\算例\feeder33\feeder33.txt'); @@ -28,14 +28,18 @@ mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; % mVolt=load('mVolt'); % mVolt=mVolt.mVolt; -% mVolt(19)=rVolt(19)*(1+sigma*4); +mVolt(2)=rVolt(2)*(1-sigma*4); %% 画Case A的图 -figure('Color',[1 1 1]); +% figure('Color',[1 1 1]); [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有 +end subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]); CaseAREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case A % CaseAREV=CaseAREV(2:end)*100; +%真实值的 plot(1:length(CaseAREV),abs(CaseAREV),'k.:','Marker','diamond'); +%测量值的 +% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond'); box off; set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) xlabel('节点号'); @@ -43,6 +47,7 @@ ylabel(' subplot(4,1,2); CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A % CaseAREA=CaseAREA(2:end)*100; +%真实值的 plot(1:length(CaseAREA),abs(CaseAREA),'k:','Marker','diamond'); box off; set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) @@ -51,7 +56,10 @@ ylabel(' subplot(4,1,3); CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A % CaseAREP=CaseAREP(2:end)*100; +%真实值的 plot(1:length(CaseAREP),abs(CaseAREP),'k:','Marker','diamond'); +%测量值的 +% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond'); box off; set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) xlabel('节点号'); @@ -59,7 +67,10 @@ ylabel(' subplot(4,1,4); CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A % CaseAREQ=CaseAREQ(2:end)*100; +%真实值的 plot(1:length(CaseAREQ),abs(CaseAREQ),'k:','Marker','diamond'); +%测量值的 +% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond'); box off; set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) xlabel('节点号'); @@ -73,7 +84,6 @@ CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5; sumCaseA_SE=sumCaseA_SE+CaseA_SE; % arrayA(1:19,I)=Vbi; % arrayA(21,I)=CaseAE*1000; - %% Case B % [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);% @@ -156,8 +166,6 @@ fprintf('Case A Case B Case C fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE); fprintf('三个Case目标值\n') fprintf('%f\t%f\t%f \n',objA,objB,objC) - -end %% 画测量值 % subplot(4,1,1); % plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram') @@ -309,3 +317,10 @@ fz=find(abs(plotGapA)==0); fz=fz(1); figure('Name','互补曲线') plot(1:fz-1,plotGapA(1:fz-1)); +figure('Name','最大不平衡量'); +%% 最大不平衡量 +% maxDismatchPQ = [0.3123e-10 0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10]; +% bar(maxDismatchPQ); +%% 计算时间 +calTime=[87 61; 211 87; 108 59;]; +bar(calTime,'stacked'); diff --git a/subOPF.m b/subOPF.m index 352cce6..76510d4 100644 --- a/subOPF.m +++ b/subOPF.m @@ -45,7 +45,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=0; plotGap=zeros(1,60); ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; -kmax=100; +kmax=100000; Precision=Precision/1; %% 加误差 %找DG @@ -90,10 +90,10 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);% uVolt(noMeasurei)=0.7*mVolt(noMeasurei); %错误数据 %mVolt(2)=5; -bigM=10; -Vbi=sparse(0.5*ones(Busnum,1)); -PDbi=sparse(0.5*ones(length(Loadi),1)); -QDbi=sparse(0.5*ones(length(Loadi),1)); +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; % 第一遍，算连续的值 fprintf('第1次迭代，算连续量。\n'); @@ -144,62 +144,68 @@ end % Gap=(Init_L*Init_Z'-Init_U*Init_W'); Gap=1000; % KK=0; +eps=1; fprintf('\n'); fprintf('第2次迭代，算离散量。\n'); -while(abs(Gap)>Precision*1) - if KK>kmax - 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); - 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 KK>0 - eps=eps*0.1; - if abs(eps)<1e-6 - eps=1e-5; +while eps>0.001 + while(abs(Gap)>Precision*10) + if KK>kmax + break; end - eps; -% if any(Vbi>0.002) -% Vbi(Vbi>0.002)=1; -% 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); + 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 KK>20 + % eps=eps*0.1; + % eps=Gap; + % if abs(eps)<1e-6 + % eps=1e-5; + % end + % eps; + % % if any(Vbi>0.002) + % % Vbi(Vbi>0.002)=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); + 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; + [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'); + KK=KK+1; 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); - 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; - [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'); - KK=KK+1; + eps=eps*0.9 + Gap=100; end %% 计算最大不平衡量 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);