From 10eb2b56103228ee2ddba358d53c63098ce19439 Mon Sep 17 00:00:00 2001 From: facat Date: Sun, 14 Apr 2013 21:58:12 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E4=BA=86=E7=AD=89=E5=BC=8F=E7=BA=A6?= =?UTF-8?q?=E6=9D=9F=E5=90=8E=E5=BE=88=E9=9A=BE=E6=94=B6=E6=95=9B?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- @SEOpti/SEOpti.m | 10 ++++++++-- @SEOpti/equ.m | 12 +++++++++--- @SEOpti/fun.m | 6 +++--- @SEOpti/init.m | 11 +++++++---- @SEOpti/nle.m | 2 +- @SEOpti/nlrhs.m | 2 +- run.asv | 42 +++++++++++++++++++++++++++--------------- run.m | 26 +++++++++++++++----------- 8 files changed, 71 insertions(+), 40 deletions(-) diff --git a/@SEOpti/SEOpti.m b/@SEOpti/SEOpti.m index 4678c90..653e12b 100644 --- a/@SEOpti/SEOpti.m +++ b/@SEOpti/SEOpti.m @@ -6,8 +6,8 @@ classdef SEOpti mBranchQ=NaN; onlyPG=NaN; onlyQG=NaN; - PG=NaN; - QG=NaN; + mPG=NaN; + mQG=NaN; sigma=NaN; Busnum=NaN; lineI=NaN; @@ -18,6 +18,12 @@ classdef SEOpti zerosInjectionIndex=NaN; cmpY=NaN; Balance=NaN; + mnle=NaN; + mnlrhs=NaN; + Y=NaN; + Yangle=NaN; + r=NaN; + c=NaN; end methods diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m index fa496d7..7c41d19 100644 --- a/@SEOpti/equ.m +++ b/@SEOpti/equ.m @@ -1,12 +1,18 @@ -function [ out_arg ] = equ(this, x ) +function [ out_arg,this ] = equ(this, x ) %EQU Summary of this function goes here % Detailed explanation goes here SEVolt=x(1:this.Busnum); SEVAngel=x(this.Busnum+1:2*this.Busnum); cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); -out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; - +YAngle=sparse(this.r,this.c,SEVAngel(this.r)-SEVAngel(this.c)-this.Yangle,this.Busnum,this.Busnum); +PQ1=diag(SEVolt)*(this.Y.*cos(YAngle))*SEVolt; +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.Busnum+this.Balance)]; +out_arg=full(out_arg); +this.mnle=zeros(2*length(this.zerosInjectionIndex)+1,1); +this.mnlrhs=this.mnle; end diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m index b09204f..284c4e0 100644 --- a/@SEOpti/fun.m +++ b/@SEOpti/fun.m @@ -1,10 +1,10 @@ function [ Objective ] = fun(this, x ) %FUN Summary of this function goes here % Detailed explanation goes here -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; % %% 支路电流 cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 @@ -18,7 +18,7 @@ Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this. %% 0注入节点 PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); %% 发电机注入功率 -Objective=Objective+(this.PG(this.onlyPG)-real(PQ(this.onlyPG)))'*(1./this.sigma^2*eye(length(this.PG(this.onlyPG))))*(this.PG(this.onlyPG)-real(PQ(this.onlyPG))); -Objective=Objective+(this.QG(this.onlyQG)-imag(PQ(this.onlyQG)))'*(1./this.sigma^2*eye(length (this.QG(this.onlyQG))))*(this.QG(this.onlyQG)-imag(PQ(this.onlyQG))); +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))); end diff --git a/@SEOpti/init.m b/@SEOpti/init.m index 5ea3240..ddcf88c 100644 --- a/@SEOpti/init.m +++ b/@SEOpti/init.m @@ -1,4 +1,4 @@ -function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP,mBranchQ ) +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c ) %INIT Summary of this function goes here % Detailed explanation goes here this.mVolt=mVolt; @@ -13,12 +13,15 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio this.cmpY=cmpY; this.onlyPG=onlyPG; this.onlyQG=onlyQG; - this.PG=PG; - this.QG=QG; + this.mPG=mPG; + this.mQG=mQG; this.Balance=Balance; this.mBranchI=mBranchI; this.mBranchP=mBranchP; this.mBranchQ=mBranchQ; - + this.Y=Y; + this.Yangle=Yangle; + this.r=r; + this.c=c; end diff --git a/@SEOpti/nle.m b/@SEOpti/nle.m index a350b73..bb545bd 100644 --- a/@SEOpti/nle.m +++ b/@SEOpti/nle.m @@ -1,7 +1,7 @@ function [ output_args ] = nle( this ) %NLE Summary of this function goes here % Detailed explanation goes here -output_args=zeros(2*length(this.zerosInjectionIndex)+1,1); +output_args=this.mnle; end diff --git a/@SEOpti/nlrhs.m b/@SEOpti/nlrhs.m index 9d1b8b2..ef7dcf4 100644 --- a/@SEOpti/nlrhs.m +++ b/@SEOpti/nlrhs.m @@ -1,7 +1,7 @@ function [ output_args ] = nlrhs( this ) %NLRHS Summary of this function goes here % Detailed explanation goes here -output_args=zeros(2*length(this.zerosInjectionIndex)+1,1); +output_args=this.mnlrhs; end diff --git a/run.asv b/run.asv index c53394f..0c69dff 100644 --- a/run.asv +++ b/run.asv @@ -2,7 +2,7 @@ clear clc yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0'); +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0'); %% 量测量 % 电压 节点电流 支路电流 节点功率 支路功率 %% @@ -39,11 +39,11 @@ mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路 rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量 %% 注入功率 -rPD=PD(PD~=0); +rPD=PD; PDi=find(PD~=0); -rQD=QD(QD~=0); +rQD=QD; QDi=find(QD~=0); -rPG=PG(PG~=0); +rPG=PG; PGi=find(PG~=0); rQG=QG(QG~=0); QGi=find(QG~=0); @@ -53,7 +53,7 @@ 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(mI)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); +measurements=length(mVolt)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); fprintf('冗余度 %f\n',measurements/stateVarCount); %% 进入状态估计计算 % SEVolt=sdpvar(length(Volt),1,'full'); @@ -77,9 +77,9 @@ zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); % zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 %% 发电机注入功率 % 先找到只有发电机的节点 -% PDQDi=union(PDi,QDi); -% onlyPG=setdiff(PGi,PDQDi); -% onlyQG=setdiff(QGi,PDQDi); +PDQDi=union(PDi,QDi); +onlyPG=setdiff(PGi,PDQDi); +onlyQG=setdiff(QGi,PDQDi); % Objective=Objective+(PG(onlyPG)-real(PQ(onlyPG)))'*(1./sigma^2*eye(length(PG(onlyPG))))*(PG(onlyPG)-real(PQ(onlyPG))); % Objective=Objective+(QG(onlyQG)-imag(PQ(onlyQG)))'*(1./sigma^2*eye(length(QG(onlyQG))))*(QG(onlyQG)-imag(PQ(onlyQG))); %% YALMIP求解 @@ -95,13 +95,25 @@ zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); -seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex ); -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2) -x0=ones(Busnum*2,1); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c ); +opts = optiset('solver','ipopt'); +opts.maxiter=2500; +opts.tolrfun=1e-3; +opts.tolafun=1e-3; +opts.warnings='warnings'; +opts.display='iter'; +x0=[ones(Busnum,1);zeros(Busnum,1)]; +[~,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) [x,fval,exitflag,info] = solve(Opt,x0); info -x %% 输出结果 -% fprintf('相对误差\n'); -% (abs(rVolt-double(SEVolt)))./(rVolt) -% MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) +SEVolt=x(1:length(Volt)); +SEVAngel=x(length(Volt)+Balance:2*length(Volt)); +fprintf('目标函数为:%f\n',fval); +fprintf('相对误差\n'); +(abs(rVolt-double(SEVolt)))./(rVolt) +MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) diff --git a/run.m b/run.m index 829b086..31ca8f8 100644 --- a/run.m +++ b/run.m @@ -2,7 +2,7 @@ clear clc yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0'); +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0'); %% 量测量 % 电压 节点电流 支路电流 节点功率 支路功率 %% @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=0.0005;% 标准差 +sigma=0.05;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -39,13 +39,13 @@ mBranchP=rBranchP.*(normrnd(0,sigma,length(rBranchP),1)+1);%支路 rBranchQ=imag((cmpV(lineI)-cmpV(lineJ)).*conj(cmpBranchI)); mBranchQ=rBranchQ.*(normrnd(0,sigma,length(rBranchQ),1)+1);%支路功率量测量 %% 注入功率 -rPD=PD(PD~=0); +rPD=PD; PDi=find(PD~=0); -rQD=QD(QD~=0); +rQD=QD; QDi=find(QD~=0); -rPG=PG(PG~=0); +rPG=PG; PGi=find(PG~=0); -rQG=QG(QG~=0); +rQG=QG; QGi=find(QG~=0); mPD=rPD.*(normrnd(0,sigma,length(rPD),1)+1); mQD=rQD.*(normrnd(0,sigma,length(rQD),1)+1); @@ -95,15 +95,19 @@ onlyQG=setdiff(QGi,PDQDi); %% Opti ToolBox Busnum=length(Volt); seOpti=SEOpti(); -seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,rBranchI,mBranchP,mBranchQ ); -nlrhs=seOpti.nlrhs(); -nle=seOpti.nle(); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,mPG,mQG,Balance,mBranchI,mBranchP,mBranchQ,Y,Yangle,r,c ); opts = optiset('solver','ipopt'); opts.maxiter=2500; +opts.tolrfun=1e-3; +opts.tolafun=1e-3; +opts.warnings='warnings'; opts.display='iter'; -% 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) x0=[ones(Busnum,1);zeros(Busnum,1)]; +[~,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) [x,fval,exitflag,info] = solve(Opt,x0); info %% 输出结果