From 72a068c4ff596111ba1e00004de63924bbe67c1e Mon Sep 17 00:00:00 2001
From: "dugg@lab-desk" <dugg@lab-desk>
Date: Sat, 24 Jan 2015 22:27:51 +0800
Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E6=94=B9=E4=BA=86=E4=B8=80=E4=B8=8B?=
 =?UTF-8?q?=E7=BB=9F=E8=AE=A1=E6=8C=87=E6=A0=87=EF=BC=8C=E6=8C=89=E4=BA=8E?=
 =?UTF-8?q?=E5=B0=94=E9=93=BF=E7=9A=84=E4=B9=A6=E4=B8=8A=E5=85=AC=E5=BC=8F?=
 =?UTF-8?q?=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
---
 OPF.m | 36 +++++++++++++++++++++++++++---------
 1 file changed, 27 insertions(+), 9 deletions(-)

diff --git a/OPF.m b/OPF.m
index bbb134f..d38f1f3 100644
--- a/OPF.m
+++ b/OPF.m
@@ -2,7 +2,10 @@ clc
 clear
 close all
 arrayA=zeros(21,10);
-for I=1:1
+sumCaseA_SE=0;
+sumCaseB_SE=0;
+sumCaseC_SE=0;
+for I=1:20
 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');
@@ -18,12 +21,12 @@ QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
 mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
 
 %加载保存的变量
-PD0=load('PD0');
-PD0=PD0.PD0;
-QD0=load('QD0');
-QD0=QD0.QD0;
-mVolt=load('mVolt');
-mVolt=mVolt.mVolt;
+% PD0=load('PD0');
+% PD0=PD0.PD0;
+% QD0=load('QD0');
+% QD0=QD0.QD0;
+% mVolt=load('mVolt');
+% mVolt=mVolt.mVolt;
 
 % mVolt(19)=rVolt(19)*(1+sigma*4);
 %% 画Case A的图
@@ -64,11 +67,13 @@ ylabel('
 %计算Case A的误差
 CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
 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(21,I)=CaseAE*1000;
 
-end
 %% 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);%
@@ -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');
 CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
 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
 [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);
@@ -135,11 +146,18 @@ ld=legend('
 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);
 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('Case A    Case B    Case C 的误差\n')
 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')