facat
/
stateestimateionyalmip-lu9-forcomparision
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

代码写好了，正在调试 

Signed-off-by: facat <dugg@21cn.com>
This commit is contained in:
facat 2013-04-17 10:04:43 +08:00
parent 19975b1f88
commit 45f7dc7dff
11 changed files with 251 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
@Opti/Geteb.m Normal file
View File

@ -0,0 +1,7 @@
function [ output_args ] = Geteb( this )
%BE Summary of this function goes here
% Detailed explanation goes here
output_args=this.eb;
end

7
@Opti/Getgle.m Normal file
View File

@ -0,0 +1,7 @@
function [ output_args ] = Getgle( this )
%GLE Summary of this function goes here
% Detailed explanation goes here
output_args=this.gle;
end

31
@Opti/Opti.m Normal file
View File

@ -0,0 +1,31 @@
classdef Opti
%OPTI Summary of this class goes here
% Detailed explanation goes here
properties
PD0=NaN;
QD0=NaN;
Volt0=NaN;
wPD=NaN;
wQD=NaN;
wVolt=NaN;
PDi=NaN;
QDi=NaN;
Y=NaN;
Angel=NaN;
r=NaN;
c=NaN;
PG=NaN;
QG=NaN;
gle=NaN;%opti toolbox大于小于等于
eb=NaN;%等式、不等式约束右边的值
Balance=NaN;
end
methods
function this=Opti()
end
end
end

27
@Opti/equ.asv Normal file
View File

@ -0,0 +1,27 @@
function [ output_args,this ] = equ( this )
%EQU Summary of this function goes here
% Detailed explanation goes here
PDi=this.PDi;
QDi=this.QDi;
Volt0=this.Volt0;
Angel=this.Angel;
r=this.r;
c=this.c;
PG=this.PG;
QG=this.QG;
PD=x(1:length(PDi));
QD=x(length(PDi)+1:length(PDi)+length(QDi));
Volt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0));
VAngel=x(length(PDi)+length(QDi)+length(Volt0):end);
PD_=zeros(PD,length(Volt0));
QD_=zeros(QD,length(Volt0));
PD_(PDi)=PD;
QD_(QDi)=QD;
busNum=length(Volt0);
VMatrix=sparse(r,c,VAngel(r)-VAngel(c)-Angel,busNum,busNum);
dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt;
dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt;
output_args=[dP;dQ];
this.gle=ze
end

30
@Opti/equ.m Normal file
View File

@ -0,0 +1,30 @@
function [ output_args,this ] = equ( this,x )
%EQU Summary of this function goes here
% Detailed explanation goes here
PDi=this.PDi;
QDi=this.QDi;
Volt0=this.Volt0;
Angel=this.Angel;
r=this.r;
c=this.c;
PG=this.PG;
QG=this.QG;
Balance=this.Balance;
Y=this.Y;
PD=x(1:length(PDi));
QD=x(length(PDi)+1:length(PDi)+length(QDi));
Volt=x(length(PDi)+length(QDi)+1:length(PDi)+length(QDi)+length(Volt0));
VAngel=x(length(PDi)+length(QDi)+length(Volt0):end);
PD_=zeros(length(Volt0),1);
QD_=zeros(length(Volt0),1);
PD_(PDi)=PD;
QD_(QDi)=QD;
busNum=length(Volt0);
VMatrix=sparse(r,c,VAngel(r)-VAngel(c)-Angel,busNum,busNum);
dP=PG-PD_-diag(Volt)*(Y.*cos(VMatrix))*Volt;
dQ=QG-QD_-diag(Volt)*(Y.*spfun(@sin,VMatrix))*Volt;
output_args=[dP;dQ];
this.gle=zeros(length(output_args),1);
this.eb=this.gle;
end

19
@Opti/init.m Normal file
View File

@ -0,0 +1,19 @@
function [ this ] = init( this,Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angel,r,c,PG,QG,Balance )
%INIT Summary of this function goes here
% Detailed explanation goes here
this.PDi=PDi;
this.QDi=QDi;
this.wPD=wPD;
this.wQD=wQD;
this.PD0=PD0;
this.QD0=QD0;
this.Y=Y;
this.Angel=Angel;
this.r=r;
this.c=c;
this.PG=PG;
this.QG=QG;
this.Balance=Balance;
this.Volt0=Volt0;
end

15
@Opti/obj.asv Normal file
View File

@ -0,0 +1,15 @@
function [ output_args ] = obj( this )
%OBJ Summary of this function goes here
% Detailed explanation goes here
wPD=this.wPD;
wQD=this.wQD;
PD=this.PD;
QD=this.QD;
wVolt=t
t4=wPD.*(PD(PD0~=0)-PD0(PD0~=0)).^2;
t5=wQD.*(QD(QD0~=0)-QD0(QD0~=0)).^2;
t6=wVolt.*(Volt-Volt0').^2;
out_arg=sum(t4)+sum(t5)+sum(t6);
end

17
@Opti/obj.m Normal file
View File

@ -0,0 +1,17 @@
function [ output_args ] = obj( this,x )
%OBJ Summary of this function goes here
% Detailed explanation goes here
PDi=this.PDi;
QDi=this.QDi;
PD=x(1:length(PDi));
QD=x(length(PDi)+1:length(PDi)+length(QDi));
wPD=this.wPD;
wQD=this.wQD;
% wVolt=this.wVolt;
PD0=this.PD0;
QD0=this.QD0;
t4=(wPD(PDi).*(PD-PD0(PDi))).^2;
t5=(wQD(QDi).*(QD-QD0(QDi))).^2;
output_args=sum(t4)+sum(t5);
end

6
OPF_Init.m
View File

@ -33,9 +33,9 @@ wQG=1*ones(size(PVi,1),1);
%randInt=randperm(size(Loadi,1)); %randInt=randperm(size(Loadi,1));
%randPDind=randInt(1:10); %randPDind=randInt(1:10);
randPDind=0; randPDind=0;
wPD=1/0.05^2*ones(sum(PD~=0),1); wPD=1/0.05*ones(length(PD),1);
wQD=1/0.05^2*ones(sum(QD~=0),1); wQD=1/0.05*ones(length(QD),1);
wVolt=1/0.05^2*ones(Busnum,1); wVolt=1/0.05*ones(Busnum,1);
%wD(randPDind)=0;%一些负荷不约束 %wD(randPDind)=0;%一些负荷不约束
%wD(7)=0; %wD(7)=0;
% wD(11)=0; % wD(11)=0;

107
Run_YALMIP.asv
View File

@ -39,52 +39,77 @@ BalVolt=Volt(Balance);
% UAngel=sdpvar(Busnum,1); % UAngel=sdpvar(Busnum,1);
% PG=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1);
% QG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1);
PD=sdpvar(Busnum,1); % PD=sdpvar(Busnum,1);
QD=sdpvar(Busnum,1); % QD=sdpvar(Busnum,1);
AngleIJ=sdpvar(Busnum,Busnum,'full'); AngleIJ=sdpvar(Busnum,Busnum,'full');
%% 加噪音 %% 加噪音
PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
%% 目标函数 %% 目标函数
Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi); % Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,Volt,Volt0,wPG,wQG,wPD,wQD,wVolt,Loadi);
%AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); %AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
%% 赋初值，可以加快求解速度。 %% 赋初值，可以加快求解速度。
assign(Volt(:),1); % assign(Volt(:),1);
assign(UAngel(:),0); % assign(UAngel(:),0);
assign(PD(:),PD0(:)); % assign(PD(:),PD0(:));
assign(QD(:),QD0(:)); % assign(QD(:),QD0(:));
%% YALMIP²¿·Ö % %% YALMIP²¿·Ö
dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; % dP=PG0-PD-diag(Volt)*Y.*cos( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt;
dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt; % dQ=QG0-QD-diag(Volt)*Y.*sin( sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum) )*Volt;
Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0; % Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0;
Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ... % Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ...
dP(setdiff(1:Busnum,Loadi))==0, ... % dP(setdiff(1:Busnum,Loadi))==0, ...
dQ(setdiff(1:Busnum,Loadi))==0, ... % dQ(setdiff(1:Busnum,Loadi))==0, ...
% dP==0, ... % % dP==0, ...
% dQ==0, ... % % dQ==0, ...
PD(PD0==0)==0, ... % PD(PD0==0)==0, ...
QD(QD0==0)==0, ... % QD(QD0==0)==0, ...
0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ... % 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ...
Volt(Balance)==BalVolt, ... % Volt(Balance)==BalVolt, ...
UAngel(Balance)==0, ... % UAngel(Balance)==0, ...
0.8*PD0<=PD<=1.2*PD0; % 0.8*PD0<=PD<=1.2*PD0;
0.8*QD0<=QD<=1.2*QD0; % 0.8*QD0<=QD<=1.2*QD0;
]; % ];
options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1'); % options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1');
sol = solvesdp(Constraints,Objective,options); % sol = solvesdp(Constraints,Objective,options);
if sol.problem == 0 % if sol.problem == 0
fprintf('Volt\n'); % fprintf('Volt\n');
dvolt=double(Volt) % dvolt=double(Volt)
fprintf('VoltAngle\n'); % fprintf('VoltAngle\n');
dVangle=double(UAngel) % dVangle=double(UAngel)
fprintf('ojb\n'); % fprintf('ojb\n');
optimalObj=double(Objective) % optimalObj=double(Objective)
sol % sol
else % else
display('Hmm, something went wrong!'); % display('Hmm, something went wrong!');
sol.info % sol.info
sol.solveroutput % sol.solveroutput
yalmiperror(sol.problem) % yalmiperror(sol.problem)
end % end
%% Opti Toolbox
Busnum=length(Volt);
PDi=find(PD~=0);
QDi=find(QD~=0);
% PD0=PD0(PDi);
% QD0=QD0(QDi);
seOpti=Opti();
seOpti=seOpti.init(PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle);
opts = optiset('solver','filtersd');
opts.maxiter=85500;
opts.maxtime=3000;
opts.maxfeval=85000;
opts.maxnodes=85000;
opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
x0=zeros(length(PDi)+length(QDi),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.obj,'ndec',length(PDi)+length(QDi),'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0);
info
fprintf('Ä¿±êº¯Êý: %f\n',fval);
toc toc

32
Run_YALMIP.m
View File

@ -39,8 +39,8 @@ BalVolt=Volt(Balance);
% UAngel=sdpvar(Busnum,1); % UAngel=sdpvar(Busnum,1);
% PG=sdpvar(Busnum,1); % PG=sdpvar(Busnum,1);
% QG=sdpvar(Busnum,1); % QG=sdpvar(Busnum,1);
PD=sdpvar(Busnum,1); % PD=sdpvar(Busnum,1);
QD=sdpvar(Busnum,1); % QD=sdpvar(Busnum,1);
AngleIJ=sdpvar(Busnum,Busnum,'full'); AngleIJ=sdpvar(Busnum,Busnum,'full');
%% ¼ÓÔëÒô %% ¼ÓÔëÒô
PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0; PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
@ -86,5 +86,31 @@ QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
% sol.solveroutput % sol.solveroutput
% yalmiperror(sol.problem) % yalmiperror(sol.problem)
% end % end
%% Opti Toolbox
Busnum=length(Volt);
PDi=find(PD~=0);
QDi=find(QD~=0);
% PD0=PD0(PDi);
% QD0=QD0(QDi);
seOpti=Opti();
seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance);
opts = optiset('solver','filtersd');
opts.maxiter=85500;
opts.maxtime=3000;
opts.maxfeval=85000;
opts.maxnodes=85000;
opts.tolrfun=1e-4;
opts.tolafun=1e-4;
opts.warnings='all';
opts.display='off';
x0=zeros(length(PDi)+length(QDi)+length(Volt)*2,1);
% [~,seOpti]=seOpti.equ(x0);
[~,seOpti]=seOpti.equ(x0);
nlrhs=seOpti.Geteb();
nle=seOpti.Getgle();
Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmix',@seOpti.equ,nlrhs,nle,'options',opts)
% Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0);
info
fprintf('Ä¿±êº¯Êý: %f\n',fval);
toc toc