把PQQG作为常数加入潮流方程中。

添加了很多统计函数。
This commit is contained in:
facat 2013-01-23 21:27:33 +08:00
parent e96ebb226d
commit 4944e38cbd
17 changed files with 256 additions and 176 deletions

View File

@ -9,6 +9,6 @@ this.currentPos=this.currentPos+1;
this.PDArray(:,this.currentPos)=PD;
this.QDArray(:,this.currentPos)=QD;
this.PGArray(this.currentPos)=PG;
% this.QGArray(this.currentPos)=QG;
this.QGArray(this.currentPos)=QG;
end

View File

@ -0,0 +1,10 @@
function [ output_args ] = MaxBranchDeviation( this, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
%% 最大支路功率偏差
% 支路功率包括线路和变压器
%%
[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchg,Branchb,k0,Volt0,Angle0);
[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchg,Branchb,k0,Volt,Angle);
real(dispLineloss0 - dispLineloss)./real(dispLineloss0)
end

View File

@ -0,0 +1,14 @@
function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
%%
% 线
%%
[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);
[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);
t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3);
t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3);
t11=dispLineloss0(:,3)>1e-5;%
t22=dispTransloss0(:,3)>1e-5;%
t3=abs([t1(t11);t2(t22)]);
output_args=max(t3(t3~=Inf));
end

View File

@ -6,6 +6,6 @@ QDMaxDev=max(abs( (this.QDArray-QD0Array)./QD0Array ),[],2);
PG0Array=repmat(PG0,this.sampleNum,1);
QG0Array=repmat(QG0,this.sampleNum,1);
PGMaxDev=max(abs( (PG0Array-this.PGArray)./PG0Array ));
% QGMaxDev=max( abs( (QG0Array-this.QGArray)./QG0Array ) );
output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev]);
QGMaxDev=max( abs( (QG0Array-this.QGArray)./QG0Array ) );
output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev;QGMaxDev]);
end

View File

@ -0,0 +1,12 @@
function [ output_args ] = PercentOfPass( this, PG0,QG0,PD0,QD0 )
%%
PD=this.PDArray(:,end);%
QD=this.QDArray(:,end);%
t1=PD-PD0;
t2=QD-QD0;
t3=abs(t1)/3*0.05;
t4=abs(t2)/3*0.05;
output_args=sum(t3<=1)+sum(t4<=1);
end

View File

@ -0,0 +1,13 @@
function [ output_args ] = SqareDeviation( this, PG0,QG0,PD0,QD0)
%% 0
PD=this.PDArray(:,end);%
QD=this.QDArray(:,end);%
t1=PD-PD0;
t2=t1./0.05;
t3=QD-QD0;
t4=t3./0.05;
output_args=sqrt(sum(t2.^2)+sum(t4.^2));
end

View File

@ -0,0 +1,19 @@
function [ output_args ] = StatBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel )
%STATBRANCHDEVIATION Summary of this function goes here
% Detailed explanation goes here
%% L2
% 线
%%
[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);
[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);
t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3);
t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3);
t11=dispLineloss0(:,3)>1e-5;%
t22=dispTransloss0(:,3)>1e-5;%
t3=[t1(t11);t2(t22)];
t4=t3.^2;
t5=sum(t4).^.5;
output_args=t5
end

View File

@ -6,15 +6,15 @@ PD0Array=repmat(PD0,1,this.sampleNum);
QD0Array=repmat(QD0,1,this.sampleNum);
PDDev=(this.PDArray-PD0Array)/0.05;
QDDev=(this.QDArray-QD0Array)/0.05;
PG0Array=repmat(PG0,this.sampleNum,1);
QG0Array=repmat(QG0,this.sampleNum,1);
PGDev=(PG0Array-this.PGArray)/0.01;
% PG0Array=repmat(PG0,this.sampleNum,1);
% QG0Array=repmat(QG0,this.sampleNum,1);
% PGDev=(PG0Array-this.PGArray)/0.01;
% QGDev=(QG0Array-this.QGArray)/0.01;
wholeMat=[PDDev;QDDev;PGDev'];
wholeMat=[PDDev;QDDev;];
t1=wholeMat.^2;
t2=sum(t1,1);
t3=t2/size(t1,1);
t4=t3.^2;
t4=t3.^.5;
output_args=sum(t4)/length(t4);
end

View File

@ -1,9 +1,7 @@
function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
KK=999;
%PU=GenU(:,2);%
PU=5*ones(length(GenU(:,2)),1);
%QU=PVQU(:,1);%
QU=5*ones(length(PVQU(:,1)),1);
VoltU=(1.1)*ones(1,Busnum);
%VoltU=10*ones(1,Busnum);
PDU=PD0(Loadi);

View File

@ -1,9 +1,6 @@
function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
KK=999;
%PL=GenL(:,2);%
PL=-5*ones(length(GenL(:,2)),1);
%QL=PVQL(:,1);%
QL=-5*ones(length(PVQL(:,1)),1);
VoltL=(0.9)*ones(1,Busnum);
%VoltL=-10*ones(1,Busnum);
PDL=PD0(Loadi);

View File

@ -1,8 +1,8 @@
%% 线
function Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,k0,Volt,Angle)
function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchi,Branchg,Branchb,k0,Volt,Angle)
%format long
fprintf('\n');
fprintf('\t\t MW\t MVar')
% fprintf('\n');
% fprintf('\t\t MW\t MVar')
cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle);
cmpVolt=cmpVolt';
y0=1i*Lineb2;
@ -16,13 +16,17 @@ deltLineS=Sij+Sji;
dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
%full(dispLineloss)
dispLineloss=sortrows(dispLineloss,-3);
full(dispLineloss)
full(dispLineloss);
%%
yij=1./(Transr+1i*Transx);
Sij=cmpVolt(Transi)./k0.*conj( ( cmpVolt(Transi)./k0- cmpVolt(Transj) ).*yij );
Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij );
deltTransS=Sij+Sji;
%%
%
deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))));
%%
dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100];
dispTransloss=sortrows(dispTransloss,-3);
full(dispTransloss)
full(dispTransloss);
end

243
OPF.asv
View File

@ -1,120 +1,137 @@
tic
clc
clear
[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('c:/newFIle.txt');
%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视9223-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';
%% 初值-即测量值
PG0=PG;
QG0=QG;
PD0=PD;
QD0=QD;
PDReal=PD;%真值
QDReal=QD;%真值
%PD0(12)=PD0(12)+0.001;
%%
PG0(Balance)=PGBal(Balance);
PG(Balance)=PGBal(Balance);
QG0(Balance)=QGBal(Balance);
QG0(PVi)=QGBal(PVi);
QG(PVi)=QGBal(PVi);
%%
[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,PD0,QD0,QD,PD);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0;
plotGap=zeros(1,60);
ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2;
kmax=60;
%% 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+normrnd(0,0.05,length(Loadi),1));
QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
%% 读变压器容量
%[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
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);
thesis=ForThesis(4,8);
for II=1:4
[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('c:/newFIle.txt');
%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视9223-1_0.5_120%.txt');
%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
%pf('c:/file31.txt');
%% 形成方程矩阵
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;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity);
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);
%% 计算功率因数
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';
%% 初值-即测量值
PG0=PG;
QG0=QG;
PD0=PD;
QD0=QD;
%%
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=KK+1;
KK=0;
plotGap=zeros(1,60);
ContrlCount=size(Loadi,1)*2+Busnum*2;
kmax=60;
%% 20120523 临时
QD_NON_ZERO=QD(PD==0 & QD~=0);
QD_NON_ZERO_IND=find(PD==0 & QD~=0);
%%
Precision=Precision/10;
%% 加误差
PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
% PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
% QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
%% 读变压器容量
%[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
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;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity);
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);
thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
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)-PDReal(Loadi))./PDReal(Loadi) );
maxPDError=max(absPDLoad(absPDLoad<10))
absQDLoad=abs( (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi) );
maxQDError=max(absQDLoad(absQDLoad<10))
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);
toc
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))
thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi));
toc

41
OPF.m
View File

@ -1,10 +1,13 @@
tic
clc
clear
%%
% 1 20130123
%%
thesis=ForThesis(4,8);
for II=1:4
[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('c:/newFIle.txt');
[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('E:\\911_2751267_2012-09-05\newFIle20.txt');
%pf('D:\Project\\\\\9223-1_0.5_120%.txt');
%pf('D:\Project\\\919.txt');
%pf('c:/file31.txt');
@ -25,19 +28,22 @@ for II=1:4
QG0=QG;
PD0=PD;
QD0=QD;
PDReal=PD;%ÕæÖµ
QDReal=QD;%ÕæÖµ
%PD0(12)=PD0(12)+0.001;
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;
%%
PVi=zeros(0,0);
PGi=zeros(0,0);
[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,PD0,QD0,QD,PD);
%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);
@ -88,8 +94,8 @@ for II=1:4
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;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity);
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%%
@ -109,10 +115,12 @@ for II=1:4
%PD=full(PD);
%% PD
% absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) );
absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) );
absPDLoad=abs( (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi) );
maxPDError=max(absPDLoad(absPDLoad<10))
absQDLoad=abs( (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi) );
Loadi(maxPDError==absPDLoad)
absQDLoad=abs( (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi) );
maxQDError=max(absQDLoad(absQDLoad<10))
Loadi(maxQDError==absQDLoad)
disp('index');
%Loadi(absPDLoad==maxPDError);
%% 线
@ -128,7 +136,12 @@ PD(Loadi)=thesis.MeanPD();
QD(Loadi)=thesis.MeanQD();
PG(Balance)=thesis.MeanPG();
QG(PVi)=thesis.MeanQG();
thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi));
thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
% 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)
thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
toc

View File

@ -15,19 +15,20 @@ tPU=sparse(GenU(:,2));%
tQU=sparse(PVQU(:,1));%
tPL=sparse(GenL(:,2));%
tQL=sparse(PVQL(:,1));%
PG(PGi)=(tPU+tPL)/2;
wPG=0*ones(size(PGi,1),1);
wQG=0*ones(size(PVi,1),0);
wPG=0;
wQG=0;
%randInt=randperm(size(Loadi,1));
%randPDind=randInt(1:10);
randPDind=0;
wPD=20*ones(size(Loadi,1),1);
wQD=0;
wPD=0/.05^2*ones(size(Loadi,1),1);
% wPD([2 7])=0;
wQD=0/.05^2*ones(size(Loadi,1),1);
% wQD([2 7])=0;
%wD(randPDind)=0;%
%wD(7)=0;
% wD(11)=0;
PD=1*PD0;
PD=0.5*PD0;
%powerFacter=0.98;
%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
QD=QD0;
QD=.5*QD0;
end

View File

@ -1,36 +1,20 @@
function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
%%
sizePGi=size(PGi,1);
sizePVi=size(PVi,1);
sizeLoadi=size(Loadi,1);
%%
dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
dg2_dPg=sparse(sizePGi,sizePVi);
dg3_dPg=sparse(sizePGi,sizeLoadi);
dg4_dPg=sparse(sizePGi,length(Loadi));
dg5_dPg=sparse(sizePGi,Busnum);
%%
dg1_dQr=sparse(sizePVi,sizePGi);
dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
dg3_dQr=sparse(sizePVi,sizeLoadi);
dg4_dQr=sparse(sizePVi,length(Loadi));
dg5_dQr=sparse(sizePVi,Busnum);
%%
dg1_dPD=sparse(size(Loadi,1),size(PGi,1));
dg2_dPD=sparse(size(Loadi,1),size(PVi,1));
dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg4_dPD=sparse(size(Loadi,1),length(Loadi));
dg5_dPD=sparse(size(Loadi,1),Busnum);
%%
dg1_dQD=sparse(size(Loadi,1),size(PGi,1));
dg2_dQD=sparse(size(Loadi,1),size(PVi,1));
dg3_dQD=sparse(length(Loadi),length(Loadi));
dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg5_dQD=sparse(size(Loadi,1),Busnum);
%%
dg1_dx=sparse(2*Busnum,sizePGi);
dg2_dx=sparse(2*Busnum,sizePVi);
dg3_dx=sparse(2*Busnum,sizeLoadi);
dg4_dx=sparse(2*Busnum,length(Loadi));
dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);

View File

@ -1,6 +1,4 @@
function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c); %

2
pf.m
View File

@ -1,4 +1,4 @@
function [kmax,Precision,Uangle,U,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(FileName)
function [kmax,Precision,Uangle,U,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)
%**************************************************************************
%
% -