1.把数据文件的电容删掉了。

2.负荷估计变量的维数与节点数不一致，只等于负荷数。
3.加了方差。

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

FormH.m

@@ -1,4 +1,4 @@
-function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle)
+function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
 %% 
 %QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
 %QD(QD~=0)=PD(QD~=0)./tan(QDcos);
@@ -6,8 +6,8 @@ function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle)
 %% 
 %PD(Loadi)=QD(Loadi)./tan(acos(0.98));
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum*3,Busnum*3);
-dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt;
+dP=PG-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt;
-dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt;
+dQ=QG-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt;
 
 Mat_H=[dP;dQ;];
 

Modification.m

@@ -14,8 +14,8 @@ Init_Y=Init_Y+AlphaD*deltY';
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 % QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 t=deltX(1:size(Loadi,1)*2);
-PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi));
+PD=PD+AlphaP*t(1:length(Loadi));
-QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2);
+QD=QD+AlphaP*t(length(Loadi)+1:length(Loadi)*2);
 t=deltX(size(Loadi,1)*2+1:ContrlCount)';
 t(Busnum*3+(Balance-1)*3+1)=0;
 t(Busnum*3+(Balance-1)*3+2)=0;

OPF.m

@@ -123,7 +123,6 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
 fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
 %% 潮流计算end
 
-
 busNum=length(nodeNum);
 Busnum=busNum;
 PQi=setxor(nodeNum,Balance);
@@ -139,14 +138,17 @@ PQi3P=zeros(length(PQi)*3,1);
 PQi3P(1:3:end)=(PQi-1)*3+1;
 PQi3P(2:3:end)=(PQi-1)*3+2;
 PQi3P(3:3:end)=(PQi-1)*3+3;
-PD3P=sparse(busNum*3,1);
+Loadi=PQi3P;
-QD3P=sparse(busNum*3,1);
+PD3P=sparse(Busnum*3,1);
+QD3P=sparse(Busnum*3,1);
 PD3P(1:3:end)=phaseASpotLoadP;
 PD3P(2:3:end)=phaseBSpotLoadP;
 PD3P(3:3:end)=phaseCSpotLoadP;
 QD3P(1:3:end)=phaseASpotLoadQ;
 QD3P(2:3:end)=phaseBSpotLoadQ;
 QD3P(3:3:end)=phaseCSpotLoadQ;
+PD3P=PD3P(Loadi);
+QD3P=QD3P(Loadi);
 QGi3P=zeros(length(QGi)*3,1);
 QGi3P(1:3:end)=(QGi-1)*3+1;
 QGi3P(2:3:end)=(QGi-1)*3+2;
@@ -168,13 +170,38 @@ Balance3P(3:3:end)=(Balance-1)*3+3;
 Vp3a((Balance-1)*3+1)=0;
 Vp3a((Balance-1)*3+2)=-120/180*pi;
 Vp3a((Balance-1)*3+3)=+120/180*pi;
-Loadi=PQi3P;
+%准备量测量和方差
-wVolt=ones(busNum*3,1);
+%真实值
-wPD=ones(length(Loadi),1);
+%三相电压幅值
-wQD=ones(length(Loadi),1);
+rVoltABCV=zeros(busNum*3,1);
-mVolt=Vp3m;
+rVoltABCV(1:3:end)=abs(VoltpABC(1,:));
-PD0=PD3P;
+rVoltABCV(2:3:end)=abs(VoltpABC(2,:));
-QD0=QD3P;
+rVoltABCV(3:3:end)=abs(VoltpABC(3,:));
+%三相电压相角
+rVoltABCA=zeros(busNum*3,1);
+rVoltABCA(1:3:end)=angle(VoltpABC(1,:));
+rVoltABCA(2:3:end)=angle(VoltpABC(2,:));
+rVoltABCA(3:3:end)=angle(VoltpABC(3,:));
+%三相负荷
+rPD3P=zeros(busNum*3,1);
+rPD3P(1:3:end)=phaseASpotLoadP;
+rPD3P(2:3:end)=phaseBSpotLoadP;
+rPD3P(3:3:end)=phaseCSpotLoadP;
+rQD3P=zeros(busNum*3,1);
+rQD3P(1:3:end)=phaseASpotLoadQ;
+rQD3P(2:3:end)=phaseBSpotLoadQ;
+rQD3P(3:3:end)=phaseCSpotLoadQ;
+rPD3P=rPD3P(Loadi);
+rQD3P=rQD3P(Loadi);
+%量测量
+sigma=0.03;
+mVoltABCV=rVoltABCV.*(1+normrnd(0,sigma,length(rVoltABCV),1));
+mPD3P=rPD3P.*(1+normrnd(0,sigma,length(rPD3P),1));
+mQD3P=rPD3P.*(1+normrnd(0,sigma,length(rQD3P),1));
+%量测方差
+wVolt=1./(abs(mVoltABCV*sigma).^2);
+wPD=1./(abs(mPD3P*sigma).^2);
+wQD=1./(abs(mQD3P*sigma).^2);
 %%
 RestraintCount=size(Loadi,1)*2; %约束条件数,放开所有QD
 Init_Z=sparse(ones(1,RestraintCount));
@@ -215,13 +242,13 @@ while(abs(Gap)>Precision)
 %     ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
     ddg=0;
     %%  开始构建deltF
-    deltF=func_deltF(wVolt,wPD,wQD,PD0,PD3P,QD3P,QD0,Volt,mVolt,Busnum,Loadi);
+    deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi);
 %     deltF=0;
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
     Mat_G=FormG(PD3P,QD3P,Loadi);
-    Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle);
+    Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
     Ly=Mat_H;
     Lz=FormLz(Mat_G,Init_L,Loadi);
     Lw=FormLw(Mat_G,Init_U,Loadi);

feeder13/data1.txt

@@ -4,7 +4,7 @@
 602	632	631	1800
 602	631	630	1700
 0
-632 30 0 0
+632 0 0 0
 0
 630 80 70 60 50 40 30
 631 40 10 11 50 30 30

func_deltdeltF.m

@@ -3,7 +3,7 @@ function deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount)
 %ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
 C=[wPD' wQD'];
 CVolt=wVolt;
-sizeC=size(C,2);
+sizeC=length(C);
 diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);
 deltdeltF=[
     diagC*2,sparse(sizeC,ContrlCount-sizeC);