facat
/
stateestimation-ipm-forcomparison
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

1.增加正态分布代码模拟用 

2.修改了权重的计算，适当处理量测值太小的情况。

Signed-off-by: facat <facat@facat.cn>
This commit is contained in:
facat 2013-08-18 12:14:29 +08:00
parent ea40031257
commit 857f181c47
3 changed files with 31 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
LevenbergMarquardt.m Normal file
View File

@ -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

8
NormalDistribution.m Normal file
View File

@ -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 %±ê×¼²î

18
run.m
View File

@ -2,9 +2,9 @@ clear
clc clc
% yalmip('clear') % yalmip('clear')
addpath('.\Powerflow') 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; %幅值 rVolt=Volt; %幅值
@ -69,8 +69,10 @@ PDQDi=union(PDi,QDi);
onlyPG=setdiff(PGi,PDQDi); onlyPG=setdiff(PGi,PDQDi);
onlyQG=setdiff(QGi,PDQDi); onlyQG=setdiff(QGi,PDQDi);
%% 计算方差 %% 计算方差
W=(diag([rVolt;rBranchP;rBranchQ;rTransP;rTransQ])*sigma).^2; measureSigma=([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma).^2;
W=sparse(1:length(W),1:length(W),400,length(W),length(W)); 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); stateVarCount=2*length(Volt);
measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ); measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ);
@ -96,6 +98,7 @@ fprintf('
%% 自己写的微分代码 %% 自己写的微分代码
% 以下都是Jacobi矩阵 % 以下都是Jacobi矩阵
SEVolt=sparse(ones(length(mVolt),1)); SEVolt=sparse(ones(length(mVolt),1));
SEVolt(Balance)=rVolt(Balance);
SEVAngle=sparse(zeros(length(mVolt),1)); SEVAngle=sparse(zeros(length(mVolt),1));
maxD=1000; maxD=1000;
Iteration=0; Iteration=0;
@ -213,10 +216,17 @@ while max(abs(g))>1e-5;
G(:,length(mVolt)+Balance)=0; G(:,length(mVolt)+Balance)=0;
G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2); G=G+sparse(length(mVolt)+Balance,length(mVolt)+Balance,1,length(mVolt)*2,length(mVolt)*2);
g(length(mVolt)+Balance)=0; 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; dX=G\-g;
dXStep=1; dXStep=1;
% dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g ); % dXStep=Armijo(z,newwordParameter,W,SEVolt,SEVAngle,dX,g );
maxD=max(abs(dX)) maxD=max(abs(dX))
fprintf('max abs g:%f\n',full(max(abs(g))));
% 更新变量 % 更新变量
SEVolt=SEVolt+dX(1:length(mVolt))*dXStep; SEVolt=SEVolt+dX(1:length(mVolt))*dXStep;
Iteration=Iteration+1; Iteration=Iteration+1;