给投节点很大的权重，就是想作为真实值用。头节点没有误差。

Signed-off-by: dmy@lab <dmy@lab.lab>

OPF.m

@@ -3,10 +3,10 @@ function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_m
 tic
 clc
 clear
-lineZ=readLineZ('feeder13\lineParameter.txt');
+lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'feeder13\data1.txt');
+    cap]=dataRead(lineZ,'./../DistributionNetwork-Power2Current/modified-feeder69\data.txt');
 % phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
 % phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
 % phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
@@ -253,6 +253,7 @@ noLoadi=[1,8,10,12];
 % noLoadi=Loadi;
 % noLoadi=[1,11];
 noLoadi=[5];
+noLoadi=[];
 % noLoadi=[5,11,19,27,37,43,54,61];
 noPQi3P=zeros(length(noLoadi)*3,1);
 noPQi3P(1:3:end)=(noLoadi-1)*3+1;
@@ -268,12 +269,19 @@ PD3PSigma=(1+normrnd(0,sigma,length(rPD3P),1));
 mPD3P=rPD3P.*PD3PSigma;
 QD3PSigma=(1+normrnd(0,sigma,length(rQD3P),1));
 mQD3P=rQD3P.*QD3PSigma;
+
+%ÓÃƽºâ½ÚµãÕæʵֵ
+mPD3P(ismember(Loadi,Balance3P))=[-real(PGQG((Balance-1)*3+1));-real(PGQG((Balance-1)*3+2));-real(PGQG((Balance-1)*3+3))];
+mQD3P(ismember(Loadi,Balance3P))=[-imag(PGQG((Balance-1)*3+1));-imag(PGQG((Balance-1)*3+2));-imag(PGQG((Balance-1)*3+3))];
+
 mPD3P(ismember(Loadi,noPQi3P))=rPD3P(ismember(Loadi,noPQi3P)).*(1+unifrnd(-0.15,0.15,length(noPQi3P),1));
 mQD3P(ismember(Loadi,noPQi3P))=rQD3P(ismember(Loadi,noPQi3P)).*(1+unifrnd(-0.15,0.15,length(noPQi3P),1));
 %量测方差
 wVolt=1./(abs(mVoltABCV*sigma/3).^2);
 wPD=1./(abs(mPD3P*sigma).^2);
+wPD(ismember(Loadi,Balance3P))=1e20;
 wQD=1./(abs(mQD3P*sigma).^2);
+wQD(ismember(Loadi,Balance3P))=1e20;
 wVolt(setdiff(1:length(wVolt),Loadi))=0;%只有负荷处才有电压量测。
 wVolt( noPQi3P)=0;
 wPD(ismember( Loadi,noPQi3P))=1./(abs(mPD3P(ismember( Loadi,noPQi3P))*0.15).^2);
@@ -305,12 +313,12 @@ while(maxD>Precision)
     %% 形成不等式约束的雅克比
     %%  开始构建ddg
     %%  开始构建deltF
-    deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi);
+%     deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi);
     %% 形成方程矩阵
     Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
     %% 开始解方程
 %         fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
-    XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
+    XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
     %%取各分量
     [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum);
     [Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
@@ -338,6 +346,10 @@ AME_QD=sum(sum(abs(rQD3P-QD3P)))/length(Loadi);
 AME_mVolt=sum(sum(abs( mVoltABCV-rVoltABCV)))/3/busNum;
 AME_mPD=sum(sum(abs(rPD3P-mPD3P)))/length(Loadi);
 AME_mQD=sum(sum(abs(rQD3P-mQD3P)))/length(Loadi);
+
+% AME_mPD=sum(sum(abs(rPD3P(1:144)-mPD3P(1:144))))/144;
+% AME_mQD=sum(sum(abs(rQD3P(1:144)-mQD3P(1:144))))/144;
+
 %返回收敛信息
 isConverged=1;
 if KK>=kmax

SolveIt.m

@@ -1,4 +1,4 @@
-function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
+function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
 % dH_dx
 dP_x=-dH_dx(Loadi,:);
 dQ_x=-dH_dx(Loadi+Busnum*3,:);

jacobian_M.m

@@ -8,8 +8,8 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
 AngleIJ=UAngel(r)-UAngel(c)-Angle;
 mat_AngleIJ=sparse(r,c,AngleIJ,Busnum*3,Busnum*3);
-H=sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+H=diag(Volt)*Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-N=-sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+N=-diag(Volt)*Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
 J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
 L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
 t1=[J,H;

test.m

@@ -9,7 +9,7 @@ AME_QDSum=0;
 AME_mVoltSum=0;
 AME_mPDSum=0;
 AME_mQDSum=0;
-totalTimeSum=1000;
+totalTimeSum=0;
 N=200;
 loopN=1;
 while 1
@@ -30,7 +30,7 @@ while 1
     AME_mVoltSum=AME_mVoltSum+AME_mVolt;
     AME_mPDSum=AME_mPDSum+AME_mPD;
     AME_mQDSum=AME_mQDSum+AME_mQD;
-    totalTimeSum=min([totalTimeSum,totalTime]);
+    totalTimeSum=totalTimeSum+totalTime;
 end
 JMeasurementSum=JMeasurementSum/N;
 AME_VoltSum=AME_VoltSum/N;
@@ -40,3 +40,4 @@ AME_QDSum=AME_QDSum/N;
 AME_mVoltSum=AME_mVoltSum/N;
 AME_mPDSum=AME_mPDSum/N;
 AME_mQDSum=AME_mQDSum/N;
+totalTimeSum=totalTimeSum/N;