改用Opti Toolbox

Signed-off-by: facat <dugg@21cn.com>
2013-04-14 19:40:35 +08:00 · 2013-04-14 19:40:35 +08:00 · 32f96a94f0
parent f50bbd3435
commit 32f96a94f0
9 changed files with 203 additions and 47 deletions
--- a/@SEOpti/SEOpti.m
+++ 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
+
--- a/@SEOpti/equ.m
+++ 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
+
--- a/@SEOpti/fun.m
+++ 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
+
--- a/@SEOpti/init.m
+++ 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
+
--- a/@SEOpti/nle.m
+++ 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
+
--- a/@SEOpti/nlrhs.m
+++ 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
+
--- a/MaxDeviation.m
+++ 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
+
--- 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)
--- 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)