stateestimateionyalmip-lu9-forcomparision/Run_YALMIP.asv

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('ieee30.dat');

%% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ
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';
%% <20><>ֵ-<2D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
PG0=PG;
QG0=QG;
PD0=PD;
QD0=QD;
PDReal=PD;%<25><>ֵ
QDReal=QD;%<25><>ֵ
%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);
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,PD,PD0,QD,randPDind,Loadi,notLoadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);

%% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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');
%% Ŀ<>꺯<EFBFBD><EABAAF>
Objective=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi);
%AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
%% <20><><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD><EFBFBD>Լӿ<D4BC><D3BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ٶȡ<D9B6>
assign(Volt(:),1);
assign(UAngel(:),0);
assign(PD(:),PD0(:));
assign(QD(:),QD0(:));
% assign(PG(:),PG0(:));
% assign(QG(:),QG0(:));
%% YALMIP<49><50><EFBFBD><EFBFBD>
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;
Constraints = [%AngleIJ-sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum)==0, ...
    dP==0, ...
    dQ==0, ...
     PG(setxor(1:Busnum, PGi) )==0, ...
     QG(setxor(1:Busnum, PVi) )==0, ...
     PD(PD0==0)==0, ...
     QD(QD0==0)==0, ...
    0.9*ones(Busnum,1)<=Volt<=1.1*ones(Busnum,1), ...
    -60*ones(Busnum,1)<=PG<=60*ones(Busnum,1), ...
    -60*ones(Busnum,1)<=QG<=60*ones(Busnum,1)
    ];
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)
    double(PG)-PG0;
    sol
else
    display('Hmm, something went wrong!');
    sol.info
    sol.solveroutput
    yalmiperror(sol.problem)
end

toc