189 lines
7.8 KiB
Matlab
189 lines
7.8 KiB
Matlab
function [totalLoss,JSM]=OPF(fileName)
|
||
tic
|
||
% clc
|
||
% clear
|
||
%% 存在问题
|
||
% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
|
||
%%
|
||
thesis=ForThesis(1,62);
|
||
for II=1:1
|
||
[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,Branchi,Branchg,Branchb,Transfork0]= ...
|
||
pf(fileName);
|
||
% pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
|
||
%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
|
||
%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
|
||
%pf('c:/file31.txt');
|
||
%% 计算功率因数
|
||
Loadi=QD~=0 | PD~=0;
|
||
PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
|
||
%%
|
||
Volt;
|
||
UAngel*180/3.1415926;
|
||
%% 通过潮流计算PG
|
||
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';
|
||
loss=sum(PGBal)-sum(PD);
|
||
%% 初值-即测量值
|
||
PG0=PG;
|
||
QG0=QG;
|
||
PD0=PD;
|
||
QD0=QD;
|
||
Volt0=Volt;
|
||
UAngel0=UAngel;
|
||
%%
|
||
PG0(Balance)=PGBal(Balance);
|
||
PG(Balance)=PGBal(Balance);
|
||
QG0(Balance)=QGBal(Balance);
|
||
QG0(PVi)=QGBal(PVi);
|
||
QG(PVi)=QGBal(PVi);
|
||
%% 真实值
|
||
RealPG=PG0;
|
||
RealQG=QG0;
|
||
RealPD=PD0;
|
||
RealQD=QD0;
|
||
%%
|
||
%PGi=zeros(1,1);
|
||
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
|
||
Gap=(Init_L*Init_Z'-Init_U*Init_W');
|
||
KK=0;
|
||
plotGap=zeros(1,60);
|
||
ContrlCount=size(Loadi,1)*2+Busnum*2;
|
||
kmax=100;
|
||
%% 20120523 临时
|
||
QD_NON_ZERO=QD(PD==0 & QD~=0);
|
||
QD_NON_ZERO_IND=find(PD==0 & QD~=0);
|
||
%%
|
||
Precision=Precision/1;
|
||
%% 加误差
|
||
% PD0(Loadi)=PD0(Loadi).*(1+0.02);
|
||
% QD0(Loadi)=QD0(Loadi).*(1+0.02);
|
||
|
||
PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
|
||
QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
|
||
% save('20PD0.mat','PD0');
|
||
% save('20QD0.mat','QD0');
|
||
% load('20PD0.mat');
|
||
% load('20QD0.mat');
|
||
% PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
|
||
% QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
|
||
setIneq=0;
|
||
|
||
%% 读变压器容量
|
||
noDataTransCapacity=0;
|
||
while(abs(Gap)>Precision)
|
||
if KK>kmax
|
||
break;
|
||
end
|
||
plotGap(KK+1)=Gap;
|
||
Init_u=Gap/2/RestraintCount*CenterA;
|
||
AngleIJMat=0;
|
||
%% 开始计算OPF
|
||
%% 形成等式约束的雅克比
|
||
deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
|
||
%% 形成不等式约束的雅克比
|
||
deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD);
|
||
%%
|
||
L_1Z=diag(Init_Z./Init_L);
|
||
U_1W=diag(Init_W./Init_U);
|
||
%% 形成海森阵
|
||
deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
|
||
%% 形成ddHy
|
||
ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
|
||
%% 开始构建ddg
|
||
ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
|
||
%% 开始构建deltF
|
||
deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
|
||
%% 形成方程矩阵
|
||
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
|
||
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
|
||
Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi);
|
||
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi);
|
||
Ly=Mat_H;
|
||
if setIneq==0
|
||
PDL=RealPD(Loadi);
|
||
% PDL(PDL>0)=0.800*PDL(PDL>0);
|
||
% PDL(PDL<0)=1.200*PDL(PDL<0);
|
||
PDL(PDL>0)=0.600*PDL(PDL>0).*unifrnd(0.75,1.25,sum(PDL>0),1);
|
||
PDL(PDL<0)=1.400*PDL(PDL<0).*unifrnd(0.75,1.25,sum(PDL<0),1);
|
||
PDL(PDL==0)=-0.400;
|
||
QDL=RealQD(Loadi);
|
||
% QDL(QDL>0)=0.800*QDL(QDL>0);
|
||
% QDL(QDL<0)=1.200*QDL(QDL<0);
|
||
QDL(QDL>0)=0.600*QDL(QDL>0).*unifrnd(0.75,1.25,sum(QDL>0),1);
|
||
QDL(QDL<0)=1.400*QDL(QDL<0).*unifrnd(0.75,1.25,sum(QDL<0),1);
|
||
QDL(QDL==0)=-0.400;
|
||
end
|
||
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,PDL,QDL);
|
||
|
||
if setIneq==0
|
||
|
||
PDU=RealPD(Loadi);
|
||
% PDU(PDU>0)=1.200*PDU(PDU>0);
|
||
% PDU(PDU<0)=0.800*PDU(PDU<0);
|
||
PDU(PDU>0)=1.400*PDU(PDU>0).*unifrnd(0.75,1.25,sum(PDU>0),1);
|
||
PDU(PDU<0)=0.600*PDU(PDU<0).*unifrnd(0.75,1.25,sum(PDU<0),1);
|
||
PDU(PDU==0)=0.400;
|
||
QDU=RealQD(Loadi);
|
||
% QDU(QDU>0)=1.200*QDU(QDU>0);
|
||
% QDU(QDU<0)=0.800*QDU(QDU<0);
|
||
QDU(QDU>0)=1.400*QDU(QDU>0).*unifrnd(0.75,1.25,sum(QDU>0),1);
|
||
QDU(QDU<0)=0.600*QDU(QDU<0).*unifrnd(0.75,1.25,sum(QDU<0),1);
|
||
QDU(QDU==0)=0.400;
|
||
|
||
end
|
||
setIneq=1;
|
||
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity,PDU,QDU);
|
||
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
|
||
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
|
||
%% 开始解方程
|
||
fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
|
||
XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
|
||
%%取各分量
|
||
[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
|
||
[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi);
|
||
Gap=(Init_L*Init_Z'-Init_U*Init_W');
|
||
KK=KK+1;
|
||
end
|
||
fprintf('迭代次数%d\n',KK);
|
||
fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)));
|
||
% DrawGap(plotGap);
|
||
%%
|
||
%Volt=full(Volt');
|
||
%PD=full(PD);
|
||
%% 统计PD误差
|
||
% absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) );
|
||
absPDLoad=abs( (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi) );
|
||
maxPDError=max(absPDLoad(absPDLoad<10))
|
||
Loadi(maxPDError==absPDLoad)
|
||
absQDLoad=abs( (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi) );
|
||
maxQDError=max(absQDLoad(absQDLoad<10))
|
||
Loadi(maxQDError==absQDLoad)
|
||
disp('index');
|
||
%Loadi(absPDLoad==maxPDError);
|
||
%% 计算总线损
|
||
totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
|
||
fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
|
||
fprintf('线损率为 %f\n',full(totalLoss/sum(PG)));
|
||
%% 计算各线损
|
||
% Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
|
||
% Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,k0,Volt,Angle)
|
||
thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
|
||
|
||
end
|
||
PD(Loadi)=thesis.MeanPD();
|
||
QD(Loadi)=thesis.MeanQD();
|
||
PG(Balance)=thesis.MeanPG();
|
||
QG(PVi)=thesis.MeanQG();
|
||
thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))%JSN
|
||
JSM=thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))%JSM
|
||
thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))%JZN
|
||
thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))%JZM
|
||
% thesis.SqareDeviation(RealPG(Balance(1)),RealQG(PVi(1)),RealPD(Loadi),RealQD(Loadi))
|
||
thesis.SqareDeviation(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
|
||
thesis.PercentOfPass(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
|
||
thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\算例\新民Ⅰ906_2729823_2012-09-06\newFIle20.txt',PD0,QD0)
|
||
thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
|
||
toc
|
||
end
|