加测试程序代码，并加入评价指标代码。

Signed-off-by: facat@lab.com <facat@lab.com>

FormLw.m

@@ -26,8 +26,8 @@ function Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P)
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
 % PDU=1.8*ones(length(Loadi),1);
 % QDU=1.8*ones(length(Loadi),1);
-PDU=1.05*rPD3P;
+PDU=1.1*rPD3P;
-QDU=1.05*rQD3P;
+QDU=1.1*rQD3P;
 t1=([PDU',QDU'])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLz.m

@@ -23,8 +23,8 @@ function Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P)
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
 % PDL=0*ones(length(Loadi),1);
 % QDL=0*ones(length(Loadi),1);
-PDL=0.95*rPD3P;
+PDL=0.9*rPD3P;
-QDL=0.95*rQD3P;
+QDL=0.9*rQD3P;
 t1=([PDL',QDL'])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

OPF.m

@@ -1,16 +1,18 @@
+
+function [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF()
 tic
 clc
 clear
-lineZ=readLineZ('feeder123\lineParameter.txt');
+lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'feeder123\data.txt');
+    cap]=dataRead(lineZ,'feeder13\data1.txt');
-phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
+% phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
-phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
+% phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
-phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
+% phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
-phaseASpotLoadQ(phaseASpotLoadQ==0)=0.002;
+% phaseASpotLoadQ(phaseASpotLoadQ==0)=0.002;
-phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002;
+% phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002;
-phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
+% phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
 
 
 
@@ -151,12 +153,12 @@ PQi3P(3:3:end)=(PQi-1)*3+3;
 Loadi=PQi3P;
 PD3P=sparse(Busnum*3,1);
 QD3P=sparse(Busnum*3,1);
-PD3P(1:3:end)=phaseASpotLoadP*0.5;
+PD3P(1:3:end)=phaseASpotLoadP*0.9;
-PD3P(2:3:end)=phaseBSpotLoadP*0.5;
+PD3P(2:3:end)=phaseBSpotLoadP*0.9;
-PD3P(3:3:end)=phaseCSpotLoadP*0.5;
+PD3P(3:3:end)=phaseCSpotLoadP*0.9;
-QD3P(1:3:end)=phaseASpotLoadQ*0.5;
+QD3P(1:3:end)=phaseASpotLoadQ*0.9;
-QD3P(2:3:end)=phaseBSpotLoadQ*0.5;
+QD3P(2:3:end)=phaseBSpotLoadQ*0.9;
-QD3P(3:3:end)=phaseCSpotLoadQ*0.5;
+QD3P(3:3:end)=phaseCSpotLoadQ*0.9;
 PD3P=PD3P(Loadi);
 QD3P=QD3P(Loadi);
 QGi3P=zeros(length(QGi)*3,1);
@@ -192,10 +194,10 @@ rVoltABCV(1:3:end)=abs(VoltpABC(1,:));
 rVoltABCV(2:3:end)=abs(VoltpABC(2,:));
 rVoltABCV(3:3:end)=abs(VoltpABC(3,:));
 %三相电压相角
-rVoltABCA=zeros(busNum*3,1);
+% rVoltABCA=zeros(busNum*3,1);
-rVoltABCA(1:3:end)=angle(VoltpABC(1,:));
+% rVoltABCA(1:3:end)=angle(VoltpABC(1,:));
-rVoltABCA(2:3:end)=angle(VoltpABC(2,:));
+% rVoltABCA(2:3:end)=angle(VoltpABC(2,:));
-rVoltABCA(3:3:end)=angle(VoltpABC(3,:));
+% rVoltABCA(3:3:end)=angle(VoltpABC(3,:));
 %三相负荷
 rPD3P=zeros(busNum*3,1);
 rPD3P(1:3:end)=phaseASpotLoadP;
@@ -216,9 +218,12 @@ mPD3P=rPD3P.*PD3PSigma;
 QD3PSigma=(1+normrnd(0,sigma,length(rQD3P),1));
 mQD3P=rPD3P.*QD3PSigma;
 %量测方差
-wVolt=1./(abs(mVoltABCV*sigma).^2)*1;
+wVolt=1./(abs(mVoltABCV*sigma).^2);
-wPD=1./(abs(mPD3P*sigma).^2)*1;
+wPD=1./(abs(mPD3P*sigma).^2);
-wQD=1./(abs(mQD3P*sigma).^2)*1;
+wQD=1./(abs(mQD3P*sigma).^2);
+wVolt(1:3:end)=0;
+wPD(1:3:end)=0;
+wQD(1:3:end)=0;
 % wVolt=1./abs(VoltSigma).^2;
 % wPD=1./abs(PD3PSigma).^2;
 % wQD=1./abs(QD3PSigma).^2;
@@ -244,15 +249,15 @@ while(abs(Gap)>Precision)
         break;
     end
     Init_u=Gap/2/RestraintCount*CenterA;
-    AngleIJMat=0;
+%     AngleIJMat=0;
     %% 开始计算OPF
     %% 形成等式约束的雅克比
     deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi);
     %% 形成不等式约束的雅克比
     deltG=func_deltG(Busnum,Loadi);
     %%
-    L_1Z=diag(Init_Z./Init_L);
+%     L_1Z=diag(Init_Z./Init_L);
-    U_1W=diag(Init_W./Init_U);
+%     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
     deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount);
 %         deltdeltF=0;
@@ -284,7 +289,26 @@ while(abs(Gap)>Precision)
     fprintf('%f\n',full(Gap));
     KK=KK+1;
 end
-(rVoltABCV-Volt)./rVoltABCV*100
+(rVoltABCV-Volt)./rVoltABCV*100;
+%% 计算统计量
+%目标函数均值
+JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2);
+%估计误差统计
+%量测量数量
+mCount=Busnum*3+length(Loadi)*3*2;
+%估计量质量
+%有PD，QD的
+SEMeasurement=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((rPD3P-PD3P)./mPD3P./sigma).^2)+sum(((rQD3P-QD3P)./mQD3P./sigma).^2);
+SEMeasurement=(SEMeasurement/mCount)^.5;
+%没有PD，QD的
+SEMeasurementNOPQ=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2);
+SEMeasurementNOPQ=(SEMeasurementNOPQ/mCount)^.5;
+%量测量质量
+SEMeasurementM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2)+sum(((rPD3P-mPD3P)./mPD3P./sigma).^2)+sum(((rQD3P-mQD3P)./mQD3P./sigma).^2);
+SEMeasurementM=(SEMeasurementM/mCount)^.5;
+%没有PD，QD的
+SEMeasurementNOPQM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2);
+SEMeasurementNOPQM=(SEMeasurementNOPQM/mCount)^.5;
 fprintf('迭代次数%d\n',KK);
 toc
-
+end

test.m

@@ -0,0 +1,22 @@
+close all
+clear
+clc
+JMeasurementSum=0;
+SEMeasurementSum=0;
+SEMeasurementNOPQSum=0;
+SEMeasurementMSum=0;
+SEMeasurementNOPQMSum=0;
+N=1000;
+for I=1:N
+    [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF();
+    JMeasurementSum=JMeasurement+JMeasurementSum;
+    SEMeasurementSum=SEMeasurement+SEMeasurementSum;
+    SEMeasurementNOPQSum=SEMeasurementNOPQ+SEMeasurementNOPQSum;
+    SEMeasurementMSum=SEMeasurementMSum+SEMeasurementM;
+    SEMeasurementNOPQMSum=SEMeasurementNOPQMSum+SEMeasurementNOPQM;
+end
+JMeasurementSum=JMeasurementSum/N;
+SEMeasurementSum=SEMeasurementSum/N;
+SEMeasurementNOPQSum=SEMeasurementNOPQSum/N;
+SEMeasurementMSum=SEMeasurementMSum/N;
+SEMeasurementNOPQMSum=SEMeasurementNOPQMSum/N;