facat
/
distributionnetwork-power2current
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
distributionnetwork-power2c.../run.m

573 lines
21 KiB
Matlab
Raw Permalink Blame History

function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0,eV2,eV0,eI2,eI0,interation1,interation2,interation0]=run()
%% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。
%
close all
clc
clear
lineZ=readLineZ('feeder104\lineParameter.txt');
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
cap]=dataRead(lineZ,'feeder104\data1.txt');
%负荷不平衡
% for I=1:length(phaseASpotLoadP)
% roll=mod(round(rand()*10),3);
% if roll==0
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
%
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
% end
% if roll==1
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
% end
% if roll==2
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
% end
% end
%% 有几相负荷为0
a=exp(1j*2*pi/3);
Tp2f=1/3*[1 1 1;
1 a a^2;
1 a^2 a];
Tp2f=sparse(Tp2f);
Tf2p=inv(Tp2f);
fsY1amp=abs(fsY1);
[r,c,fsY1ang]=find(fsY1);
fsY1ang=angle(fsY1ang);
Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP;
Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ;
busNum=length(phaseASpotLoadP);
%给序电压赋初值
Vmf1=sparse(ones(busNum,1));
Vaf1=sparse(zeros(busNum,1));
%先求解正序的
PQi=nodeNum;
PG=sparse(busNum,1);
QG=sparse(busNum,1);
QGi=[Balance];
PD=Pabc/3;
QD=Qabc/3;
Loadi=find(PD~=0);
maxD=100000;% 最大不平衡量
EPS=1e-6;
k=0;
kmax=20;
fsY11=fsY1;
fsY00=fsY0;
fsY22=fsY2;
Vf2=sparse(busNum,1);
If2=sparse(busNum,1);
Vf0=sparse(busNum,1);
If0=sparse(busNum,1);
%准备序矩阵
%平衡节点置0置1
fsY2(Balance,:)=0;
fsY2(:,Balance)=0;
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%平衡节点置0置1
fsY0(Balance,:)=0;
fsY0(:,Balance)=0;
fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%%LU分解
[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
%算初始补偿功率
tic
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);
xs=1;
phaseASpotLoadP=phaseASpotLoadP*xs;
phaseBSpotLoadP=phaseBSpotLoadP*xs;
phaseCSpotLoadP=phaseCSpotLoadP*xs;
phaseASpotLoadQ=phaseASpotLoadQ*xs;
phaseBSpotLoadQ=phaseBSpotLoadQ*xs;
phaseCSpotLoadQ=phaseCSpotLoadQ*xs;
while(k<=kmax && maxD> EPS)
k=k+1;
%把补偿电容看作负荷
SA=VoltpA.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,1),busNum,1));
SB=VoltpB.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
%先不要电容
SA=0;
SB=0;
SC=0;
iterPD=PD+real(SA+SB+SC)/3;
iterQD=QD+imag(SA+SB+SC)/3;
iterPhaseASpotLoadP=phaseASpotLoadP+real(SA);
iterPhaseBSpotLoadP=phaseBSpotLoadP+real(SB);
iterPhaseCSpotLoadP=phaseCSpotLoadP+real(SC);
iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA);
iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB);
iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC);
%%做最小二乘法
[dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ...
PQi,PG,QG,QGi,iterPD,iterQD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%不平衡量
maxD=max(abs([dP;dQ;]));
jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
[dV, dVangle]=Solv(busNum,jaco,dP,dQ);%解出修正量
[Vmf1, Vaf1]=Modify(Vmf1,Vaf1,dV,dVangle,1);
fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
%转换为三相电压
VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%用Tp2f\ 代替Tf2p*
VoltpA=conj(VoltpABC(1,:)');
CurpA=-conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
VoltpB=conj(VoltpABC(2,:)');
CurpB=-conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
VoltpC=conj(VoltpABC(3,:)');
CurpC=-conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
f012=Tp2f*conj([CurpA';CurpB';CurpC']);
If0=conj(f012(1,:)');
If1=conj(f012(2,:)');
If2=conj(f012(3,:)');
If0(Balance)=0;
If2(Balance)=0;
%Vf0=fsY0\If0;
Vf0=fsY0Q*(fsY0U\(fsY0L\(fsY0P*(fsY0R\If0))));
%Vf2=fsY2\If2;
Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
fprintf('迭代时间%f\n',toc);
end
FortiscueToc=toc;
fprintf('Fortiscue法计算时间 %f\n',FortiscueToc);
Vf1=Vmf1.*exp(1j*Vaf1);
%% 计算三相不平衡度
%电压不平衡度
eV2=mean(abs(Vf2)./abs(Vf1))/length(Vf1)/3;
eV0=mean(abs(Vf0)./abs(Vf1))/length(Vf1)/3;
%电流不平衡度
eI2=mean(abs(If2(find(If2)))./abs(If1(find(If1))))/length(find(If1))/3;
eI0=mean(abs(If0(find(If0)))./abs(If1(find(If1))))/length(find(If1))/3;
%%
(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率
conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
IpABC=Tf2p*conj([If0';If1';If2']);
%转换回三相电压
VoltpABC=Tf2p*conj([ Vf0'; Vf1'; Vf2']);
rVoltpABC=VoltpABC;
rVoltpA=VoltpA;
rVoltpB=VoltpB;
rVoltpC=VoltpC;
disp([' A B C'])
full(abs(VoltpABC'))
fprintf('节点号对应\n');
disp([setIJ,nodeNum ])
%%检查反推回去的功率是否满足
ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
% clear PD QD PG QG;
%准备量测量
rIf0=If0;
rIf1=If1;
rIf2=If2;
sigma=0.03;
% iterPhaseASpotLoadP=phaseASpotLoadP;
% iterPhaseBSpotLoadP=phaseBSpotLoadP;
% iterPhaseCSpotLoadP=phaseCSpotLoadP;
% iterPhaseASpotLoadQ=phaseASpotLoadQ;
% iterPhaseBSpotLoadQ=phaseBSpotLoadQ;
% 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));
%
mVoltpA=abs(rVoltpA).*(1+normrnd(0,sigma/1,length(VoltpA),1));
mVoltpB=abs(rVoltpB).*(1+normrnd(0,sigma/1,length(VoltpB),1)).*exp(1j*-120/180*pi);
mVoltpC=abs(rVoltpC).*(1+normrnd(0,sigma/1,length(VoltpC),1)).*exp(1j*+120/180*pi);
%
% mVoltpA=sparse(ones(busNum,1));
% mVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
% mVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
%% 开始进入状态估计
tic
%转换得到徐电压的量测量
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);
%全部转换为负荷电流
mCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA);
mCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB);
mCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC);
%转换为序电流
mf012=Tp2f*conj([mCurpA';mCurpB';mCurpC']);
%把三序电流分离出来
mIf0=conj(mf012(1,:)');
mIf1=conj(mf012(2,:)');%负荷电流
mIf2=conj(mf012(3,:)');
mIf0=-mIf0;
mIf1=-mIf1;%mIf1是注入电流,相当于发电机电流
mIf2=-mIf2;
%% 设定量测误差
%% 先算正序的
%平衡节点电流
BalI1r=real(-sum(rIf1));
BalI1i=imag(-sum(rIf1));
%电压
%制作量测量
% mfsY11=fsY11;
% mfsY11(:,Balance)=0;
% mfsY11(Balance,:)=0;
% mfsY11=mfsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV1=inv(mfsY11)*(mIf1)+1;
% sigmaV1=normrnd(0,sigma,length(Loadi),1);
% V1measurement=rV1(Loadi).*(1+sigmaV1);
V1measurement=mfV1(Loadi);
% wV1r=abs(real( rV1(Loadi).*sigmaV1 ));
% wV1i=abs(imag( rV1(Loadi).*sigmaV1 ));
% sigmaI1=normrnd(0,sigma,length(Loadi),1);
% I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流
I1measurement=mIf1(Loadi);
% wI1r=abs( real(mIf1(Loadi).*sigmaI1) );
% wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
%正序测量误差的第二种形式
% wI1r=abs(real(I1measurement)).*sigma;
% wI1i=abs(imag(I1measurement)).*sigma;
IArW=calIWr(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma,mphaseASpotLoadP,mphaseASpotLoadQ );
IAiW=calIWi(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma,mphaseASpotLoadP,mphaseASpotLoadQ );
IBrW=calIWr(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma,mphaseBSpotLoadP,mphaseBSpotLoadQ );
IBiW=calIWi(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma,mphaseBSpotLoadP,mphaseBSpotLoadQ );
ICrW=calIWr(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma,mphaseCSpotLoadP,mphaseCSpotLoadQ );
ICiW=calIWi(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma,mphaseCSpotLoadP,mphaseCSpotLoadQ );
% wI1r=sqrt((IArW+ICrW+IBrW)/9 );
wI1r=sqrt(1/9*(IArW+0.5^2.*IBrW+0.8668^2.*IBiW+0.5^2.*ICrW+0.866^2.*ICiW));
wI1r=wI1r(Loadi);
% wI1i=sqrt((IAiW+ICiW+IBiW)/9 );
wI1i=sqrt(1/9*(IAiW+0.5^2.*IBiW+0.8668^2.*IBrW+0.5^2.*ICiW+0.866^2.*ICrW));
wI1i=wI1i(Loadi);
% wV1r=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV1) );
wV1r=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV1.^2) );
% wV1i=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV1) );
wV1i=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV1.^2) );
wV1r=wV1r(Loadi);
wV1i=wV1i(Loadi);
% [ 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);
% fprintf('目标值 %f\n',full(f));
%% 算负序的
BalI2r=real(-sum(rIf2));
BalI2i=imag(-sum(rIf2));
%电压
%制作量测量
% mfsY22=fsY22;
% mfsY22(:,Balance)=0;
% mfsY22(Balance,:)=0;
% mfsY22=mfsY22+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV2=inv(mfsY22)*(mIf2);
% sigmaV2=normrnd(0,sigma,length(Loadi),1);
% V2measurement=rV2(Loadi).*(1+sigmaV2);
V2measurement=mfV2(Loadi);
% wV2r=abs(real( rV2(Loadi).*sigmaV2 ));
% wV2i=abs(imag( rV2(Loadi).*sigmaV2 ));
% sigmaI2=normrnd(0,sigma,length(Loadi),1);
% I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流
I2measurement=mIf2(Loadi);
% wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
% wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
%负序误差的第二种形式
% wI2r=abs(real(I2measurement)).*sigma;
% wI2i=abs(imag(I2measurement)).*sigma;
% wI2r=sqrt((IArW+ICrW+IBrW)/9 );
wI2r=sqrt(1/9*(IArW+0.5^2.*IBrW+0.8668^2.*IBiW+0.5^2.*ICrW+0.866^2.*ICiW));
wI2r=wI2r(Loadi);
% wI2i=sqrt((IAiW+ICiW+IBiW)/9 );
wI2i=sqrt(1/9*(IAiW+0.5^2.*IBiW+0.8668^2.*IBrW+0.5^2.*ICiW+0.866^2.*ICrW));
wI2i=wI2i(Loadi);
% wV2r=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV2) );
wV2r=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV2.^2) );
% wV2i=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV2) );
wV2i=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV2.^2) );
wV2r=wV2r(Loadi);
wV2i=wV2i(Loadi);
% [ 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);
% fprintf('目标值 %f\n',full(f));
%% 算零序
BalI0r=real(-sum(rIf0));
BalI0i=imag(-sum(rIf0));
%电压
%制作量测量
% mfsY00=fsY00;
% mfsY00(:,Balance)=0;
% mfsY00(Balance,:)=0;
% mfsY00=mfsY00+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
% rV0=inv(mfsY00)*(mIf0);
% sigmaV0=normrnd(0,sigma,length(Loadi),1);
% V0measurement=rV0(Loadi).*(1+sigmaV0);
V0measurement=mfV0(Loadi);
% wV0r=abs(real( rV0(Loadi).*sigmaV0 ));
% wV0i=abs(imag( rV0(Loadi).*sigmaV0 ));
% sigmaI0=normrnd(0,sigma,length(Loadi),1);
% I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流
I0measurement=mIf0(Loadi);
% wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
% wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
%零序误差的第二种形式
%
% wI0r=abs(real(I0measurement)).*sigma;
% wI0i=abs(imag(I0measurement)).*sigma;
% wI0r=sqrt((IArW+ICrW+IBrW)/9 );
wI0r=sqrt(1/9*(IArW+IBrW+ICrW));
wI0r=wI0r(Loadi);
% wI0i=sqrt((IAiW+ICiW+IBiW)/9 );
wI0i=sqrt(1/9*(IAiW+IBiW+ICiW));
wI0i=wI0i(Loadi);
% wV0r=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV0) );
wV0r=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV0.^2) );
% wV0i=sqrt(2*( (abs(mVoltpA).*sigma/3).^2+(abs(mVoltpB).*sigma/3).^2+(abs(mVoltpC).*sigma/3).^2 )/9.*abs(mfV0) );
wV0i=sqrt( (4/9/4*(abs(mVoltpA)+abs(mVoltpC)+abs(mVoltpC)).^2).*abs(mfV0.^2) );
wV0r=wV0r(Loadi);
wV0i=wV0i(Loadi);
% matlabpool local 3
%% 用三相方差确定三序方差
%% 不同量测配置
noLoadi=[1,8,10,12];
% noLoadi=[1:9:max(Loadi)];
% noLoadi=Loadi(1:6:end);
% noLoadi=[1,19];
% noLoadi=[4,3,2,5,6,7,11,14,22];
% noLoadi=Loadi;
noLoadi=[8];
% noLoadi=[5];
% noLoadi=[Loadi(5)];
% noLoadi=Loadi(1:3);
wI0r(ismember(Loadi,noLoadi))=1e100000;
wI0i(ismember(Loadi,noLoadi))=1e100000;
wI1r(ismember(Loadi,noLoadi))=1e100000;
wI1i(ismember(Loadi,noLoadi))=1e100000;
wI2r(ismember(Loadi,noLoadi))=1e100000;
wI2i(ismember(Loadi,noLoadi))=1e100000;
wV0r(ismember(Loadi,noLoadi))=1e100000;
wV0i(ismember(Loadi,noLoadi))=1e100000;
wV1r(ismember(Loadi,noLoadi))=1e100000;
wV1i(ismember(Loadi,noLoadi))=1e100000;
wV2r(ismember(Loadi,noLoadi))=1e100000;
wV2i(ismember(Loadi,noLoadi))=1e100000;
%考虑负荷曲线不准确的情况
guessCurpA=conj((phaseASpotLoadP+1j*phaseASpotLoadQ)./mVoltpA).*(1+unifrnd(-0.15,0.15,length(phaseASpotLoadQ),1));
guessCurpB=conj((phaseBSpotLoadP+1j*phaseBSpotLoadQ)./mVoltpB).*(1+unifrnd(-0.15,0.15,length(phaseBSpotLoadQ),1));
guessCurpC=conj((phaseCSpotLoadP+1j*phaseCSpotLoadQ)./mVoltpC).*(1+unifrnd(-0.15,0.15,length(phaseCSpotLoadQ),1));
guessIf012=Tp2f*conj([-guessCurpA';-guessCurpB';-guessCurpC']);
guessIf0=conj(guessIf012(1,:)');
guessIf1=conj(guessIf012(2,:)');
guessIf2=conj(guessIf012(3,:)');
guessIf0=guessIf0(Loadi);
guessIf1=guessIf1(Loadi);
guessIf2=guessIf2(Loadi);
totalTime1=0;%正序时间
totalTime2=0;%负序时间
totalTime0=0;%零序时间
toc
% parpool(2);
% matlabpool local 3;
V1r1=0;
V1i=0;
I1r=0;
I1i=0;
isConverged1=0;
calcuTime=0;
totalTime1=0;
V2r=0;
V2i=0;
I2r=0;
I2i=0;
isConverged2=0;
totalTime2=0;
V0r=0;
V0i=0;
I0r=0;
I0i=0;
isConverged0=0;
totalTime0=0;
tic
for II=1:3
dI_F=100000000000000;
if II==1
% fprintf('正序\n');
[ V1r,V1i,I1r,I1i,isConverged1,calcuTime,KK ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1,dI_F,2,guessIf1,rIf1,noLoadi );
totalTime1=calcuTime;
interation1=KK;
fprintf('正序时间%f\n',totalTime1);
V1r1=V1r;
end
% end
% while 0
if II==2
fprintf('负序\n');
[ V2r,V2i,I2r,I2i,isConverged2,calcuTime,KK ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0,dI_F,3,guessIf2,rIf2,noLoadi );
totalTime2=calcuTime;
interation2=KK;
fprintf('负序时间%f\n',totalTime2);
end
if II==3
fprintf('零序\n');
[ V0r,V0i,I0r,I0i,isConverged0,calcuTime,KK ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0,dI_F,1,guessIf0,rIf0,noLoadi );
totalTime0=calcuTime;
interation0=KK;
fprintf('零序时间%f\n',totalTime0);
end
end
fprintf('总时间 %f\n',toc);
% matlabpool close
% f=sum(([real(I0measurement);imag(I0measurement)]-[I0r;I0i]).^2)+sum((real(rV0)-V0r).^2)+sum((imag(rV0)-V0i).^2);
% fprintf('目标值 %f\n',full(f));
%% 把三序合成三相
% 三相电压
SEVf0=V0r+1j*V0i;
SEVf1=V1r+1j*V1i;
SEVf2=V2r+1j*V2i;
SEVoltpABC=Tf2p*conj([ SEVf0'; SEVf1'; SEVf2']);
% SEVoltpABC2=Tf2p*conj([ rV0'; rV1'; rV2']);
% 三序电流
SEIf0=I0r+1j*I0i;
SEIf1=I1r+1j*I1i;
SEIf2=I2r+1j*I2i;
SEIpABC=full(Tf2p*conj([SEIf0';SEIf1';SEIf2']));
%看一下差多少
VError=(full(abs(VoltpABC))-abs(SEVoltpABC))./abs(VoltpABC)*100;
VError=reshape(VError,size(VError,1)*size(VError,2),1);
% barPlot( VError, 10,['相对误差%'],['分布密度'],['电压幅值']);
% 三相负荷
rThreeLoad=[ phaseASpotLoadP'+1j*phaseASpotLoadQ';
phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
];
mThreeLoad=[ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
];
% rThreeLoad=rThreeLoad(:,setxor(1:size(SEVoltpABC,2),Balance));
SEThreeLoad=SEVoltpABC(:,Loadi).*conj(-SEIpABC);
% phaseLoadPError=real(rThreeLoad-SEThreeLoad)./real(rThreeLoad)*100;
% phaseLoadQError=imag(rThreeLoad-SEThreeLoad)./imag(rThreeLoad)*100;
% phaseLoadPError=reshape(phaseLoadPError,size(phaseLoadPError,1)*size(phaseLoadPError,2),1);
% phaseLoadQError=reshape(phaseLoadQError,size(phaseLoadQError,1)*size(phaseLoadQError,2),1);
% figure()
% barPlot( phaseLoadPError, 10,['相对误差%'],['分布密度'],['有功负荷误差']);
% figure()
% barPlot( phaseLoadQError, 10,'相对误差%','分布密度','无功负荷误差');
%% 计算统计量
%目标函数均值
rVoltABCV=conj([VoltpA';VoltpB';VoltpC']);
mVoltABCV=conj([mVoltpA';mVoltpB';mVoltpC']);
rPD3P=real([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
]);
rPD3P=rPD3P(:,Loadi);
mPD3P=real([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
]);
mPD3P=mPD3P(:,Loadi);
PD3P=real(SEThreeLoad);
rQD3P=imag([ phaseASpotLoadP'+1j*phaseASpotLoadQ';
phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
]);
rQD3P=rQD3P(:,Loadi);
mQD3P=imag([ mphaseASpotLoadP'+1j*mphaseASpotLoadQ';
mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ';
mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ';
]);
mQD3P=mQD3P(:,Loadi);
QD3P=imag(SEThreeLoad);
JMeasurement=sum(sum((( abs(mVoltABCV)-abs(SEVoltpABC) )./ abs(mVoltABCV)./sigma).^2))+sum(sum(((mPD3P-PD3P)./mPD3P./sigma).^2))+sum(sum(((mQD3P-QD3P)./mQD3P./sigma).^2));
%估计误差统计
%量测量数量
Busnum=busNum;
mCount=Busnum*3+length(Loadi)*3*2;
%估计量质量
AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(SEVoltpABC))))/size(rVoltABCV,1)/size(rVoltABCV,2);
AME_VAngle=sum(sum(abs( angle(rVoltABCV)-angle(SEVoltpABC))))/size(rVoltABCV,1)/size(rVoltABCV,2);
AME_PD=sum(sum(abs(real(SEThreeLoad-rThreeLoad(:,Loadi)))))/length(Loadi)/3;
AME_QD=sum(sum(abs(imag(SEThreeLoad-rThreeLoad(:,Loadi)))))/length(Loadi)/3;
AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(SEVoltpABC))));
AME_VAngle=sum(sum(abs( angle(rVoltABCV)-angle(SEVoltpABC))));
AME_PD=sum(sum(abs(real(SEThreeLoad-rThreeLoad(:,Loadi)))));
AME_QD=sum(sum(abs(imag(SEThreeLoad-rThreeLoad(:,Loadi)))));
%计算与量测值的
AME_mVolt=sum(sum(abs( abs(mVoltABCV)-abs(rVoltABCV))))/size(rVoltABCV,1)/size(rVoltABCV,2);
AME_mPD=sum(sum(abs(real(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))))/length(Loadi)/3;
AME_mQD=sum(sum(abs(imag(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))))/length(Loadi)/3;
AME_mVolt=sum(sum(abs( abs(mVoltABCV)-abs(rVoltABCV))));
AME_mPD=sum(sum(abs(real(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))));
AME_mQD=sum(sum(abs(imag(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))));
%无偏估计
% AME_Volt=abs(rVoltABCV)-abs(SEVoltpABC);
% AME_VAngle= ( (( angle(rVoltABCV)-angle(SEVoltpABC))));
% AME_PD= ( ((real(SEThreeLoad-rThreeLoad(:,Loadi)))));
% AME_QD= ( ((imag(SEThreeLoad-rThreeLoad(:,Loadi)))));
% %计算与量测值的
% AME_mVolt= ( (( abs(mVoltABCV)-abs(rVoltABCV))));
% AME_mPD= ( ((real(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))));
% AME_mQD= ( ((imag(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi)))));
%返回收敛信息
isConverged=isConverged1*isConverged2*isConverged0;
% phaseASpotLoadP=sparse(Loadi,1,real(SEThreeLoad(1,:)),69,1);
% phaseBSpotLoadP=sparse(Loadi,1,real(SEThreeLoad(2,:)),69,1);
% phaseCSpotLoadP=sparse(Loadi,1,real(SEThreeLoad(3,:)),69,1);
% phaseASpotLoadQ=sparse(Loadi,1,imag(SEThreeLoad(1,:)),69,1);
% phaseBSpotLoadQ=sparse(Loadi,1,imag(SEThreeLoad(2,:)),69,1);
% phaseCSpotLoadQ=sparse(Loadi,1,imag(SEThreeLoad(3,:)),69,1);
% ub=checkSSatisfied(Balance,phaseABCY,SEVoltpABC, ...
% phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
% phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
end
Reference in New Issue View Git Blame Copy Permalink