diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m index d1e95d5..3e286bb 100644 --- a/@SEOpti/equ.m +++ b/@SEOpti/equ.m @@ -10,9 +10,10 @@ PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); % PQ2=diag(SEVolt)*(this.Y.*sin(YAngle))*SEVolt; out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; % out_arg=[PQ1(this.zerosInjectionIndex);PQ2(this.zerosInjectionIndex);]; +% out_arg=[out_arg;x(this.Balance+this.Busnum)]; out_arg=[out_arg;x(this.Balance+this.Busnum)]; out_arg=full(out_arg); -this.mnle=zeros(2*length(this.zerosInjectionIndex)+1,1); +this.mnle=zeros(length(out_arg),1); this.mnlrhs=this.mnle; end diff --git a/@SEOpti/fun.asv b/@SEOpti/fun.asv index f23e6b8..e3682ed 100644 --- a/@SEOpti/fun.asv +++ b/@SEOpti/fun.asv @@ -4,9 +4,8 @@ function [ Objective ] = fun(this, x ) Objective=0; SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); -Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt-this.mVolt);%电压 -% Objective=0; -% %% 支路电流 +Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%电压 +%% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 @@ -17,12 +16,16 @@ SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2); Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2); %% 变压器 -newwordParameter=th -transP=TransPower( newwordParameter,cmpY,Volt,VAngel ); +newwordParameter=this.newwordParameter; +cmpY=this.cmpY; +transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); +Objective=Objective+sum((transP-this.mTransP).^2/this.sigma^2); +transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel ); +Objective=Objective+sum((transQ-this.mTransQ).^2/this.sigma^2); %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 -% Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.onlyPG)-real(PQ(this.onlyPG))); +Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigmas.onlyPG)-real(PQ(this.onlyPG))); % Objective=Objective+(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length(this.mQG(this.onlyQG))))*(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG))); end diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m index 47fee30..ebe5058 100644 --- a/@SEOpti/fun.m +++ b/@SEOpti/fun.m @@ -4,29 +4,30 @@ function [ Objective ] = fun(this, x ) Objective=0; SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); -Objective=(SEVolt-this.mVolt)'*(1./this.sigma^2*eye(length(this.mVolt)))*(SEVolt-this.mVolt);%电压 -% Objective=0; -% %% 支路电流 +Objective=sum(((SEVolt-this.mVolt)./this.sigma).^2);%电压 +%% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -Objective=Objective+sum((SEBranchI-this.mBranchI.^2).^2./this.sigma^2);%%电流 +% (SEBranchI-this.mBranchI).^2./(this.sigma^2) +Objective=Objective+sum((SEBranchI-this.mBranchI).^2./(this.sigma^2));%%电流 %% 线路功率 SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); -Objective=Objective+sum((SEBranchP-this.mBranchP).^2/this.sigma^2); -Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/this.sigma^2); +% (SEBranchP-this.mBranchP).^2/(this.sigma^2)./(this.mBranchP.^2); +Objective=Objective+sum((SEBranchP-this.mBranchP).^2/(this.sigma^2)); +Objective=Objective+sum((SEBranchQ-this.mBranchQ).^2/(this.sigma^2)); %% 变压器 newwordParameter=this.newwordParameter; cmpY=this.cmpY; transP=TransPower( newwordParameter,cmpY,SEVolt,SEVAngel ); -Objective=Objective+sum((transP-this.mTransP).^2/this.sigma^2); +Objective=Objective+sum((transP-this.mTransP).^2/(this.sigma^2)); transQ=TransReactivePower( newwordParameter,cmpY,SEVolt,SEVAngel ); -Objective=Objective+sum((transQ-this.mTransQ).^2/this.sigma^2); +Objective=Objective+sum((transQ-this.mTransQ).^2/(this.sigma^2)); %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 -% Objective=Objective+(this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.mPG(this.onlyPG))))*(this.mPG(this.onlyPG)-real(PQ(this.onlyPG))); -% Objective=Objective+(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length(this.mQG(this.onlyQG))))*(this.mQG(this.onlyQG)-imag(PQ(this.onlyQG))); +% Objective=Objective+sum(((this.mPG(this.onlyPG)-real(PQ(this.onlyPG)))/this.sigma).^2); +% Objective=Objective+sum(((this.mQG(this.onlyQG)-imag(PQ(this.onlyQG)))./this.sigma).^2); end diff --git a/TransPower.m b/TransPower.m index 41fd53c..e0b924c 100644 --- a/TransPower.m +++ b/TransPower.m @@ -7,6 +7,6 @@ cmpSEV=Volt.*exp(1j*VAngel); % % cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt)); cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ)); SETransP=real((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI)); -output_args=SETransP; +output_args=full(SETransP); end diff --git a/TransReactivePower.m b/TransReactivePower.m index 7ca1d2f..bbffe3c 100644 --- a/TransReactivePower.m +++ b/TransReactivePower.m @@ -7,6 +7,6 @@ cmpSEV=Volt.*exp(1j*VAngel); % % cmpY=sparse(r,c,Y.*(1j*YAngel),length(Volt),length(Volt)); cmpSETransI=(cmpSEV(transI)-cmpSEV(transJ)).*cmpY(sub2ind(size(cmpY),transI,transJ)); SETransQ=imag((cmpSEV(transI)-cmpSEV(transJ)).*conj(cmpSETransI)); -output_args=SETransQ; +output_args=full(SETransQ); end diff --git a/mBranchI.mat b/mBranchI.mat new file mode 100644 index 0000000..502780e Binary files /dev/null and b/mBranchI.mat differ diff --git a/mBranchP.mat b/mBranchP.mat new file mode 100644 index 0000000..1ba197c Binary files /dev/null and b/mBranchP.mat differ diff --git a/mBranchQ.mat b/mBranchQ.mat new file mode 100644 index 0000000..d3655ba Binary files /dev/null and b/mBranchQ.mat differ diff --git a/mPG.mat b/mPG.mat new file mode 100644 index 0000000..c08e25a Binary files /dev/null and b/mPG.mat differ diff --git a/mQG.mat b/mQG.mat new file mode 100644 index 0000000..6df8a65 Binary files /dev/null and b/mQG.mat differ diff --git a/mTransP.mat b/mTransP.mat new file mode 100644 index 0000000..57d0292 Binary files /dev/null and b/mTransP.mat differ diff --git a/mTransQ.mat b/mTransQ.mat new file mode 100644 index 0000000..f7a3b7a Binary files /dev/null and b/mTransQ.mat differ diff --git a/mVolt.mat b/mVolt.mat new file mode 100644 index 0000000..492b786 Binary files /dev/null and b/mVolt.mat differ diff --git a/run.asv b/run.asv index 982b7a1..c5863ff 100644 --- a/run.asv +++ b/run.asv @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=1;% 标准差 +sigma=0.01;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -30,7 +30,7 @@ lineJ=newwordParameter.line.lineJ; lineR=newwordParameter.line.lineR; lineX=newwordParameter.line.lineX; lineB2=newwordParameter.line; -transI= + cmpBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流 rBranchI=abs(cmpBranchI);% 支路电流幅值 mBranchI=rBranchI.*(normrnd(0,sigma,length(rBranchI),1)+1);%支路电流量测量 @@ -52,10 +52,6 @@ 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); -%% 冗余度计算 -stateVarCount=2*length(Volt); -measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); -fprintf('冗余度 %f\n',measurements/stateVarCount); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); @@ -69,9 +65,31 @@ mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1); PDQDi=union(PDi,QDi); onlyPG=setdiff(PGi,PDQDi); onlyQG=setdiff(QGi,PDQDi); +%% 冗余度计算 +stateVarCount=2*length(Volt); +measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ); +fprintf('冗余度 %f\n',measurements/stateVarCount); %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); +%% save +% save('mVolt','mVolt'); +% save('mPG','mPG'); +% save('mQG','mQG'); +% save('mBranchI','mBranchI'); +% save('mBranchP','mBranchP'); +% save('mBranchQ','mBranchQ'); +% save('mTransP','mTransP'); +% save('mTransQ','mTransQ'); +%% load +load('mVolt'); +load('mPG'); +load('mQG'); +load('mBranchI'); +load('mBranchP'); +load('mBranchQ'); +load('mTransP'); +load('mTransQ'); seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ); opts = optiset('solver','ipopt'); opts.maxiter=85500; @@ -82,12 +100,13 @@ opts.tolrfun=1e-4; opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; -x0=[ones(Busnum,1);zeros(Busnum,1)]; +% x0=[ones(Busnum,1);zeros(Busnum,1)]; +x0=[rVolt;rVAngle]; [~,seOpti]=seOpti.equ(x0); nlrhs=seOpti.nlrhs(); nle=seOpti.nle(); -% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) +Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info %% 输出结果 @@ -98,3 +117,4 @@ fprintf(' (abs(rVolt-double(SEVolt)))./(rVolt) MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel ) +seOpti.fun(x); \ No newline at end of file diff --git a/run.m b/run.m index 791e6a0..830d046 100644 --- a/run.m +++ b/run.m @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=1;% 标准差 +sigma=0.03;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -52,10 +52,6 @@ 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); -%% 冗余度计算 -stateVarCount=2*length(Volt); -measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); -fprintf('冗余度 %f\n',measurements/stateVarCount); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); @@ -69,9 +65,31 @@ mTransQ=rTransQ.*(normrnd(0,sigma,length(rTransQ),1)+1); PDQDi=union(PDi,QDi); onlyPG=setdiff(PGi,PDQDi); onlyQG=setdiff(QGi,PDQDi); +%% 冗余度计算 +stateVarCount=2*length(Volt); +measurements=length(mVolt)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG)+length(mTransP)+length(mTransQ); +fprintf('冗余度 %f\n',measurements/stateVarCount); %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); +%% save +% save('mVolt','mVolt'); +% save('mPG','mPG'); +% save('mQG','mQG'); +% save('mBranchI','mBranchI'); +% save('mBranchP','mBranchP'); +% save('mBranchQ','mBranchQ'); +% save('mTransP','mTransP'); +% save('mTransQ','mTransQ'); +%% load +load('mVolt'); +load('mPG'); +load('mQG'); +load('mBranchI'); +load('mBranchP'); +load('mBranchQ'); +load('mTransP'); +load('mTransQ'); seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,mTransP,mTransQ); opts = optiset('solver','ipopt'); opts.maxiter=85500; @@ -82,12 +100,13 @@ opts.tolrfun=1e-4; opts.tolafun=1e-4; opts.warnings='all'; opts.display='off'; -x0=[ones(Busnum,1);zeros(Busnum,1)]; +x0=[ones(Busnum,1);-0.5*ones(Busnum,1)]; +% x0=[rVolt;rVAngel]; [~,seOpti]=seOpti.equ(x0); nlrhs=seOpti.nlrhs(); nle=seOpti.nle(); -% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) +Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) +% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'options',opts) [x,fval,exitflag,info] = solve(Opt,x0); info %% 输出结果 @@ -98,3 +117,4 @@ fprintf(' (abs(rVolt-double(SEVolt)))./(rVolt) MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) plotAndComparison( rVolt,rVAngel,SEVolt,SEVAngel ) +seOpti.fun(x); \ No newline at end of file