From 32f96a94f09bf8ff733d4c998352ca72cf502373 Mon Sep 17 00:00:00 2001 From: facat Date: Sun, 14 Apr 2013 19:40:35 +0800 Subject: [PATCH] =?UTF-8?q?=E6=94=B9=E7=94=A8Opti=20Toolbox?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- @SEOpti/SEOpti.m | 33 ++++++++++++++++++++++ @SEOpti/equ.m | 13 +++++++++ @SEOpti/fun.m | 25 +++++++++++++++++ @SEOpti/init.m | 21 ++++++++++++++ @SEOpti/nle.m | 7 +++++ @SEOpti/nlrhs.m | 7 +++++ MaxDeviation.m | 8 ++++++ run.asv | 72 +++++++++++++++++++++++++++++++----------------- run.m | 64 ++++++++++++++++++++++++++++-------------- 9 files changed, 203 insertions(+), 47 deletions(-) create mode 100644 @SEOpti/SEOpti.m create mode 100644 @SEOpti/equ.m create mode 100644 @SEOpti/fun.m create mode 100644 @SEOpti/init.m create mode 100644 @SEOpti/nle.m create mode 100644 @SEOpti/nlrhs.m create mode 100644 MaxDeviation.m diff --git a/@SEOpti/SEOpti.m b/@SEOpti/SEOpti.m new file mode 100644 index 0000000..883f5d3 --- /dev/null +++ b/@SEOpti/SEOpti.m @@ -0,0 +1,33 @@ +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; + end + + methods + function [this]=SEOpti() + + end + + end + + +end + diff --git a/@SEOpti/equ.m b/@SEOpti/equ.m new file mode 100644 index 0000000..88eb2ca --- /dev/null +++ b/@SEOpti/equ.m @@ -0,0 +1,13 @@ +function [ out_arg ] = 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); +% 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)]; +end + diff --git a/@SEOpti/fun.m b/@SEOpti/fun.m new file mode 100644 index 0000000..01ea66f --- /dev/null +++ b/@SEOpti/fun.m @@ -0,0 +1,25 @@ +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))); +end + diff --git a/@SEOpti/init.m b/@SEOpti/init.m new file mode 100644 index 0000000..1b233e4 --- /dev/null +++ b/@SEOpti/init.m @@ -0,0 +1,21 @@ +function [ this ] = init(this,Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY,onlyPG,onlyQG,PG,QG,Balance ) +%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; + +end + diff --git a/@SEOpti/nle.m b/@SEOpti/nle.m new file mode 100644 index 0000000..a350b73 --- /dev/null +++ b/@SEOpti/nle.m @@ -0,0 +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); + +end + diff --git a/@SEOpti/nlrhs.m b/@SEOpti/nlrhs.m new file mode 100644 index 0000000..9d1b8b2 --- /dev/null +++ b/@SEOpti/nlrhs.m @@ -0,0 +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); + +end + diff --git a/MaxDeviation.m b/MaxDeviation.m new file mode 100644 index 0000000..ad6f95a --- /dev/null +++ b/MaxDeviation.m @@ -0,0 +1,8 @@ +function [output_arg]=MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) +t1=[rVolt;rVAngel]; +t2=[SEVolt;SEVAngel]; +t3=(t1(t1~=0)-double(t2(t1~=0)))./t1(t1~=0); +t4=abs(t3); +output_arg=max(t4); +end + diff --git a/run.asv b/run.asv index 3e7cdaa..c53394f 100644 --- a/run.asv +++ b/run.asv @@ -1,5 +1,8 @@ +clear +clc +yalmip('clear') addpath('.\Powerflow') -[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee30.dat', '0'); +[~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee4.dat', '0'); %% 量测量 % 电压 节点电流 支路电流 节点功率 支路功率 %% @@ -12,6 +15,7 @@ sigma=0.05;% %电压幅值 rVolt=Volt; %幅值 mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量 +rVAngel=Vangle; %% 电流 %注入电流 cmpY=Y.*exp(1j*sparse(r,c,Yangle,length(Y),length(Y)));%复数导纳矩阵 @@ -52,36 +56,52 @@ stateVarCount=2*length(Volt); measurements=length(mVolt)+length(mI)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); fprintf('冗余度 %f\n',measurements/stateVarCount); %% 进入状态估计计算 -SEVolt=sdpvar(length(Volt),1); -SEVAngel=sdpvar(length(Vangle),1); -Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%%电压 +% SEVolt=sdpvar(length(Volt),1,'full'); +% SEVAngel=sdpvar(length(Vangle),1,'full'); +% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压 %% 支路电流 -cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 -cmpSEBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流 -SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 +% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 +% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流 +% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 +% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 %% 支路功率 -SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); -Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); +% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); +% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); +% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); +% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); -PQ=diag(SEVolt)*conj(cmpY*SEVolt); -zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 -zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 +% PQ=diag(cmpSEV)*cmpY*cmpSEV; +% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 +% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 %% 发电机注入功率 % 先找到只有发电机的节点 -PDQDi=union(PDi,QDi); -onlyPG=setdiff(PGi,PDQDi); -onlyQG=setdiff(QGi,PDQDi); -PG-real(PQ()) +% 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求解 -Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; -% Constraints=[[zeros(length(c)) A' -eye(length(lbounds))]*x==-c;[A zeros(length(b)) zeros(length(b),length(lbounds))]*x<=b;10>=x>=0]; -options = sdpsettings('verbose',2,'solver','ipopt'); -solvesdp(Constraints,Objective,options) -double(Objective) -fprintf('相对误差\n'); -(abs(rVolt-double(SEVolt)))./(rVolt) +% assign(SEVolt,1); +% assign(SEVAngel,1); +% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; +% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0]; +% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1); +% solvesdp([],Objective,options) +% double(Objective) +% double(SEVolt) +% double(SEVAngel) +%% 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); +[x,fval,exitflag,info] = solve(Opt,x0); +info +x +%% 输出结果 +% fprintf('相对误差\n'); +% (abs(rVolt-double(SEVolt)))./(rVolt) +% MaxDeviation(rVolt,SEVolt,rVAngel,SEVAngel) diff --git a/run.m b/run.m index f860291..a7b264d 100644 --- a/run.m +++ b/run.m @@ -1,3 +1,6 @@ +clear +clc +yalmip('clear') addpath('.\Powerflow') [~, ~, ~, ~,Volt,Vangle,Y,Yangle,r,c,newwordParameter,PG,QG,PD,QD,Balance]=pf('ieee14.dat', '0'); %% 量测量 @@ -12,6 +15,7 @@ sigma=0.05;% %电压幅值 rVolt=Volt; %幅值 mVolt=rVolt.*(normrnd(0,sigma,length(Volt),1)+1);%电压量测量 +rVAngel=Vangle; %% 电流 %注入电流 cmpY=Y.*exp(1j*sparse(r,c,Yangle,length(Y),length(Y)));%复数导纳矩阵 @@ -52,37 +56,55 @@ stateVarCount=2*length(Volt); measurements=length(mVolt)+length(mI)+length(mBranchI)+length(mBranchP)+length(mBranchQ)+length(mPG)+length(mQG); fprintf('冗余度 %f\n',measurements/stateVarCount); %% 进入状态估计计算 -SEVolt=sdpvar(length(Volt),1); -SEVAngel=sdpvar(length(Vangle),1); -Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%%电压 +% SEVolt=sdpvar(length(Volt),1,'full'); +% SEVAngel=sdpvar(length(Vangle),1,'full'); +% Objective=(SEVolt-mVolt)'*(1./sigma^2*eye(length(mVolt)))*(SEVolt-mVolt);%电压 %% 支路电流 -cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 -cmpSEBranchI=(cmpV(lineI)-cmpV(lineJ))./(lineR+1j*lineX);%复数支路电流 -SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 -Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 +% cmpSEV=SEVolt.*exp(1j*SEVAngel); %复数电压 +% cmpSEBranchI=(cmpSEV(lineI)-cmpSEV(lineJ))./(lineR+1j*lineX);%复数支路电流 +% SEBranchI=abs(cmpSEBranchI);% 支路电流幅值 +% Objective=Objective+(SEBranchI-mBranchI)'*(1./sigma^2*eye(length(mBranchI)))*(SEBranchI-mBranchI);%%电流 %% 支路功率 -SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); -Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); -Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); +% SEBranchP=real((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); +% SEBranchQ=imag((cmpSEV(lineI)-cmpSEV(lineJ)).*conj(SEBranchI)); +% Objective=Objective+(SEBranchP-mBranchP)'*(1./sigma^2*eye(length(mBranchP)))*(SEBranchP-mBranchP); +% Objective=Objective+(SEBranchQ-mBranchQ)'*(1./sigma^2*eye(length(mBranchQ)))*(SEBranchQ-mBranchQ); %% 0注入节点 zerosInjectionIndex=1:length(Volt); zerosInjectionIndex=zerosInjectionIndex( ~(PD~=0|QD~=0|PG~=0|QG~=0) ); -PQ=diag(SEVolt)*conj(cmpY*SEVolt); -zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 -zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 +% PQ=diag(cmpSEV)*cmpY*cmpSEV; +% zeroInjP=real(PQ(zerosInjectionIndex));%% 0注入节点 +% zeroInjQ=imag(PQ(zerosInjectionIndex));%% 0注入节点 %% 发电机注入功率 % 先找到只有发电机的节点 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))); +% 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求解 -Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; -% Constraints=[[zeros(length(c)) A' -eye(length(lbounds))]*x==-c;[A zeros(length(b)) zeros(length(b),length(lbounds))]*x<=b;10>=x>=0]; -options = sdpsettings('verbose',2,'solver','ipopt'); -solvesdp(Constraints,Objective,options) -double(Objective) +% assign(SEVolt,1); +% assign(SEVAngel,1); +% Constraints=[SEVAngel(Balance)==0,zeroInjP==0,zeroInjP==0]; +% Constraints=[SEVAngel(Balance)==0,PQ(zerosInjectionIndex)==0]; +% options = sdpsettings('verbose',2,'solver','filtersd','usex0',1); +% solvesdp([],Objective,options) +% double(Objective) +% double(SEVolt) +% double(SEVAngel) +%% Opti ToolBox +Busnum=length(Volt); +seOpti=SEOpti(); +seOpti=seOpti.init(Busnum, mVolt,sigma,newwordParameter,zerosInjectionIndex,cmpY ,onlyPG,onlyQG,PG,QG,Balance); +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); +[x,fval,exitflag,info] = solve(Opt,x0); +%% 输出结果 +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)