From ae2981658ef4d4a9167faebd18c7bb7dc917b005 Mon Sep 17 00:00:00 2001 From: facat Date: Sun, 14 Apr 2013 20:33:17 +0800 Subject: [PATCH] =?UTF-8?q?=E4=B8=8D=E6=98=AF=E7=94=A8=E7=94=B5=E6=B5=81?= =?UTF-8?q?=E4=BD=9C=E4=B8=BA=E9=87=8F=E6=B5=8B=E9=87=8F=EF=BC=8C=E5=9B=A0?= =?UTF-8?q?=E4=B8=BAabs=E5=8F=AF=E5=AF=BC=E4=B8=8D=E5=A5=BD=E5=A4=84?= =?UTF-8?q?=E7=90=86=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- @SEOpti/SEOpti.m | 37 ++++++++++++++++++------------------- @SEOpti/equ.m | 1 - @SEOpti/fun.m | 39 +++++++++++++++++++-------------------- @SEOpti/init.asv | 24 ++++++++++++++++++++++++ @SEOpti/init.m | 5 ++++- run.m | 15 ++++++++++----- 6 files changed, 75 insertions(+), 46 deletions(-) create mode 100644 @SEOpti/init.asv diff --git a/@SEOpti/SEOpti.m b/@SEOpti/SEOpti.m index 883f5d3..4678c90 100644 --- a/@SEOpti/SEOpti.m +++ b/@SEOpti/SEOpti.m @@ -1,24 +1,23 @@ classdef SEOpti properties - mVolt=0; - mVAngel=0; - mBranchI=0; - mBranchP=0; - mBranchQ=0; - onlyPG=0; - onlyQG=0; - PG=0; - QG=0; - sigma=0; - Busnum=0; - lineI=0; - lineJ=0; - lineR=0; - lineX=0; - lineB2=0; - zerosInjectionIndex=0; - cmpY=0; - Balance=0; + mVolt=NaN; + mBranchI=NaN; + mBranchP=NaN; + mBranchQ=NaN; + onlyPG=NaN; + onlyQG=NaN; + PG=NaN; + QG=NaN; + sigma=NaN; + Busnum=NaN; + lineI=NaN; + lineJ=NaN; + lineR=NaN; + lineX=NaN; + lineB2=NaN; + zerosInjectionIndex=NaN; + cmpY=NaN; + Balance=NaN; end methods diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m index 88eb2ca..fa496d7 100644 --- a/@SEOpti/equ.m +++ b/@SEOpti/equ.m @@ -5,7 +5,6 @@ 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); -% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0]; out_arg=[real(PQ(this.zerosInjectionIndex));imag(PQ(this.zerosInjectionIndex));]; out_arg=[out_arg;x(this.Busnum+this.Balance)]; diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m index 01ea66f..b09204f 100644 --- a/@SEOpti/fun.m +++ b/@SEOpti/fun.m @@ -1,25 +1,24 @@ function [ Objective ] = fun(this, x ) %FUN Summary of this function goes here % Detailed explanation goes here - 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);%电压 - %% 支路电流 - cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 - cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 - SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 - Objective=Objective+(SEBranchI-this.mBranchI)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI);%%电流 - %% 支路功率 - SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(SEBranchI)); - SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(SEBranchI)); - Objective=Objective+(SEBranchP-this.mBranchP)'*(1./this.sigma^2*eye(length(this.mBranchP)))*(SEBranchP-this.mBranchP); - Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this.mBranchQ)))*(SEBranchQ-this.mBranchQ); -% %% 0注入节点 - PQ=diag(cmpSEV)*conj(this.cmpY*cmpSEV); -% zeroInjP=real(PQ(this.zerosInjectionIndex));%% 0注入节点 -% zeroInjQ=imag(PQ(this.zerosInjectionIndex));%% 0注入节点 - %% 发电机注入功率 - 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=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);%电压 +% %% 支路电流 +cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 +cmpSEBranchI=(cmpSEV(this.lineI)-cmpSEV(this.lineJ))./(this.lineR+1j*this.lineX);%复数支路电流 +% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 +% Objective=Objective+(SEBranchI-this.mBranchI)'*(1./this.sigma^2*eye(length(this.mBranchI)))*(SEBranchI-this.mBranchI);%%电流 +%% 支路功率 +SEBranchP=real((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); +SEBranchQ=imag((cmpSEV(this.lineI)-cmpSEV(this.lineJ)).*conj(cmpSEBranchI)); +Objective=Objective+(SEBranchP-this.mBranchP)'*(1./this.sigma^2*eye(length(this.mBranchP)))*(SEBranchP-this.mBranchP); +Objective=Objective+(SEBranchQ-this.mBranchQ)'*(1./this.sigma^2*eye(length(this.mBranchQ)))*(SEBranchQ-this.mBranchQ); +%% 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))); end diff --git a/@SEOpti/init.asv b/@SEOpti/init.asv new file mode 100644 index 0000000..0d22066 --- /dev/null +++ b/@SEOpti/init.asv @@ -0,0 +1,24 @@ +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP ) +%INIT Summary of this function goes here +% Detailed explanation goes here + this.mVolt=mVolt; + this.sigma=sigma; + this.Busnum=Busnum; + this.lineI=newwordParameter.line.lineI; + this.lineJ=newwordParameter.line.lineJ; + this.lineR=newwordParameter.line.lineR; + this.lineX=newwordParameter.line.lineX; + this.lineB2=newwordParameter.line; + this.zerosInjectionIndex=zerosInjectionIndex; + this.cmpY=cmpY; + this.onlyPG=onlyPG; + this.onlyQG=onlyQG; + this.PG=PG; + this.QG=QG; + this.Balance=Balance; + this.mBranchI=mBranchI; + this.mBranchP=mBranchP; + this.mBranchQ=mBranchQ; + +end + diff --git a/@SEOpti/init.m b/@SEOpti/init.m index 1b233e4..5ea3240 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 ) +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,mBranchI,mBranchP,mBranchQ ) %INIT Summary of this function goes here % Detailed explanation goes here this.mVolt=mVolt; @@ -16,6 +16,9 @@ function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectio this.PG=PG; this.QG=QG; this.Balance=Balance; + this.mBranchI=mBranchI; + this.mBranchP=mBranchP; + this.mBranchQ=mBranchQ; end diff --git a/run.m b/run.m index a7b264d..829b086 100644 --- a/run.m +++ b/run.m @@ -10,7 +10,7 @@ addpath('.\Powerflow') % 电压 相角 %% %% 开始生成量测量 -sigma=0.05;% 标准差 +sigma=0.0005;% 标准差 %% 电压 %电压幅值 rVolt=Volt; %幅值 @@ -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'); @@ -95,12 +95,17 @@ 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); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance,rBranchI,mBranchP,mBranchQ ); nlrhs=seOpti.nlrhs(); nle=seOpti.nle(); -Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle) -x0=ones(Busnum*2,1); +opts = optiset('solver','ipopt'); +opts.maxiter=2500; +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)]; [x,fval,exitflag,info] = solve(Opt,x0); +info %% 输出结果 SEVolt=x(1:length(Volt)); SEVAngel=x(length(Volt)+Balance:2*length(Volt));