facat
/
pes2014-debug2-laplace
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

加错误数据识别 

Signed-off-by: facat <dmy@dmy-PC.(none)>
This commit is contained in:
unknown 2012-07-25 15:08:54 +08:00
parent f9dc83b862
commit 798e1e88ac
8 changed files with 242 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
FormH.m
View File

@ -1,7 +1,7 @@
function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND) function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND)
%% %%
QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt'); QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
QD(QD~=0)=PD(QD~=0)./tan(QDcos); %QD(QD~=0)=PD(QD~=0)./tan(QDcos);
%QD(QD_NON_ZERO_IND)=QD_NON_ZERO; %QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
%% %%
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);

4
Modification.m
View File

@ -1,8 +1,8 @@
function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD]=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,Loadi) function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD]=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,Loadi)
AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU); AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
fprintf('AlphaP %f\n',full(AlphaP)); %fprintf('AlphaP %f\n',full(AlphaP));
AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW); AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
fprintf('AlphaD %f\n',full(AlphaD)); %fprintf('AlphaD %f\n',full(AlphaD));
Init_Z=Init_Z+AlphaD*deltZ'; Init_Z=Init_Z+AlphaD*deltZ';
Init_L=Init_L+AlphaP*deltL'; Init_L=Init_L+AlphaP*deltL';

18
OPF.m
View File

@ -1,9 +1,9 @@
tic tic
clear 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]= ... [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]= ...
pf('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919.txt'); pf('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\½ðºþ924£¨2-1£©_0.5_85%.txt');
%pf('c:/file31.txt'); %pf('c:/file31.txt');
%pf('ieee10471PG.dat'); %pf('ieee118PG.dat');
%% ¼ÆË㹦ÂÊÒòÊý %% ¼ÆË㹦ÂÊÒòÊý
atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0)) atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0))
@ -13,13 +13,16 @@ UAngel*180/3.1415926;
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
%% ³õÖµ %% ³õÖµ-¼´²âÁ¿Öµ
PG0=PG; PG0=PG;
PD0=PD; PD0=PD;
PDReal=PD;%ÕæÖµ
PD0(12)=PD0(12)+0.001;
%% %%
PG0(Balance)=PGBal(Balance); PG0(Balance)=PGBal(Balance);
%% %%
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD); [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
Gap=(Init_L*Init_Z'-Init_U*Init_W'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
plotGap=zeros(1,50); plotGap=zeros(1,50);
@ -78,6 +81,11 @@ DrawGap(plotGap);
%Volt=full(Volt'); %Volt=full(Volt');
%PD=full(PD); %PD=full(PD);
%% ͳ¼ÆPDÎó²î %% ͳ¼ÆPDÎó²î
maxPDError=max(abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) )); % absPDLoad=abs( (PD(Loadi)-PD0(Loadi))./PD0(Loadi) );
absPDLoad=abs( (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi) );
maxPDError=max(absPDLoad)
disp('index')
Loadi(absPDLoad==maxPDError)
toc toc

106
OPFWrongDataCheck.asv Normal file
View File

@ -0,0 +1,106 @@
clc
clear
for II=1:16
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]= ...
pf('D:\Project\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt');
%pf('c:/file31.txt');
%pf('ieee118PG.dat');
%% 计算功率因数
atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
Volt;
UAngel*180/3.1415926;
%% 通过潮流计算PG
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
%% 初值-即测量值
PG0=PG;
PD0=PD;
PDReal=PD;%真值
%%
PG0(Balance)=PGBal(Balance);
%%
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
PD0(Loadi(II))=PD0(Loadi(II))*(1+0.018);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0;
plotGap=zeros(1,50);
ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
kmax=60;
%% 20120523 临时
QD_NON_ZERO=QD(PD==0 & QD~=0);
QD_NON_ZERO_IND=find(PD==0 & QD~=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);
%%
L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U);
%% 形成海森阵
deltdeltF=func_deltdeltF(PVi,wG,wD,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);
%% 开始构建deltF
deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,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,Loadi);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%% 开始解方程
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]=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,Loadi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;
end
fprintf('迭代次数%d\n',KK);
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,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);
%disp('index')
LoadiArray=Loadi(absPDLoad==maxPDError);
if length(LoadiArray)>1
disp('没找出')
Loadi(II)
end
if length(LoadiArray)==1
if LoadiArray~=Loadi(II)
disp('没找出')
Loadi(II)
end
end
toc;
end

106
OPFWrongDataCheck.m Normal file
View File

@ -0,0 +1,106 @@
clc
clear
for II=1:53
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]= ...
pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
%pf('c:/file31.txt');
%pf('ieee118PG.dat');
%% 计算功率因数
atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
Volt;
UAngel*180/3.1415926;
%% 通过潮流计算PG
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
%% 初值-即测量值
PG0=PG;
PD0=PD;
PDReal=PD;%真值
%%
PG0(Balance)=PGBal(Balance);
%%
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
PD0(Loadi(II))=PD0(Loadi(II))*(1+0.086);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0;
plotGap=zeros(1,50);
ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
kmax=60;
%% 20120523 临时
QD_NON_ZERO=QD(PD==0 & QD~=0);
QD_NON_ZERO_IND=find(PD==0 & QD~=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);
%%
L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U);
%% 形成海森阵
deltdeltF=func_deltdeltF(PVi,wG,wD,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);
%% 开始构建deltF
deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,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,Loadi);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
Ly=Mat_H;
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
%% 开始解方程
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]=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,Loadi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;
end
fprintf('迭代次数%d\n',KK);
ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,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);
%disp('index')
LoadiArray=Loadi(absPDLoad==maxPDError);
if length(LoadiArray)>1
disp('没找出')
Loadi(II)
end
if length(LoadiArray)==1
if LoadiArray~=Loadi(II)
disp('没找出')
Loadi(II)
end
end
toc;
end

12
OPF_Init.m
View File

@ -1,5 +1,5 @@
function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,HH) function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD)
Loadi=find(QD~=0); Loadi=find(QD~=0 & PD~=0);
%Loadi=[1:Busnum]'; %Loadi=[1:Busnum]';
RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数 RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数
t_Bal_volt=Volt(Balance); t_Bal_volt=Volt(Balance);
@ -17,11 +17,11 @@ 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;
wG=0.1*ones(size(PGi,1),1); wG=1*ones(size(PGi,1),1);
randInt=randperm(size(Loadi,1)); randInt=randperm(size(Loadi,1));
randPDind=randInt(1:HH); randPDind=randInt(1:10);
wD=0.1*ones(size(Loadi,1),1); wD=1*ones(size(Loadi,1),1);
wD(randPDind)=0;%一些负荷不约束 %wD(randPDind)=0;%一些负荷不约束
%wD(7)=0; %wD(7)=0;
% wD(11)=0; % wD(11)=0;
PD=1*PD0; PD=1*PD0;

4
SolveIt.m
View File

@ -25,6 +25,10 @@ aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,Co
deltG(t+Balance,:)=0; deltG(t+Balance,:)=0;
%% %%
dxdy=aa\yy; dxdy=aa\yy;
%% KLU
%spy(aa)
%dxdy = klu(aa,'\',full(yy));
%%
dX=dxdy(1:ContrlCount); dX=dxdy(1:ContrlCount);
dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
dL=Lz+deltG'*dX; dL=Lz+deltG'*dX;

8
pf.m
View File

@ -6,7 +6,7 @@ function [kmax,Precision,Uangle,U,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Li
% 程序编者： % 程序编者：
% 编制时间：2010.12 % 编制时间：2010.12
%************************************************************************** %**************************************************************************
clc; %clc;
tic; tic;
%% 读取数据文件 %% 读取数据文件
[Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,... [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
@ -15,18 +15,18 @@ tic;
[GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,... [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
Transforx,Transfork0,Branchi,Branchb); Transforx,Transfork0,Branchi,Branchb);
[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量 [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量
disp('迭代次数i 最大不平衡量'); %disp('迭代次数i 最大不平衡量');
%% 循环体计算 %% 循环体计算
for i = 0:kmax for i = 0:kmax
[Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵 [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵
% disp('第一次雅克比'); % disp('第一次雅克比');
m = max(abs(PQ)); m = max(abs(PQ));
m=full(m); m=full(m);
fprintf(' %u %.8f \n',i,m); %fprintf(' %u %.8f \n',i,m);
if m > Precision %判断不平衡量是否满足精度要求 if m > Precision %判断不平衡量是否满足精度要求
[Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量 [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量
else else
disp(['收敛，迭代次数为',num2str(i),'次']); %disp(['收敛，迭代次数为',num2str(i),'次']);
break %若满足精度要求，则计算收敛 break %若满足精度要求，则计算收敛
end end
end end