From 857f181c471255fd96888e29d5699cd39aa96723 Mon Sep 17 00:00:00 2001
From: facat <facat@facat.cn>
Date: Sun, 18 Aug 2013 12:14:29 +0800
Subject: [PATCH] =?UTF-8?q?1.=E5=A2=9E=E5=8A=A0=E6=AD=A3=E6=80=81=E5=88=86?=
 =?UTF-8?q?=E5=B8=83=E4=BB=A3=E7=A0=81=E6=A8=A1=E6=8B=9F=E7=94=A8=202.?=
 =?UTF-8?q?=E4=BF=AE=E6=94=B9=E4=BA=86=E6=9D=83=E9=87=8D=E7=9A=84=E8=AE=A1?=
 =?UTF-8?q?=E7=AE=97=EF=BC=8C=E9=80=82=E5=BD=93=E5=A4=84=E7=90=86=E9=87=8F?=
 =?UTF-8?q?=E6=B5=8B=E5=80=BC=E5=A4=AA=E5=B0=8F=E7=9A=84=E6=83=85=E5=86=B5?=
 =?UTF-8?q?=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>
---
 LevenbergMarquardt.m |  9 +++++++++
 NormalDistribution.m |  8 ++++++++
 run.m                | 18 ++++++++++++++----
 3 files changed, 31 insertions(+), 4 deletions(-)
 create mode 100644 LevenbergMarquardt.m
 create mode 100644 NormalDistribution.m

diff --git a/LevenbergMarquardt.m b/LevenbergMarquardt.m
new file mode 100644
index 0000000..46fe336
--- /dev/null
+++ b/LevenbergMarquardt.m
@@ -0,0 +1,9 @@
+function [ output_args ] = LevenbergMarquardt( J,x )
+t=0.1;
+A=J'*J;
+u=t*max(diag(A));
+dh=(A+u*eye(size(A)) )\g;
+
+
+end
+
diff --git a/NormalDistribution.m b/NormalDistribution.m
new file mode 100644
index 0000000..186f461
--- /dev/null
+++ b/NormalDistribution.m
@@ -0,0 +1,8 @@
+sim=100000;
+%add=0;
+s=zeros(sim,1);
+for i=1:sim
+    s(i)=6*(1+normrnd(0,0.2));
+end
+mu=sum(s)/sim %平均值
+sum((s-mu).^2/sim)^.5 %标准差
diff --git a/run.m b/run.m
index 21a8167..389e569 100644
--- a/run.m
+++ b/run.m
@@ -2,9 +2,9 @@ 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('ieee118.dat', '0');
 %% 开始生成量测量
-sigma=0.09;% 标准差
+sigma=0.05;% 标准差
 %% 电压
 %电压幅值
 rVolt=Volt; %幅值
@@ -69,8 +69,10 @@ PDQDi=union(PDi,QDi);
 onlyPG=setdiff(PGi,PDQDi);
 onlyQG=setdiff(QGi,PDQDi);
 %% 计算方差
-W=(diag([rVolt;rBranchP;rBranchQ;rTransP;rTransQ])*sigma).^2;
-W=sparse(1:length(W),1:length(W),400,length(W),length(W));
+measureSigma=([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma).^2;
+measureSigma(measureSigma<1e-4)=mean(measureSigma);
+W=diag(1./measureSigma) ;
+% W=sparse(1:length(W),1:length(W),400,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);
@@ -96,6 +98,7 @@ fprintf('
 %% 自己写的微分代码
 % 以下都是Jacobi矩阵
 SEVolt=sparse(ones(length(mVolt),1));
+SEVolt(Balance)=rVolt(Balance);
 SEVAngle=sparse(zeros(length(mVolt),1));
 maxD=1000;
 Iteration=0;
@@ -213,10 +216,17 @@ while max(abs(g))>1e-5;
     G(:,length(mVolt)+Balance)=0;
     G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2);
     g(length(mVolt)+Balance)=0;
+    % 平衡节点电压恒定;
+    G(Balance,:)=0;
+    G(:,Balance)=0;
+    G=G+sparse(Balance,Balance,1,length(mVolt)*2,length(mVolt)*2);
+    g(Balance)=0;
+    % 求解修正量
     dX=G\-g;
     dXStep=1;
 %     dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
     maxD=max(abs(dX))
+    fprintf('max abs g:%f\n',full(max(abs(g))));
     % 更新变量
     SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
     Iteration=Iteration+1;