stateestimateionyalmip-lu9-.../Run_YALMIP.m

91 lines
2.9 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

clc
clear
yalmip('clear')
tic
[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, ...
Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ...
pf('E:/算例/柳金926_21671693_2012-09-06/newFIle20.txt');
%% 潮流等式
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
%% 初值-即测量值
PG0=PG;
QG0=QG;
PD0=PD;
QD0=QD;
PDReal=PD;%真值
QDReal=QD;%真值
%PD0(12)=PD0(12)+0.001;
PG0(Balance)=PGBal(Balance);
QG0(Balance)=QGBal(Balance);
QG0(PVi)=QGBal(PVi);
PG(Balance)=PGBal(Balance);
QG(PVi)=QGBal(PVi);
Volt0=Volt;
% PF=1;
% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
% dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
% dPD=abs(dP./PD);
% dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
% dQD=abs(dQ./QD);
% maxdPQ=max([dPD(dPD<10);dQD(dQD<10)]);
[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);
% Volt=sdpvar(Busnum,1);
% UAngel=sdpvar(Busnum,1);
% PG=sdpvar(Busnum,1);
% QG=sdpvar(Busnum,1);
PD=sdpvar(Busnum,1);
QD=sdpvar(Busnum,1);
AngleIJ=sdpvar(Busnum,Busnum,'full');
%% 加噪音
PD0=(1+normrnd(0,0.05,length(PD0),1)).*PD0;
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);
%AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
%% 赋初值,可以加快求解速度。
% assign(Volt(:),1);
% assign(UAngel(:),0);
% assign(PD(:),PD0(:));
% assign(QD(:),QD0(:));
% %% YALMIP部分
% 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;
% Loadi=PD0~=0 | QD0~=0 |PG0~=0|QG0~=0;
% Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ...
% dP(setdiff(1:Busnum,Loadi))==0, ...
% dQ(setdiff(1:Busnum,Loadi))==0, ...
% % dP==0, ...
% % dQ==0, ...
% PD(PD0==0)==0, ...
% QD(QD0==0)==0, ...
% 0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ...
% Volt(Balance)==BalVolt, ...
% UAngel(Balance)==0, ...
% 0.8*PD0<=PD<=1.2*PD0;
% 0.8*QD0<=QD<=1.2*QD0;
% ];
% options = sdpsettings('verbose',2,'showprogress',1,'debug',0,'solver','ipopt','usex0','1');
% sol = solvesdp(Constraints,Objective,options);
% if sol.problem == 0
% fprintf('Volt\n');
% dvolt=double(Volt)
% fprintf('VoltAngle\n');
% dVangle=double(UAngel)
% fprintf('ojb\n');
% optimalObj=double(Objective)
% sol
% else
% display('Hmm, something went wrong!');
% sol.info
% sol.solveroutput
% yalmiperror(sol.problem)
% end
toc