From 62389acec5a4eb802f0d5487ca8d1a49434b440f Mon Sep 17 00:00:00 2001
From: facat <facat@facat.cn>
Date: Fri, 30 Aug 2013 16:11:16 +0800
Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E6=94=B9=E4=BA=86=E6=9D=83=E9=87=8D?=
 =?UTF-8?q?=E4=BC=B0=E8=AE=A1=E7=BB=93=E6=9E=9C=E5=BE=88=E5=A5=BD=EF=BC=8C?=
 =?UTF-8?q?=E4=BD=86=E6=98=AF=E5=AF=B9=E4=BA=8E1047=E8=8A=82=E7=82=B9?=
 =?UTF-8?q?=EF=BC=8C=E6=94=B6=E6=95=9B=E6=80=A7=E4=B8=8D=E5=A5=BD=E3=80=82?=
 =?UTF-8?q?=E5=87=86=E5=A4=87=E7=94=A8LM=E6=B3=95=E3=80=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Signed-off-by: facat <facat@facat.cn>
---
 plotAndComparison.m |  2 ++
 run.m               | 76 +++++++++++++++++++++++----------------------
 2 files changed, 41 insertions(+), 37 deletions(-)

diff --git a/plotAndComparison.m b/plotAndComparison.m
index 7eda615..7c870b9 100644
--- a/plotAndComparison.m
+++ b/plotAndComparison.m
@@ -4,9 +4,11 @@ subplot(2,1,1);
 plot(x,rVolt,'r-');
 hold on
 plot(x,SEVolt);
+legend('real','se');
 subplot(2,1,2);
 plot(x,rVAngel,'r-');
 hold on
 plot(x,SEVAngel);
+legend('real','se');
 end
 
diff --git a/run.m b/run.m
index 2f9ea8c..fa012de 100644
--- a/run.m
+++ b/run.m
@@ -4,7 +4,7 @@ clc
 addpath('.\Powerflow')
 [~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('s1047.dat', '0');
 %% 开始生成量测量
-sigma=0.05;% 标准差
+sigma=0.03;% 标准差
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -70,9 +70,9 @@ PDQDi=union(PDi,QDi);
 onlyPG=setdiff(PGi,PDQDi);
 onlyQG=setdiff(QGi,PDQDi);
 %% 计算方差
-measureSigma=([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma).^2;
-measureSigma(measureSigma<1e-4)=mean(measureSigma);
-W=diag(1./measureSigma) ;
+measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
+measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
+W=diag(1./measureSigma.^2) ;
 % W=sparse(1:length(W),1:length(W),400,length(W),length(W));
 %% 冗余度计算
 stateVarCount=2*length(Volt);
@@ -107,7 +107,7 @@ optimalCondition=100;
 eps=1e-4;
 % while max(abs(g))>1e-5;
 % while maxD>1e-5
-while max(abs(optimalCondition))>eps
+while max(abs(maxD))>eps
     % 电压
     dV_dV=sparse(1:length(mVolt),1:length(mVolt),1,length(mVolt),length(mVolt));%电压量测量的微分
     dV_dTyta=sparse(length(mVolt),length(mVolt));
@@ -200,33 +200,33 @@ while max(abs(optimalCondition))>eps
         ,length(transI),length(mVolt));%变压器
     %% 考虑等式约束
     % 等式约束的Jacobi
-    r=newwordParameter.r;
-    c=newwordParameter.c;
-    Yangle=newwordParameter.Yangle;
-    VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ;
-    YdotSin=Y.* ( spfun(@sin,VAngleIJ) );
-    YdotCos=Y.* ( spfun (@cos, VAngleIJ ) );
-    diag_Volt_YdotCos=diag(SEVolt)*YdotCos;
-    diag_Volt_YdotSin=diag(SEVolt)*YdotSin;
-    YdotCosVolt=YdotCos*Volt;
-    YdotSinVolt=YdotSin*Volt;
-    diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt;
-    diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt;
-    diag_YdotSinVolt_=diag(YdotSinVolt);
-    diag_YdotCosVolt_=diag(YdotCosVolt);
-    dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次
-    dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta
-    dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
-    dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
-    % C 是等式约束 c 的Jacobi
-    C=[dPdV dPdTyta;
-        dQdV dQdTyta];
-    C=C(zerosInjectionIndex,:);
-    % 形成等式约束 c
-    nodeP=diag_Volt_YdotCosVolt;
-    nodeQ=diag_Volt_YdotSinVolt;
-    nodePQ=[nodeP;nodeQ];
-    c=nodePQ(zerosInjectionIndex);
+%     r=newwordParameter.r;
+%     c=newwordParameter.c;
+%     Yangle=newwordParameter.Yangle;
+%     VAngleIJ=sparse(r,c,SEVAngle(r)-SEVAngle(c) -Yangle,length(mVolt),length(mVolt)) ;
+%     YdotSin=Y.* ( spfun(@sin,VAngleIJ) );
+%     YdotCos=Y.* ( spfun (@cos, VAngleIJ ) );
+%     diag_Volt_YdotCos=diag(SEVolt)*YdotCos;
+%     diag_Volt_YdotSin=diag(SEVolt)*YdotSin;
+%     YdotCosVolt=YdotCos*Volt;
+%     YdotSinVolt=YdotSin*Volt;
+%     diag_Volt_YdotCosVolt=diag_Volt_YdotCos*Volt;
+%     diag_Volt_YdotSinVolt=diag_Volt_YdotSin*Volt;
+%     diag_YdotSinVolt_=diag(YdotSinVolt);
+%     diag_YdotCosVolt_=diag(YdotCosVolt);
+%     dPdTyta=diag_Volt_YdotSin*diag(SEVolt)-diag_YdotSinVolt_*diag(SEVolt); % 简化第三次
+%     dQdTyta=-diag_Volt_YdotCos*diag(SEVolt)+diag_YdotCosVolt_*diag(SEVolt);%dQ/dThyta
+%     dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
+%     dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
+%     % C 是等式约束 c 的Jacobi
+%     C=[dPdV dPdTyta;
+%         dQdV dQdTyta];
+%     C=C(zerosInjectionIndex,:);
+%     % 形成等式约束 c
+%     nodeP=diag_Volt_YdotCosVolt;
+%     nodeQ=diag_Volt_YdotSinVolt;
+%     nodePQ=[nodeP;nodeQ];
+%     c=nodePQ(zerosInjectionIndex);
     %% 进入迭代
     H=[dV_dV,dV_dTyta;
         dLPij_dVi+dLPij_dVj,dLPij_dThetai+dLPij_dThetaj ;
@@ -245,9 +245,12 @@ while max(abs(optimalCondition))>eps
     g=-H'*W*(z-h);
    
     % 形成大的求解矩阵
-    a=[G C';
-        C zeros(size(C,1),size(C,1))];
-    b=[-g;-c];
+%     a=[G C';
+%         C zeros(size(C,1),size(C,1))];
+%     b=[-g;-c];
+    a=G;
+    b=-g;
+
     % 平衡节点相角恒定
     a(length(mVolt)+Balance,:)=0;
     a(:,length(mVolt)+Balance)=0;
@@ -268,8 +271,7 @@ while max(abs(optimalCondition))>eps
     Iteration=Iteration+1;
     SEVAngle=SEVAngle+dX(length(mVolt)+1:length(mVolt)*2)*dXStep;
     lamda=-dX(length(mVolt)*2+1:end);
-    optimalCondition=[-g+C'*lamda;
-        c];
+    optimalCondition=-g;
     optimalCondition(Balance)=0;
     optimalCondition(Balance+length(mVolt))=0;
 end