提交到gogs

Signed-off-by: dmy@lab <dmy@lab.lab>

StatDeviation.m

@@ -0,0 +1,13 @@
+function StatDeviation(rVolt,SEVolt,rVAngel,SEVAngel,sigma)
+t1=[rVolt;rVAngel];
+t2=[SEVolt;SEVAngel];
+t3=sum((t2-t1).^2)/length(t1);
+t3=t3.^0.5;
+fprintf('״̬Á¿Í³¼ÆÎó²î %f\n',full(t3));
+t1=[rVolt];
+t2=[SEVolt];
+t3=sum( ( (t2-t1)./(rVolt.*sigma) ).^2 )/length(t1);
+t3=t3.^0.5;
+fprintf('Á¿²âÁ¿Í³¼ÆÎó²î %f\n',full(t3));
+end
+

plotAndComparison.m

@@ -1,14 +1,18 @@
-function plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel )
+function plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel,mVolt )
 x=1:length(rVolt);
 subplot(2,1,1);
 plot(x,rVolt,'r-');
+title('µçѹ')
 hold on
 plot(x,SEVolt);
-legend('real','se');
+hold on
+plot(x,mVolt,'g');
+legend('actual','estimate','measure');
 subplot(2,1,2);
 plot(x,rVAngel,'r-');
 hold on
 plot(x,SEVAngel);
-legend('real','se');
+title('Ïà½Ç')
+legend('actual','estimate');
 end
 

run.m

@@ -2,13 +2,13 @@ clear
 clc
 % yalmip('clear')
 addpath('.\Powerflow')
-[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('s1047.dat', '0');
+[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\ËãÀý\lipuy2.dat', '0');
 %% 节点总数
 busNum=length(Volt);
 %% 变量总数
 varNum=busNum*2;
 %% 开始生成量测量
-sigma=0.03;% ±ê×¼²î
+sigma=0.0000001;% ±ê×¼²î
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -81,7 +81,7 @@ onlyQG=setdiff(QGi,PDQDi);
 measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
 measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
 W=sparse(diag(1./measureSigma.^2)) ;
-% W=sparse(1:length(W),1:length(W),1,length(W),length(W));
+W=sparse(1:length(W),1:length(W),1,length(W),length(W));
 %% 冗余度计算
 stateVarCount=2*length(Volt);
 measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
@@ -107,23 +107,23 @@ fprintf('
 %% 自己写的微分代码
 % 初始化一些值
 %SEVolt=sparse(ones(length(mVolt),1));
-SEVolt=rVolt;
+SEVolt=mVolt;
 % SEVolt=rVolt;
 % SEVolt(1:2:end)=1.2;
 % load('SEVolt');
 % load('SEVAngle');
 SEVolt(Balance)=rVolt(Balance);
-% SEVAngle=-0.1*sparse(ones(length(mVolt),1));
+SEVAngle=-0.3*sparse(ones(length(mVolt),1));
-% SEVAngle(Balance)=0;
+SEVAngle(Balance)=0;
-SEVAngle=rVAngel;
+% SEVAngle=mVAngel;
 % SEVAngle=rVAngel;
 % SEVolt2=    load('SEVolt2');
 % SEVolt2=SEVolt2.SEVolt;
 % SEVAngle2=    load('SEVAngle2');
 % SEVAngle2=SEVAngle2.SEVAngle;
-maxD=1000;
+maxD=1;
 Iteration=0;
-optimalCondition=100;
+optimalCondition=1000;
 eps=1e-4;
 mu=0;
 v=2;
@@ -131,7 +131,8 @@ ojbFunDecrease=1000;% Ŀ
 % 以下都是Jacobi矩阵
 % while max(abs(g))>1e-5;
 % while maxD>1e-5
-while max(abs(optimalCondition))>eps
+%while maxD>1e-5 && max(abs(optimalCondition))/maxD>1e2
+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));
@@ -444,20 +445,28 @@ while max(abs(optimalCondition))>eps
     %海森矩阵对角，即
     % [ d2h_dV2     0
     %      0     d2h_dT2]
-    hessian14=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2 ...
+%     hessian14=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2 ...
-        +d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ...
+%         +d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ...
-        +d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ...
+%         +d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ...
-        +d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ...;
+%         +d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ...;
-        +d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj'+d2TPij_dVj2 ...
+%         +d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj'+d2TPij_dVj2 ...
-        +d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTidTj'+d2TPij_dTj2 ...
+%         +d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTidTj'+d2TPij_dTj2 ...
-        +d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ...
+%         +d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ...
-        +d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2;
+%         +d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2;
+hessian14=d2LPij_dVi2+d2LPij_dVidVj+d2LPij_dVidVj'+d2LPij_dVj2;
+hessian14=hessian14+d2LPij_dTi2+d2LPij_dTidTj+d2LPij_dTidTj'+d2LPij_dTj2 ;
+hessian14=hessian14+d2LQij_dVi2+d2LQij_dVidVj+d2LQij_dVidVj'+d2LQij_dVj2 ;
+hessian14=hessian14+d2LQij_dTi2+d2LQij_dTidTj+d2LQij_dTidTj'+d2LQij_dTj2 ;
+hessian14=hessian14+d2TPij_dVi2+d2TPij_dVidVj+d2TPij_dVidVj'+d2TPij_dVj2 ;
+hessian14=hessian14+d2TPij_dTi2+d2TPij_dTidTj+d2TPij_dTidTj'+d2TPij_dTj2 ;
+hessian14=hessian14+d2TQij_dVi2+d2TQij_dVidVj+d2TQij_dVidVj'+d2TQij_dVj2 ;
+hessian14=hessian14+d2TQij_dTi2+d2TQij_dTidTj+d2TQij_dTidTj'+d2TQij_dTj2;
     hessian2=d2LPij_dVidTi+d2LPij_dVidTj+d2LPij_dVjdTi+d2LPij_dVjdTj ...
         +d2LQij_dVidTi+d2LQij_dVidTj+d2LQij_dVjdTi+d2LQij_dVjdTj ...
         +d2TPij_dVidTi+d2TPij_dVidTj+d2TPij_dVjdTi+d2TPij_dVjdTj ...
         +d2TQij_dVidTi+d2TQij_dVidTj+d2TQij_dVjdTi+d2TQij_dVjdTj;
     hessian3=hessian2';
-    if Iteration<=-1
+    if Iteration<=4
          G=H'*W*H;
     else    
         G=H'*W*H+hessian14+hessian2+hessian3;
@@ -496,9 +505,9 @@ fprintf('Ŀ
 fprintf('相对误差\n');
 (abs(rVolt-double(SEVolt)))./(rVolt);
 MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle)
-plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
+plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle,mVolt )
 % 检查最优性条件
-
+StatDeviation(rVolt,SEVolt,rVAngel,SEVAngle,sigma);
 if any(abs(optimalCondition)>eps)
     fprintf('最优性条件不满足\n')
 else