准备改成标准的WLS方法。

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('.\..\DistributionNetwork-Power2Current\modified-feeder69\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,'.\..\DistributionNetwork-Power2Current\modified-feeder69\data.txt');
+    cap]=dataRead(lineZ,'./../DistributionNetwork-Power2Current/modified-feeder69\data-ijÏàû¸ººÉ.txt');
 % phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
 % phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
 % phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
@@ -14,6 +14,37 @@ lineZ=readLineZ('.\..\DistributionNetwork-Power2Current\modified-feeder69\linePa
 % phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002;
 % phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
 
+%¸ººÉ²»Æ½ºâ
+for I=1:length(phaseASpotLoadP)
+    roll=mod(round(rand()*10),3);
+    if roll==0
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
+        
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
+    end
+    if roll==1
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
+        
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
+    end
+    if roll==2
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
+        
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
+    end
+end
 
 
 %% 潮流计算begin
@@ -134,7 +165,6 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
 PGQG=CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
 fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
 %% 潮流计算end
-
 busNum=length(nodeNum);
 Busnum=busNum;
 % PQi=setxor(nodeNum,Balance);
@@ -215,7 +245,8 @@ rQD3P=rQD3P(Loadi);
 noLoadi=[1,8,10,12];
 % noLoadi=Loadi;
 % noLoadi=[1,11];
-noLoadi=[];
+noLoadi=[5];
+% noLoadi=[5,11,19,27,37,43,54,61];
 noPQi3P=zeros(length(noLoadi)*3,1);
 noPQi3P(1:3:end)=(noLoadi-1)*3+1;
 noPQi3P(2:3:end)=(noLoadi-1)*3+2;
@@ -223,23 +254,23 @@ noPQi3P(3:3:end)=(noLoadi-1)*3+3;
 % noPQi3P=Loadi;
 %量测量
 sigma=0.03;
-VoltSigma=(1+normrnd(0,sigma,length(rVoltABCV),1));
+VoltSigma=(1+normrnd(0,sigma/3,length(rVoltABCV),1));
 mVoltABCV=rVoltABCV.*VoltSigma;
+% mVoltABCV(noPQi3P)=rVoltABCV(noPQi3P).*(1+normrnd(0,0.10,length(noPQi3P),1));
 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,noPQi3P))=mPD3P(ismember(Loadi,noPQi3P)).*(1+normrnd(0,0.15,length(noPQi3P),1));
-mQD3P(ismember(Loadi,noPQi3P))=mQD3P(ismember(Loadi,noPQi3P)).*(1+normrnd(0,0.15,length(noPQi3P),1));
+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).^2);
-wPD=1./(abs(mPD3P*.15).^2);
-wQD=1./(abs(mQD3P*.15).^2);
-
+wVolt=1./(abs(mVoltABCV*sigma/3).^2);
+wPD=1./(abs(mPD3P*sigma).^2);
+wQD=1./(abs(mQD3P*sigma).^2);
 wVolt(setdiff(1:length(wVolt),Loadi))=0;%只有负荷处才有电压量测。
 wVolt( noPQi3P)=0;
-wPD(ismember( Loadi,noPQi3P))=0;
-wQD(ismember( Loadi,noPQi3P))=0;
+wPD(ismember( Loadi,noPQi3P))=1./(abs(mPD3P(ismember( Loadi,noPQi3P))*0.15).^2);
+wQD(ismember( Loadi,noPQi3P))=1./(abs(mQD3P(ismember( Loadi,noPQi3P))*0.15).^2);
 %%
 % RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数
 % RestraintCount=size(Loadi,1)*2; %约束条件数
@@ -251,10 +282,12 @@ kmax=200;
 Precision=1e-5;
 CenterA=0.1;
 %% 加误差
-Volt=1*Vp3m;
-UAngel=Vp3a;
+
+
 maxD=100;
 tic
+Volt=1*Vp3m;
+UAngel=Vp3a;
 while(maxD>Precision)
     if KK>kmax
         break;
@@ -288,14 +321,14 @@ JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD
 Busnum=busNum;
 mCount=Busnum*3+length(Loadi)*3*2;
 %估计量质量
-AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(Volt))));
-AME_VAngle=sum(sum(abs( rVoltABCA-UAngel)));
-AME_PD=sum(sum(abs(rPD3P-PD3P)));
-AME_QD=sum(sum(abs(rQD3P-QD3P)));
+AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(Volt))))/3/busNum;
+AME_VAngle=sum(sum(abs( rVoltABCA-UAngel)))/3/busNum;
+AME_PD=sum(sum(abs(rPD3P-PD3P)))/length(Loadi);
+AME_QD=sum(sum(abs(rQD3P-QD3P)))/length(Loadi);
 %计算与量测值的
-AME_mVolt=sum(sum(abs( mVoltABCV-rVoltABCV)));
-AME_mPD=sum(sum(abs(rPD3P-mPD3P)));
-AME_mQD=sum(sum(abs(rQD3P-mQD3P)));
+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);
 %返回收敛信息
 isConverged=1;
 if KK>=kmax

SolveIt.m

@@ -8,8 +8,7 @@ function XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,m
 %     ];
 aa=[
     deltF'*diag([wPD ;wQD ;wVolt])*deltF deltH'
-    deltH ones(size(deltH,1),size(deltH',2))
-    
+    deltH zeros(size(deltH,1),size(deltH',2))
 ];
 % yy=[Lx;-Ly];
 % t3=2*wPD.*(PD-PD0);