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;