facat
/
pes2014-debug2-laplace
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-01-24 22:27:51 +08:00
parent 779c2d0e16
commit 72a068c4ff
1 changed files with 27 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
OPF.m
View File

@ -2,7 +2,10 @@ clc
clear clear
close all close all
arrayA=zeros(21,10); arrayA=zeros(21,10);
for I=1:1 sumCaseA_SE=0;
sumCaseB_SE=0;
sumCaseC_SE=0;
for I=1:20
close 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]= ... [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'); pf('E:\算例\feeder33\feeder33.txt');
@ -18,12 +21,12 @@ QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
%加载保存的变量 %加载保存的变量
PD0=load('PD0'); % PD0=load('PD0');
PD0=PD0.PD0; % PD0=PD0.PD0;
QD0=load('QD0'); % QD0=load('QD0');
QD0=QD0.QD0; % QD0=QD0.QD0;
mVolt=load('mVolt'); % mVolt=load('mVolt');
mVolt=mVolt.mVolt; % mVolt=mVolt.mVolt;
% mVolt(19)=rVolt(19)*(1+sigma*4); % mVolt(19)=rVolt(19)*(1+sigma*4);
%% 画Case A的图 %% 画Case A的图
@ -64,11 +67,13 @@ ylabel('
%计算Case A的误差 %计算Case A的误差
CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132); CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
notZeros=find(PD0~=0);
CaseA_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseA_SE=sumCaseA_SE+CaseA_SE;
% arrayA(1:19,I)=Vbi; % arrayA(1:19,I)=Vbi;
% arrayA(21,I)=CaseAE*1000; % arrayA(21,I)=CaseAE*1000;
end
%% Case B %% 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,[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);% [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);%
@ -94,6 +99,12 @@ CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
plot(1:length(CaseBREQ),abs(CaseBREQ),'b:','Marker','square'); plot(1:length(CaseBREQ),abs(CaseBREQ),'b:','Marker','square');
CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132); CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
noMeasurei=[2,3,5,20,24,27,28,10,11,12,13];
Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]);
CaseB_SE=sum(((PD0(Measurei)-PD(Measurei))./(PD0(Measurei)*sigma)).^2)+sum(((QD0(Measurei)-QD(Measurei))./(QD0(Measurei)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseB_SE=CaseB_SE+sum(((PD0(noMeasurei)-PD(noMeasurei))./(PD0(noMeasurei)*20*sigma)).^2)+sum(((QD0(noMeasurei)-QD(noMeasurei))./(QD0(noMeasurei)*20*sigma)).^2);
CaseB_SE=(CaseB_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseB_SE=sumCaseB_SE+CaseB_SE;
%% Case C %% Case C
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);%
subplot(4,1,1); subplot(4,1,1);
@ -135,11 +146,18 @@ ld=legend('
set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]); set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]);
CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132); CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
CaseC_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseC_SE=sumCaseC_SE+CaseC_SE;
% fprintf('目标函数值 %.2f\n',full(obj)); % fprintf('目标函数值 %.2f\n',full(obj));
fprintf('Case A Case B Case C 的误差\n') fprintf('Case A Case B Case C 的误差\n')
fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE); fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE);
fprintf('三个Case目标值\n') fprintf('三个Case目标值\n')
fprintf('%f\t%f\t%f \n',objA,objB,objC) fprintf('%f\t%f\t%f \n',objA,objB,objC)
end
%% 画测量值 %% 画测量值
% subplot(4,1,1); % subplot(4,1,1);
% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram') % plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')