1.加了电流幅值量测公式

2.在状态估计中加了电流量测

LineCurrent.m

@@ -0,0 +1,7 @@
+function [ output_args ] = LineCurrent( Linei,Linej,Liner,Linex,Volt,VAngle )
+cmpV=Volt.*exp(1j*VAngle);
+GB=1./(Liner+1j*Linex);
+cmpI=( cmpV(Linei)-cmpV(Linej) ).*GB;
+output_args=abs(cmpI);
+end
+

run.m

@@ -2,7 +2,7 @@ clear
 clc
 % yalmip('clear')
 addpath('.\Powerflow')
-[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4-DN.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\算例\feeder33\feeder33ieee.txt', '0');
 % 'E:\算例\feeder33\feeder33ieee.txt'
 %% 开始生成量测量
 sigma=0.03;% 标准差
@@ -28,6 +28,8 @@ lineX=newwordParameter.line.lineX;
 lineB2=newwordParameter.line.lineB2;
 lineG=real(1./(lineR+1j*lineX));
 lineB=imag(1./(lineR+1j*lineX));
+lineY=abs(1./(lineR+1j*lineX));
+lineYa=angle(1./(lineR+1j*lineX));
 cmpBranchI=BranchI( cmpV,lineI,lineJ,lineR,lineX );%复数支路电流
 rBranchI=abs(cmpBranchI);% 支路电流幅值
 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量
@@ -72,7 +74,7 @@ onlyPG=setdiff(PGi,PDQDi);
 onlyQG=setdiff(QGi,PDQDi);
 %% 计算方差
 % measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
-measureSigma=abs(([rVolt;rPD(PDi);rQD(QDi);].*sigma));
+measureSigma=abs(([rVolt;rBranchI;rPD(PDi);rQD(QDi);].*sigma));
 measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
 W=sparse(diag(1./measureSigma.^2)) ;
 % W=eye(length(W));
@@ -115,8 +117,8 @@ v=2;
 ojbFunDecrease=1000;% 目标函数下降
 % 以下都是Jacobi矩阵
 % while max(abs(g))>1e-5;
-% while maxD>1e-5
+while maxD>1e-5
-while max(abs(optimalCondition))>eps
+% while max(abs(optimalCondition))>eps
     % 电压
     dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
     dV_dTyta=sparse(length(mVolt),length(mVolt));
@@ -207,6 +209,32 @@ while max(abs(optimalCondition))>eps
         transG.*cos(SEVAngle(transI)-SEVAngle(transJ)) +transB.*sin(SEVAngle(transI)-SEVAngle(transJ))...
         ) ...
         ,length(transI),length(mVolt));%变压器
+    %电流幅值实部
+    dLIrij_dVi=sparse(1:length(lineI),lineI, ...
+        lineY.*cos(SEVAngle(lineI)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIrij_dVj=sparse(1:length(lineI),lineI, ...
+        -lineY.*cos(SEVAngle(lineJ)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIrij_dThetai=sparse(1:length(lineI),lineI, ...
+        -SEVolt(lineI).*lineY.*sin(SEVAngle(lineI)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIrij_dThetaj=sparse(1:length(lineI),lineI, ...
+        SEVolt(lineJ).*lineY.*sin(SEVAngle(lineJ)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    %电流幅值虚部
+    dLIiij_dVi=sparse(1:length(lineI),lineI, ...
+        lineY.*sin(SEVAngle(lineI)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIiij_dVj=sparse(1:length(lineI),lineI, ...
+        -lineY.*sin(SEVAngle(lineJ)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIiij_dThetai=sparse(1:length(lineI),lineI, ...
+        SEVolt(lineI).*lineY.*cos(SEVAngle(lineI)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
+    dLIiij_dThetaj=sparse(1:length(lineI),lineI, ...
+        -SEVolt(lineJ).*lineY.*cos(SEVAngle(lineJ)+lineYa)...
+        ,length(lineI),length(mVolt));%线路的
     %% 考虑注入功率
     % 等式约束的Jacobi
     r=newwordParameter.r;
@@ -289,7 +317,12 @@ while max(abs(optimalCondition))>eps
 %         dLQij_dVi+dLQij_dVj,dLQij_dThetai+dLQij_dThetaj ;
 %         dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
 %         dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
+%     H=[dV_dV,dV_dTyta;
+%         dPdV(PDi,:),dPdTyta(PDi,:);
+%         dQdV(QDi,:),dQdTyta(QDi,:)];%jacobi
+    
     H=[dV_dV,dV_dTyta;
+        2*dLIrij_dVi+2*dLIiij_dVi+2*dLIrij_dVj+2*dLIiij_dVj,2*dLIrij_dThetai+2*dLIiij_dThetai+2*dLIrij_dThetaj+2*dLIiij_dThetaj
         dPdV(PDi,:),dPdTyta(PDi,:);
         dQdV(QDi,:),dQdTyta(QDi,:)];%jacobi
     
@@ -303,12 +336,12 @@ while max(abs(optimalCondition))>eps
     SEPD=diag(SEVolt)*Y.* ( spfun (@cos, VAngleIJ ) )*SEVolt;
     SEQD=diag(SEVolt)*Y.* ( spfun(@sin,VAngleIJ) )*SEVolt;
     
-    h=[SEVolt;SEPD(PDi);SEQD(QDi);];
+    h=[SEVolt;abs(SEBranchI).^2;SEPD(PDi);SEQD(QDi);];
     
 %     h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
 %     z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
 
-    z=[mVolt;-mPD(PDi);-mQD(QDi)];
+    z=[mVolt;mBranchI.^2;-mPD(PDi);-mQD(QDi)];
 
 
     G=H'*W*H;
@@ -384,11 +417,20 @@ fprintf('Ŀ
 fprintf('相对误差\n');
 (abs(rVolt-double(SEVolt)))./(rVolt);
 MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle)
-plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
+% plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
 % 检查最优性条件
 
 if any(abs(optimalCondition)>eps)
     fprintf('最优性条件不满足\n')
 else
     fprintf('最优性条件满足\n')
-end
+end
+
+
+Linei=newwordParameter.line.lineI;
+Linej=newwordParameter.line.lineJ;
+Liner=newwordParameter.line.lineR;
+Linex=newwordParameter.line.lineX;
+SELineCurrent=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngle );
+RealLineCurrent=LineCurrent( Linei,Linej,Liner,Linex,rVolt,rVAngel );
+sqrt(sum((SELineCurrent-RealLineCurrent).^2)/length(SELineCurrent))