facat
/
distributionnetwork-power2current
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

换了一种计算序量测量的方式。 

用最大误差作为量测误差。

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2014-11-09 21:52:19 +08:00
parent 77dd602671
commit 92469a3004
1 changed files with 90 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

146
run.m
View File

@ -127,33 +127,46 @@ fprintf('
%% 开始进入状态估计 %% 开始进入状态估计
% clear PD QD PG QG; % clear PD QD PG QG;
%准备量测量 %准备量测量
sigma=0.03;
iterPhaseASpotLoadP=phaseASpotLoadP; iterPhaseASpotLoadP=phaseASpotLoadP;
iterPhaseBSpotLoadP=phaseBSpotLoadP; iterPhaseBSpotLoadP=phaseBSpotLoadP;
iterPhaseCSpotLoadP=phaseCSpotLoadP; iterPhaseCSpotLoadP=phaseCSpotLoadP;
iterPhaseASpotLoadQ=phaseASpotLoadQ; iterPhaseASpotLoadQ=phaseASpotLoadQ;
iterPhaseBSpotLoadQ=phaseBSpotLoadQ; iterPhaseBSpotLoadQ=phaseBSpotLoadQ;
iterPhaseCSpotLoadQ=phaseCSpotLoadQ; iterPhaseCSpotLoadQ=phaseCSpotLoadQ;
%加误差
mphaseASpotLoadP=phaseASpotLoadP.*(1+normrnd(0,sigma,length(phaseASpotLoadP),1));
mphaseBSpotLoadP=phaseBSpotLoadP.*(1+normrnd(0,sigma,length(phaseBSpotLoadP),1));
mphaseCSpotLoadP=phaseCSpotLoadP.*(1+normrnd(0,sigma,length(phaseCSpotLoadP),1));
mphaseASpotLoadQ=phaseASpotLoadQ.*(1+normrnd(0,sigma,length(phaseASpotLoadQ),1));
mphaseBSpotLoadQ=phaseBSpotLoadQ.*(1+normrnd(0,sigma,length(phaseBSpotLoadQ),1));
mphaseCSpotLoadQ=phaseCSpotLoadQ.*(1+normrnd(0,sigma,length(phaseCSpotLoadQ),1));
%电压量测量
mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1));
mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1));
mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1));
%转换得到徐电压的量测量
fV012=Tp2f*conj([mVoltpA';mVoltpB';mVoltpC']);
mfV0=conj(fV012(1,:)');
mfV1=conj(fV012(2,:)');
mfV2=conj(fV012(3,:)');
% VoltpA=sparse(ones(busNum,1));
% VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
% VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
%全部转换为负荷电流 %全部转换为负荷电流
VoltpA=sparse(ones(busNum,1)); mCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA);
VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); mCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB);
VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi); mCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC);
CurpA=conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
CurpB=conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
CurpC=conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
%转换为序电流 %转换为序电流
f012=Tp2f*conj([CurpA';CurpB';CurpC']); mf012=Tp2f*conj([mCurpA';mCurpB';mCurpC']);
%把三序电流分离出来 %把三序电流分离出来
If0=conj(f012(1,:)'); mIf0=conj(mf012(1,:)');
If1=conj(f012(2,:)');%负荷电流 mIf1=conj(mf012(2,:)');%负荷电流
If2=conj(f012(3,:)'); mIf2=conj(mf012(3,:)');
%试着算一下正序电流 mIf0=-mIf0;
% fsY11*V1; mIf1=-mIf1;%mIf1是注入电流,相当于发电机电流
%形成负荷序电流的测量值 mIf2=-mIf2;
%% 设定量测误差 %% 设定量测误差
sigma=0.03;
mIf0=-If0;
mIf1=-If1;%mIf1是注入电流,相当于发电机电流
mIf2=-If2;
%% 先算正序的 %% 先算正序的
%平衡节点电流 %平衡节点电流
fprintf('正序\n'); fprintf('正序\n');
@ -161,19 +174,26 @@ BalI1r=real(-sum(mIf1));
BalI1i=imag(-sum(mIf1)); BalI1i=imag(-sum(mIf1));
%电压 %电压
%制作量测量 %制作量测量
mfsY11=fsY11; % mfsY11=fsY11;
mfsY11(:,Balance)=0; % mfsY11(:,Balance)=0;
mfsY11(Balance,:)=0; % mfsY11(Balance,:)=0;
mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); % mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
rV1=inv(mfsY11)*(mIf1)+1; % rV1=inv(mfsY11)*(mIf1)+1;
sigmaV1=normrnd(0,sigma,length(Loadi),1); % sigmaV1=normrnd(0,sigma,length(Loadi),1);
V1measurement=rV1(Loadi).*(1+sigmaV1); % V1measurement=rV1(Loadi).*(1+sigmaV1);
wV1r=abs(real( rV1(Loadi).*sigmaV1 )); V1measurement=mfV1(Loadi);
wV1i=abs(imag( rV1(Loadi).*sigmaV1 )); % wV1r=abs(real( rV1(Loadi).*sigmaV1 ));
sigmaI1=normrnd(0,sigma,length(Loadi),1); % wV1i=abs(imag( rV1(Loadi).*sigmaV1 ));
I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流 % sigmaI1=normrnd(0,sigma,length(Loadi),1);
wI1r=abs( real(mIf1(Loadi).*sigmaI1) ); % I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流
wI1i=abs( imag(mIf1(Loadi).*sigmaI1) ); I1measurement=mIf1(Loadi);
% wI1r=abs( real(mIf1(Loadi).*sigmaI1) );
% wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
%正序测量误差的第二种形式
wV1r=abs(real(V1measurement)).*sigma;
wV1i=abs(imag(V1measurement)).*sigma;
wI1r=abs(real(I1measurement)).*sigma;
wI1i=abs(imag(I1measurement)).*sigma;
% [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 ); % [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2); % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
% fprintf('目标值 %f\n',full(f)); % fprintf('目标值 %f\n',full(f));
@ -183,19 +203,26 @@ BalI2r=real(-sum(mIf2));
BalI2i=imag(-sum(mIf2)); BalI2i=imag(-sum(mIf2));
%电压 %电压
%制作量测量 %制作量测量
mfsY22=fsY22; % mfsY22=fsY22;
mfsY22(:,Balance)=0; % mfsY22(:,Balance)=0;
mfsY22(Balance,:)=0; % mfsY22(Balance,:)=0;
mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); % mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
rV2=inv(mfsY22)*(mIf2); % rV2=inv(mfsY22)*(mIf2);
sigmaV2=normrnd(0,sigma,length(Loadi),1); % sigmaV2=normrnd(0,sigma,length(Loadi),1);
V2measurement=rV2(Loadi).*(1+sigmaV2); % V2measurement=rV2(Loadi).*(1+sigmaV2);
wV2r=abs(real( rV2(Loadi).*sigmaV2 )); V2measurement=mfV2(Loadi);
wV2i=abs(imag( rV2(Loadi).*sigmaV2 )); % wV2r=abs(real( rV2(Loadi).*sigmaV2 ));
sigmaI2=normrnd(0,sigma,length(Loadi),1); % wV2i=abs(imag( rV2(Loadi).*sigmaV2 ));
I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流 % sigmaI2=normrnd(0,sigma,length(Loadi),1);
wI2r=abs( real(mIf2(Loadi).*sigmaI2) ); % I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流
wI2i=abs( imag(mIf2(Loadi).*sigmaI2) ); I2measurement=mIf2(Loadi);
% wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
% wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
%负序误差的第二种形式
wV2r=abs(real(V2measurement)).*sigma;
wV2i=abs(imag(V2measurement)).*sigma;
wI2r=abs(real(I2measurement)).*sigma;
wI2i=abs(imag(I2measurement)).*sigma;
% [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 ); % [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
% f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2); % f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
% fprintf('目标值 %f\n',full(f)); % fprintf('目标值 %f\n',full(f));
@ -205,19 +232,26 @@ BalI0r=real(-sum(mIf0));
BalI0i=imag(-sum(mIf0)); BalI0i=imag(-sum(mIf0));
%电压 %电压
%制作量测量 %制作量测量
mfsY00=fsY00; % mfsY00=fsY00;
mfsY00(:,Balance)=0; % mfsY00(:,Balance)=0;
mfsY00(Balance,:)=0; % mfsY00(Balance,:)=0;
mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); % mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
rV0=inv(mfsY00)*(mIf0); % rV0=inv(mfsY00)*(mIf0);
sigmaV0=normrnd(0,sigma,length(Loadi),1); % sigmaV0=normrnd(0,sigma,length(Loadi),1);
V0measurement=rV0(Loadi).*(1+sigmaV0); % V0measurement=rV0(Loadi).*(1+sigmaV0);
wV0r=abs(real( rV0(Loadi).*sigmaV0 )); V0measurement=mfV0(Loadi);
wV0i=abs(imag( rV0(Loadi).*sigmaV0 )); % wV0r=abs(real( rV0(Loadi).*sigmaV0 ));
sigmaI0=normrnd(0,sigma,length(Loadi),1); % wV0i=abs(imag( rV0(Loadi).*sigmaV0 ));
I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流 % sigmaI0=normrnd(0,sigma,length(Loadi),1);
wI0r=abs( real(mIf0(Loadi).*sigmaI0) ); % I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流
wI0i=abs( imag(mIf0(Loadi).*sigmaI0) ); I0measurement=mIf0(Loadi);
% wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
% wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
%零序误差的第二种形式
wV0r=abs(real(V0measurement)).*sigma;
wV0i=abs(imag(V0measurement)).*sigma;
wI0r=abs(real(I0measurement)).*sigma;
wI0i=abs(imag(I0measurement)).*sigma;
% matlabpool local 3 % matlabpool local 3
tic tic
for II=1:3 for II=1:3