facat
/
stateestimation-self-derivatives-measurementdistribution
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

1.添加生成Laplace分布的函数 

2.考虑负荷和线路功率量测

Signed-off-by: dmy@lab <dmy@lab.lab>
This commit is contained in:
dmy@lab 2015-03-30 17:50:51 +08:00
parent d3146da2b9
commit d8597b5a64
3 changed files with 76 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
LineCurrent.m Normal file
View File

@ -0,0 +1,7 @@
function [ output_args ] = LineCurrent( Linei,Linej,Liner,Linex,Volt,VAngle )
cmpV=Volt.*exp(1j*VAngle);
GB=1./(Liner+1j*Linex);
cmpI=( cmpV(Linei)-cmpV(Linej) ).*GB;
output_args=abs(cmpI);
end

30
laprnd.m Normal file
View File

@ -0,0 +1,30 @@
function y = laprnd(m, n, mu, sigma)
%LAPRND generate i.i.d. laplacian random number drawn from laplacian distribution
% with mean mu and standard deviation sigma.
% mu : mean
% sigma : standard deviation
% [m, n] : the dimension of y.
% Default mu = 0, sigma = 1.
% For more information, refer to
% http://en.wikipedia.org./wiki/Laplace_distribution
% Author : Elvis Chen (bee33@sjtu.edu.cn)
% Date : 01/19/07
%Check inputs
if nargin < 2
error('At least two inputs are required');
end
if nargin == 2
mu = 0; sigma = 1;
end
if nargin == 3
sigma = 1;
end
% Generate Laplacian noise
u = rand(m, n)-0.5;
b = sigma / sqrt(2);
y = mu - b * sign(u).* log(1- 2* abs(u));

60
run.m
View File

@ -5,18 +5,21 @@ addpath('.\Powerflow')
[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\算例\feeder33\feeder33ieee.txt', '0'); [~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\算例\feeder33\feeder33ieee.txt', '0');
% 'E:\算例\feeder33\feeder33ieee.txt' % 'E:\算例\feeder33\feeder33ieee.txt'
%% 开始生成量测量 %% 开始生成量测量
sigma=0.01;% ±ê×¼²î sigma=0.03;% 标准差
loop=1; loop=1;
VoltAAE=0; VoltAAE=0;
VAngleAAE=0; VAngleAAE=0;
LineIndex=[1:16];
while 1 while 1
%% 电压 %% 电压
%电压幅值 %电压幅值
rVolt=Volt; %幅值 rVolt=Volt; %幅值
BalanceVolt=Volt(Balance); BalanceVolt=Volt(Balance);
% mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%µçѹÁ¿²âÁ¿ mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量
%均匀分布 %均匀分布
mVolt=rVolt.*(unifrnd(-3*sigma,3*sigma,length(rVolt),1)+1); % mVolt=rVolt.*(unifrnd(-3*sigma,3*sigma,length(rVolt),1)+1);
%Laplace分布
% mVolt=rVolt.*(laprnd(length(rVolt),1,0,sigma)+1);
rVAngel=Vangle; rVAngel=Vangle;
%% 电流 %% 电流
%注入电流 %注入电流
@ -66,19 +69,23 @@ while 1
% 正态分布 % 正态分布
% mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1); mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1);
% mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1); mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1);
% mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1); mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1);
% mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1); mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1);
% 均匀分布 % 均匀分布
mPD=rPD.*(unifrnd(-3*sigma,3*sigma,length(rPD),1)+1); % mPD=rPD.*(unifrnd(-3*sigma,3*sigma,length(rPD),1)+1);
mQD=rQD.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1); % mQD=rQD.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1);
mPG=rPG.*(unifrnd(-3*sigma,3*sigma,length(rPG),1)+1); % mPG=rPG.*(unifrnd(-3*sigma,3*sigma,length(rPG),1)+1);
mQG=rQG.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1); % mQG=rQG.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1);
% PD0(Loadi)=RealPD(Loadi).*(1+unifrnd(-3*sigma,3*sigma,length(Loadi),1));
% QD0(Loadi)=RealQD(Loadi).*(1+unifrnd(-3*sigma,3*sigma,length(Loadi),1)); %Laplace 分布
% mVolt=rVolt.*(1+unifrnd(-3*sigma,3*sigma,1,length(rVolt))); % mPD=rPD.*(laprnd(length(rPD),1,0,sigma)+1);
% mQD=rQD.*(laprnd(length(rQD),1,0,sigma)+1);
% mPG=rPG.*(laprnd(length(rPG),1,0,sigma)+1);
% mQG=rQG.*(laprnd(length(rQD),1,0,sigma)+1);
%% 0注入节点 %% 0注入节点
zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=1:length(Volt);
zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) );
@ -90,7 +97,7 @@ while 1
onlyQG=setdiff(QGi,PDQDi); onlyQG=setdiff(QGi,PDQDi);
%% 计算方差 %% 计算方差
% measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma)); % measureSigma=abs(([rVolt;rBranchP;rBranchQ;rTransP;rTransQ].*sigma));
measureSigma=abs(([rVolt;rPD(PDi);rQD(QDi);].*sigma)); measureSigma=abs(([rVolt;rPD(PDi);rQD(QDi);rBranchP(LineIndex);rBranchQ(LineIndex);rTransP;rTransQ].*sigma));
% measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6)); % measureSigma(measureSigma<1e-6)=mean(measureSigma(measureSigma>1e-6));
W=sparse(diag(1./measureSigma.^2)) ; W=sparse(diag(1./measureSigma.^2)) ;
% W=eye(length(W)); % W=eye(length(W));
@ -309,7 +316,12 @@ while 1
% dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi % dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi
H=[dV_dV,dV_dTyta; H=[dV_dV,dV_dTyta;
dPdV(PDi,:),dPdTyta(PDi,:); dPdV(PDi,:),dPdTyta(PDi,:);
dQdV(QDi,:),dQdTyta(QDi,:)];%jacobi dQdV(QDi,:),dQdTyta(QDi,:);
dLPij_dVi(LineIndex,:)+dLPij_dVj(LineIndex,:),dLPij_dThetai(LineIndex,:)+dLPij_dThetaj(LineIndex,:) ;
dLQij_dVi(LineIndex,:)+dLQij_dVj(LineIndex,:),dLQij_dThetai(LineIndex,:)+dLQij_dThetaj(LineIndex,:) ;
dTPij_dVi+dTPij_dVj,dTPij_dThetai+dTPij_dThetaj;
dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj
];%jacobi
SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流 SEBranchI=BranchI( SEVolt.*exp(1j*SEVAngle),lineI,lineJ,lineR,lineX );%复数支路电流
SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 ); SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 );
@ -321,12 +333,12 @@ while 1
SEPD=diag(SEVolt)*Y.* ( spfun (@cos, VAngleIJ ) )*SEVolt; SEPD=diag(SEVolt)*Y.* ( spfun (@cos, VAngleIJ ) )*SEVolt;
SEQD=diag(SEVolt)*Y.* ( spfun(@sin,VAngleIJ) )*SEVolt; SEQD=diag(SEVolt)*Y.* ( spfun(@sin,VAngleIJ) )*SEVolt;
h=[SEVolt;SEPD(PDi);SEQD(QDi);]; h=[SEVolt;SEPD(PDi);SEQD(QDi);SEBranchP(LineIndex);SEBranchQ(LineIndex);SETransP;SETransQ];
% h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ]; % h=[SEVolt;SEBranchP;SEBranchQ;SETransP;SETransQ];
% z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ]; % z=[mVolt;mBranchP;mBranchQ;mTransP;mTransQ];
z=[mVolt;-mPD(PDi);-mQD(QDi)]; z=[mVolt;-mPD(PDi);-mQD(QDi);mBranchP(LineIndex);mBranchQ(LineIndex);mTransP;mTransQ];
G=H'*W*H; G=H'*W*H;
@ -397,7 +409,7 @@ while 1
end end
VoltAAE=VoltAAE+sum(abs((SEVolt-rVolt)./rVolt)); VoltAAE=VoltAAE+sum(abs((SEVolt-rVolt)./rVolt));
VAngleAAE=VAngleAAE+sum(abs((SEVAngle(2:end)-rVAngel(2:end))./rVAngel(2:end))); VAngleAAE=VAngleAAE+sum(abs((SEVAngle(2:end)-rVAngel(2:end))./rVAngel(2:end)));
if loop>=500 if loop>=1
break break
end end
loop=loop+1; loop=loop+1;
@ -413,9 +425,15 @@ fprintf('
MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle) MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle)
plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle ) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle )
% 检查最优性条件 % 检查最优性条件
if any(abs(optimalCondition)>eps) if any(abs(optimalCondition)>eps)
fprintf('最优性条件不满足\n') fprintf('最优性条件不满足\n')
else else
fprintf('最优性条件满足\n') fprintf('最优性条件满足\n')
end end
Linei=newwordParameter.line.lineI;
Linej=newwordParameter.line.lineJ;
Liner=newwordParameter.line.lineR;
Linex=newwordParameter.line.lineX;
SELineCurrent=LineCurrent( Linei,Linej,Liner,Linex,SEVolt,SEVAngle );
RealLineCurrent=LineCurrent( Linei,Linej,Liner,Linex,rVolt,rVAngel );
sqrt(sum((SELineCurrent-RealLineCurrent).^2)/length(SELineCurrent))