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

保持PG，QG不变。收敛了。

Browse Source 
Signed-off-by: facat <dugg@21cn.com>
This commit is contained in:
facat
2013-04-17 11:55:04 +08:00
parent fe0baaa029
commit 825696caa0
6 changed files with 29 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
@Opti/obj.m
View File

@@ -13,5 +13,6 @@ QD0=this.QD0;
t4=(wPD(PDi).*(PD-PD0(PDi))).^2; t4=(wPD(PDi).*(PD-PD0(PDi))).^2;
t5=(wQD(QDi).*(QD-QD0(QDi))).^2; t5=(wQD(QDi).*(QD-QD0(QDi))).^2;
output_args=sum(t4)+sum(t5); output_args=sum(t4)+sum(t5);
output_args=full(output_args);
end end

4
OPF_Init.m
View File

@@ -26,8 +26,8 @@ tPU=sparse(GenU(:,2));%
tQU=sparse(PVQU(:,1));% 无功上限 tQU=sparse(PVQU(:,1));% 无功上限
tPL=sparse(GenL(:,2));% 发电机有功下限 tPL=sparse(GenL(:,2));% 发电机有功下限
tQL=sparse(PVQL(:,1));% 无功下限 tQL=sparse(PVQL(:,1));% 无功下限
PG(PGi)=(tPU+tPL)/2; % PG(PGi)=(tPU+tPL)/2;
QG(PVi)=(tQU+tQL)/2; % QG(PVi)=(tQU+tQL)/2;
wPG=1*ones(size(PGi,1),1); wPG=1*ones(size(PGi,1),1);
wQG=1*ones(size(PVi,1),1); wQG=1*ones(size(PVi,1),1);
%randInt=randperm(size(Loadi,1)); %randInt=randperm(size(Loadi,1));

BIN
PD0.mat Normal file
View File

Binary file not shown.

BIN
QD0.mat Normal file
View File

Binary file not shown.

26
Run_YALMIP.asv
View File

@@ -5,7 +5,7 @@ tic
[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB, ... [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB, ...
Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL, ... Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL, ...
Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ... Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ...
pf('E:/算例/柳金Ⅰ926_21671693_2012-09-06/newFIle20.txt'); pf('E:/算例/东际911_2751267_2012-09-05/newFIle20.txt');
%% 潮流等式 %% 潮流等式
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
@@ -32,6 +32,10 @@ Volt0=Volt;
% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt'; % dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
% dQD=abs(dQ./QD); % dQD=abs(dQ./QD);
% maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]); % maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]);
xVolt=Volt;
xUAngel=UAngel;
% VMatrix=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
% dP=PG-PD-diag(xVolt)*(Y.*cos(VMatrix))*xVolt';
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD); [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
%% 定义变量 %% 定义变量
BalVolt=Volt(Balance); BalVolt=Volt(Balance);
@@ -44,6 +48,7 @@ BalVolt=Volt(Balance);
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;
save('PD0','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);
@@ -93,7 +98,7 @@ QDi=find(QD~=0);
% PD0=PD0(PDi); % PD0=PD0(PDi);
% QD0=QD0(QDi); % QD0=QD0(QDi);
seOpti=Opti(); seOpti=Opti();
seOpti=seOpti.init(PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle); seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance);
opts = optiset('solver','filtersd'); opts = optiset('solver','filtersd');
opts.maxiter=85500; opts.maxiter=85500;
opts.maxtime=3000; opts.maxtime=3000;
@@ -102,13 +107,16 @@ opts.maxnodes=85000;
opts.tolrfun=1e-4; opts.tolrfun=1e-4;
opts.tolafun=1e-4; opts.tolafun=1e-4;
opts.warnings='all'; opts.warnings='all';
opts.display='off'; opts.display='iter';
x0=zeros(length(PDi)+length(QDi),1); % x0=[zeros(length(PDi)+length(QDi),1); ...
% [~,seOpti]=seOpti.equ(x0); % ones(length(Volt),1); ...
% nlrhs=seOpti.nlrhs(); % zeros(length(Volt),1)];
% nle=seOpti.nle(); x0=[PD(PDi);QD(QDi);xVolt';xUAngel'];
% Opt = opti('fun',@seOpti.fun,'ndec',length(Volt)*2,'nlmix',@seOpti.equ,nlrhs,nle,'options',opts) [~,seOpti]=seOpti.equ(x0);
Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts) 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); [x,fval,exitflag,info] = solve(Opt,x0);
info info
fprintf('目标函数: %f\n',fval); fprintf('目标函数: %f\n',fval);

12
Run_YALMIP.m
View File

@@ -47,8 +47,12 @@ BalVolt=Volt(Balance);
% 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; load('PD0');
QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0; load('QD0');
% PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
% save('PD0','PD0');
% QD0=(1+normrnd(0,0.05,length(QD0),1)).*QD0;
% save('QD0','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);
@@ -98,7 +102,7 @@ QDi=find(QD~=0);
% QD0=QD0(QDi); % QD0=QD0(QDi);
seOpti=Opti(); seOpti=Opti();
seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance); seOpti=seOpti.init(Volt0,PDi,QDi,wPD,wQD,PD0,QD0,Y,Angle,r,c,PG,QG,Balance);
opts = optiset('solver','filtersd'); opts = optiset('solver','ipopt');
opts.maxiter=85500; opts.maxiter=85500;
opts.maxtime=3000; opts.maxtime=3000;
opts.maxfeval=85000; opts.maxfeval=85000;
@@ -118,5 +122,7 @@ Opt = opti('fun',@seOpti.obj,'ndec',length(Volt)*2+length(PDi)+length(QDi),'nlmi
% Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts) % Opt = opti('fun',@seOpti.obj,'ndec',length(PDi)+length(QDi),'options',opts)
[x,fval,exitflag,info] = solve(Opt,x0); [x,fval,exitflag,info] = solve(Opt,x0);
info info
fval=seOpti.obj(x);
fprintf('目标函数: %f\n',fval); fprintf('目标函数: %f\n',fval);
toc toc