如果虚拟电压方差太小，就会不容易收敛。要限制方差太小。

Signed-off-by: dugg@lab-desk <dugg@lab-desk>

IPMLoop.m

@@ -1,4 +1,4 @@
-function [ V1r,V1i,I1r,I1i ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref )
+function [ V1r,V1i,I1r,I1i,isConverged ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref )
 %把每个序的循环写在这个函数中。其实也就是内点法循环。
 V1r=1*ones(busNum,1);
 V1i=0*ones(busNum,1);
@@ -23,7 +23,7 @@ Init_U=1*sparse(ones(RestraintCount,1));
 Init_Y=sparse(2*busNum,1);%等式约束乘子
 Gap=(Init_L'*Init_Z-Init_U'*Init_W);
 kmax=60;
-while(abs(Gap)>0.000001)
+while(abs(Gap)>0.00001)
     if KK>kmax
         break;
     end
@@ -84,7 +84,11 @@ while(abs(Gap)>0.000001)
     fprintf('Gap %f\n',full(Gap));
     KK=KK+1;
 end
-f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
+% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
-
+if abs(Gap)<0.00001
+    isConverged=1;
+else
+    isConverged=0;
+end
 end
 

run.m

@@ -1,4 +1,4 @@
-function [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=run()
+function [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM,isConverged]=run()
 %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。
 %
 close all
@@ -146,9 +146,9 @@ mphaseASpotLoadQ=phaseASpotLoadQ.*(1+normrnd(0,sigma,length(phaseASpotLoadQ),1))
 mphaseBSpotLoadQ=phaseBSpotLoadQ.*(1+normrnd(0,sigma,length(phaseBSpotLoadQ),1));
 mphaseCSpotLoadQ=phaseCSpotLoadQ.*(1+normrnd(0,sigma,length(phaseCSpotLoadQ),1));
 %电压量测量
-mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1));
+% mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1));
-mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1));
+% mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1));
-mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1));
+% mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1));
 %
 mVoltpA=abs(VoltpA).*(1+normrnd(0,sigma,length(VoltpA),1));
 mVoltpB=abs(VoltpB).*(1+normrnd(0,sigma,length(VoltpB),1)).*exp(1j*-120/180*pi);
@@ -181,7 +181,6 @@ mIf2=-mIf2;
 %% 设定量测误差
 %% 先算正序的
 %平衡节点电流
-fprintf('ÕýÐò\n');
 BalI1r=real(-sum(rIf1));
 BalI1i=imag(-sum(rIf1));
 %电压
@@ -204,13 +203,13 @@ I1measurement=mIf1(Loadi);
 %正序测量误差的第二种形式
 wV1r=abs(real(V1measurement)).*sigma;
 wV1i=abs(imag(V1measurement)).*sigma;
+wV1i(abs(wV1i)<1e-5)=1e10;
 wI1r=abs(real(I1measurement)).*sigma;
 wI1i=abs(imag(I1measurement)).*sigma;
 % [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
 % fprintf('目标值 %f\n',full(f));
 %% 算负序的
-fprintf('¸ºÐò\n');
 BalI2r=real(-sum(rIf2));
 BalI2i=imag(-sum(rIf2));
 %电压
@@ -239,7 +238,6 @@ wI2i=abs(imag(I2measurement)).*sigma;
 % f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
 % fprintf('目标值 %f\n',full(f));
 %% 算零序
-fprintf('ÁãÐò\n');
 BalI0r=real(-sum(rIf0));
 BalI0i=imag(-sum(rIf0));
 %电压
@@ -269,13 +267,16 @@ tic
 for II=1:3
     
     if II==1
-        [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
+        fprintf('正序\n');
+        [ V1r,V1i,I1r,I1i,isConverged1 ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
     end
     if II==2
-        [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
+        fprintf('负序\n');
+        [ V2r,V2i,I2r,I2i,isConverged2 ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
     end
     if II==3
-        [ V0r,V0i,I0r,I0i ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 );
+        fprintf('零序\n');
+        [ V0r,V0i,I0r,I0i,isConverged0 ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 );
     end
     
 end
@@ -359,4 +360,6 @@ SEMeasurementM=(SEMeasurementM/mCount)^.5;
 %没有PD，QD的
 SEMeasurementNOPQM=sum(sum((( abs(rVoltABCV)-abs(mVoltABCV))./abs(mVoltABCV)./sigma).^2));
 SEMeasurementNOPQM=(SEMeasurementNOPQM/mCount)^.5;
+%返回收敛信息
+isConverged=isConverged1*isConverged2*isConverged0;
 end

test.m

@@ -8,7 +8,11 @@ SEMeasurementMSum=0;
 SEMeasurementNOPQMSum=0;
 N=1000;
 for I=1:N
-    [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=run();
+    [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM,isConverged]=run();
+    if isConverged==0
+        I=I-1;
+        continue
+    end
 % [JMeasurement]=run();
     JMeasurementSum=JMeasurement+JMeasurementSum;
     SEMeasurementSum=SEMeasurement+SEMeasurementSum;