修正一个小错误，原来没有用潮流计算的值作为头节点的值就进行状态估计了。

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

CalPGQG.m

@@ -0,0 +1,25 @@
+function [ output_args ] = CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ )
+%%利用Fortiscue方法计算得到了相电压和相角后，反推负荷功率，检查是否和给定的功率一致。
+busNum=size(phaseABCY,1)/3;
+VoltpA=VoltpABC(1,:);
+VoltpB=VoltpABC(2,:);
+VoltpC=VoltpABC(3,:);
+Vp3=sparse(busNum*3,1);
+Vp3(1:3:end)=VoltpA;
+Vp3(2:3:end)=VoltpB;
+Vp3(3:3:end)=VoltpC;
+Ip3=phaseABCY*Vp3;
+Sp3=Vp3.*conj(Ip3);
+%VoltpABC.*conj(IpABC);
+pLoadABC=sparse(length(phaseASpotLoadP)*3,1);
+pLoadABC(1:3:end)=phaseASpotLoadP+1j*phaseASpotLoadQ;
+pLoadABC(2:3:end)=phaseBSpotLoadP+1j*phaseBSpotLoadQ;
+pLoadABC(3:3:end)=phaseCSpotLoadP+1j*phaseCSpotLoadQ;
+ck=Sp3+pLoadABC;
+% ck(3*(Balance-1)+1:3*(Balance-1)+3)=0;
+% if any(abs(ck)>1e-5)
+%     fprintf('反推回的功率不匹配。潮流方程约束不满足。\n');
+% end
+output_args=ck;
+end
+

FormLw.m

@@ -1,4 +1,4 @@
-function Lw=FormLw(Mat_G,Init_U,Loadi)
+function Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P)
 % KK=999;
 % %PU=GenU(:,2);%发电机有功上界
 % 
@@ -24,8 +24,10 @@ function Lw=FormLw(Mat_G,Init_U,Loadi)
 % QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end));
 % PF=0.85;
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
-PDU=1.8*ones(length(Loadi),1);
+% PDU=1.8*ones(length(Loadi),1);
-QDU=1.8*ones(length(Loadi),1);
+% QDU=1.8*ones(length(Loadi),1);
+PDU=1.05*rPD3P;
+QDU=1.05*rQD3P;
 t1=([PDU',QDU'])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLz.m

@@ -1,4 +1,4 @@
-function Lz=FormLz(Mat_G,Init_L,Loadi)
+function Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P)
 % KK=999;
 % 
 % VoltL=(0.9)*ones(1,Busnum);
@@ -21,8 +21,10 @@ function Lz=FormLz(Mat_G,Init_L,Loadi)
 % QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end));
 % QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end));
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
-PDL=0*ones(length(Loadi),1);
+% PDL=0*ones(length(Loadi),1);
-QDL=0*ones(length(Loadi),1);
+% QDL=0*ones(length(Loadi),1);
+PDL=0.95*rPD3P;
+QDL=0.95*rQD3P;
 t1=([PDL',QDL'])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

OPF.m

@@ -120,6 +120,7 @@ disp([setIJ,nodeNum ])
 ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
     phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
+PGQG=CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
 fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
 %% 潮流计算end
 
@@ -141,32 +142,36 @@ 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;
+PD3P(1:3:end)=phaseASpotLoadP*0.5;
-PD3P(2:3:end)=phaseBSpotLoadP;
+PD3P(2:3:end)=phaseBSpotLoadP*0.5;
-PD3P(3:3:end)=phaseCSpotLoadP;
+PD3P(3:3:end)=phaseCSpotLoadP*0.5;
-QD3P(1:3:end)=phaseASpotLoadQ;
+QD3P(1:3:end)=phaseASpotLoadQ*0.5;
-QD3P(2:3:end)=phaseBSpotLoadQ;
+QD3P(2:3:end)=phaseBSpotLoadQ*0.5;
-QD3P(3:3:end)=phaseCSpotLoadQ;
+QD3P(3:3:end)=phaseCSpotLoadQ*0.5;
 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;
 QGi3P(3:3:end)=(QGi-1)*3+3;
-PGA=sum(PD3P(1:3:end))*1.03;
-PGB=sum(PD3P(2:3:end))*1.03;
-PGC=sum(PD3P(3:3:end))*1.03;
-QGA=sum(QD3P(1:3:end))*1.03;
-QGB=sum(QD3P(2:3:end))*1.03;
-QGC=sum(QD3P(3:3:end))*1.03;
-PG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[PGA,PGB,PGC],Busnum*3,1);
-QG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[QGA,QGB,QGC],Busnum*3,1);
 Vp3m=abs(Vp3);
 Vp3a=angle(Vp3);
 Balance3P=zeros(length(Balance)*3,1);
 Balance3P(1:3:end)=(Balance-1)*3+1;
 Balance3P(2:3:end)=(Balance-1)*3+2;
 Balance3P(3:3:end)=(Balance-1)*3+3;
+
+% PGA=sum(PD3P(1:3:end))*1.03;
+% PGB=sum(PD3P(2:3:end))*1.03;
+% PGC=sum(PD3P(3:3:end))*1.03;
+% QGA=sum(QD3P(1:3:end))*1.03;
+% QGB=sum(QD3P(2:3:end))*1.03;
+% QGC=sum(QD3P(3:3:end))*1.03;
+% PG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[PGA,PGB,PGC],Busnum*3,1);
+% QG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[QGA,QGB,QGC],Busnum*3,1);
+PG3P=real(PGQG);
+QG3P=imag(PGQG);
+
 Vp3a((Balance-1)*3+1)=0;
 Vp3a((Balance-1)*3+2)=-120/180*pi;
 Vp3a((Balance-1)*3+3)=+120/180*pi;
@@ -195,13 +200,19 @@ rPD3P=rPD3P(Loadi);
 rQD3P=rQD3P(Loadi);
 %量测量
 sigma=0.03;
-mVoltABCV=rVoltABCV.*(1+normrnd(0,sigma,length(rVoltABCV),1));
+VoltSigma=(1+normrnd(0,sigma,length(rVoltABCV),1));
-mPD3P=rPD3P.*(1+normrnd(0,sigma,length(rPD3P),1));
+mVoltABCV=rVoltABCV.*VoltSigma;
-mQD3P=rPD3P.*(1+normrnd(0,sigma,length(rQD3P),1));
+PD3PSigma=(1+normrnd(0,sigma,length(rPD3P),1));
+mPD3P=rPD3P.*PD3PSigma;
+QD3PSigma=(1+normrnd(0,sigma,length(rQD3P),1));
+mQD3P=rPD3P.*QD3PSigma;
 %量测方差
-wVolt=1./(abs(mVoltABCV*sigma).^2);
+wVolt=1./(abs(mVoltABCV*sigma).^2)*0;
-wPD=1./(abs(mPD3P*sigma).^2);
+wPD=1./(abs(mPD3P*sigma).^2)*1;
-wQD=1./(abs(mQD3P*sigma).^2);
+wQD=1./(abs(mQD3P*sigma).^2)*1;
+% wVolt=1./abs(VoltSigma).^2;
+% wPD=1./abs(PD3PSigma).^2;
+% wQD=1./abs(QD3PSigma).^2;
 %%
 RestraintCount=size(Loadi,1)*2; %约束条件数,放开所有QD
 Init_Z=sparse(ones(1,RestraintCount));
@@ -250,10 +261,10 @@ while(abs(Gap)>Precision)
     Mat_G=FormG(PD3P,QD3P,Loadi);
     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);
+    Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P);
-    Lw=FormLw(Mat_G,Init_U,Loadi);
+    Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
+%     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
 %     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,Busnum,Loadi);
@@ -264,6 +275,7 @@ while(abs(Gap)>Precision)
     fprintf('%f\n',full(Gap));
     KK=KK+1;
 end
+(rVoltABCV-Volt)./rVoltABCV*100
 fprintf('迭代次数%d\n',KK);
 toc