From d8597b5a64d2613dcccfebd5bf9ceaa06e5b100e Mon Sep 17 00:00:00 2001 From: "dmy@lab" Date: Mon, 30 Mar 2015 17:50:51 +0800 Subject: [PATCH] =?UTF-8?q?1.=E6=B7=BB=E5=8A=A0=E7=94=9F=E6=88=90Laplace?= =?UTF-8?q?=E5=88=86=E5=B8=83=E7=9A=84=E5=87=BD=E6=95=B0=202.=E8=80=83?= =?UTF-8?q?=E8=99=91=E8=B4=9F=E8=8D=B7=E5=92=8C=E7=BA=BF=E8=B7=AF=E5=8A=9F?= =?UTF-8?q?=E7=8E=87=E9=87=8F=E6=B5=8B?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dmy@lab --- LineCurrent.m | 7 ++++++ laprnd.m | 30 ++++++++++++++++++++++++++ run.m | 60 +++++++++++++++++++++++++++++++++------------------ 3 files changed, 76 insertions(+), 21 deletions(-) create mode 100644 LineCurrent.m create mode 100644 laprnd.m diff --git a/LineCurrent.m b/LineCurrent.m new file mode 100644 index 0000000..d47624f --- /dev/null +++ b/LineCurrent.m @@ -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 + diff --git a/laprnd.m b/laprnd.m new file mode 100644 index 0000000..ac7a10f --- /dev/null +++ b/laprnd.m @@ -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)); diff --git a/run.m b/run.m index 49c8aa8..c3be1c0 100644 --- a/run.m +++ b/run.m @@ -5,18 +5,21 @@ addpath('.\Powerflow') [~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('E:\算例\feeder33\feeder33ieee.txt', '0'); % 'E:\算例\feeder33\feeder33ieee.txt' %% 开始生成量测量 -sigma=0.01;% 标准差 +sigma=0.03;% 标准差 loop=1; VoltAAE=0; VAngleAAE=0; +LineIndex=[1:16]; while 1 %% 电压 %电压幅值 rVolt=Volt; %幅值 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; %% 电流 %注入电流 @@ -66,19 +69,23 @@ while 1 % 正态分布 -% mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1); -% mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1); -% mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1); -% mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1); + mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1); + mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1); + mPG=rPG.*(normrnd(0,sigma,length(rPG),1)+1); + mQG=rQG.*(normrnd(0,sigma,length(rQG),1)+1); % 均匀分布 - mPD=rPD.*(unifrnd(-3*sigma,3*sigma,length(rPD),1)+1); - mQD=rQD.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1); - mPG=rPG.*(unifrnd(-3*sigma,3*sigma,length(rPG),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)); -% mVolt=rVolt.*(1+unifrnd(-3*sigma,3*sigma,1,length(rVolt))); +% mPD=rPD.*(unifrnd(-3*sigma,3*sigma,length(rPD),1)+1); +% mQD=rQD.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1); +% mPG=rPG.*(unifrnd(-3*sigma,3*sigma,length(rPG),1)+1); +% mQG=rQG.*(unifrnd(-3*sigma,3*sigma,length(rQD),1)+1); + +%Laplace 分布 +% 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注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); @@ -90,7 +97,7 @@ while 1 onlyQG=setdiff(QGi,PDQDi); %% 计算方差 % 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)); W=sparse(diag(1./measureSigma.^2)) ; % W=eye(length(W)); @@ -309,7 +316,12 @@ while 1 % dTQij_dVi+dTQij_dVj,dTQij_dThetai+dTQij_dThetaj];%jacobi H=[dV_dV,dV_dTyta; 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 );%复数支路电流 SEBranchP=BranchP( SEVolt.*exp(1j*SEVAngle),SEBranchI,lineI,lineB2 ); @@ -321,12 +333,12 @@ while 1 SEPD=diag(SEVolt)*Y.* ( spfun (@cos, 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]; % 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; @@ -397,7 +409,7 @@ while 1 end VoltAAE=VoltAAE+sum(abs((SEVolt-rVolt)./rVolt)); VAngleAAE=VAngleAAE+sum(abs((SEVAngle(2:end)-rVAngel(2:end))./rVAngel(2:end))); - if loop>=500 + if loop>=1 break end loop=loop+1; @@ -413,9 +425,15 @@ fprintf(' MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngle) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngle ) % 检查最优性条件 - if any(abs(optimalCondition)>eps) fprintf('最优性条件不满足\n') else fprintf('最优性条件满足\n') -end \ No newline at end of file +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)) \ No newline at end of file