From 798e1e88ac37ccac6fadf29c581547c68a7e85ea Mon Sep 17 00:00:00 2001 From: unknown Date: Wed, 25 Jul 2012 15:08:54 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E9=94=99=E8=AF=AF=E6=95=B0=E6=8D=AE?= =?UTF-8?q?=E8=AF=86=E5=88=AB?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- FormH.m | 2 +- Modification.m | 4 +- OPF.m | 18 +++++-- OPFWrongDataCheck.asv | 106 ++++++++++++++++++++++++++++++++++++++++++ OPFWrongDataCheck.m | 106 ++++++++++++++++++++++++++++++++++++++++++ OPF_Init.m | 12 ++--- SolveIt.m | 4 ++ pf.m | 8 ++-- 8 files changed, 242 insertions(+), 18 deletions(-) create mode 100644 OPFWrongDataCheck.asv create mode 100644 OPFWrongDataCheck.m diff --git a/FormH.m b/FormH.m index 6d5cb72..88fee3d 100644 --- a/FormH.m +++ b/FormH.m @@ -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) %% 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; %% AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); diff --git a/Modification.m b/Modification.m index 65f89ef..2e4b6f0 100644 --- a/Modification.m +++ b/Modification.m @@ -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) 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); -fprintf('AlphaD %f\n',full(AlphaD)); +%fprintf('AlphaD %f\n',full(AlphaD)); Init_Z=Init_Z+AlphaD*deltZ'; Init_L=Init_L+AlphaP*deltL'; diff --git a/OPF.m b/OPF.m index c5802c7..c67ef47 100644 --- a/OPF.m +++ b/OPF.m @@ -1,9 +1,9 @@ tic 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]= ... -pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); +pf('D:\Project\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt'); %pf('c:/file31.txt'); -%pf('ieee10471PG.dat'); +%pf('ieee118PG.dat'); %% 计算功率因数 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); PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt'; -%% 初值 +%% 初值-即测量值 PG0=PG; PD0=PD; + +PDReal=PD;%真值 +PD0(12)=PD0(12)+0.001; %% 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'); KK=0; plotGap=zeros(1,50); @@ -78,6 +81,11 @@ DrawGap(plotGap); %Volt=full(Volt'); %PD=full(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 diff --git a/OPFWrongDataCheck.asv b/OPFWrongDataCheck.asv new file mode 100644 index 0000000..04bce85 --- /dev/null +++ b/OPFWrongDataCheck.asv @@ -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 + diff --git a/OPFWrongDataCheck.m b/OPFWrongDataCheck.m new file mode 100644 index 0000000..16b6b1a --- /dev/null +++ b/OPFWrongDataCheck.m @@ -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 + diff --git a/OPF_Init.m b/OPF_Init.m index 0ca0a4f..aca5a39 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -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) -Loadi=find(QD~=0); +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 & PD~=0); %Loadi=[1:Busnum]'; RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数 t_Bal_volt=Volt(Balance); @@ -17,11 +17,11 @@ tPL=sparse(GenL(:,2));% tQL=sparse(PVQL(:,1));% 无功下限 PG(PGi)=(tPU+tPL)/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)); -randPDind=randInt(1:HH); -wD=0.1*ones(size(Loadi,1),1); -wD(randPDind)=0;%一些负荷不约束 +randPDind=randInt(1:10); +wD=1*ones(size(Loadi,1),1); +%wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0; PD=1*PD0; diff --git a/SolveIt.m b/SolveIt.m index 23d81cc..4d0662c 100644 --- a/SolveIt.m +++ b/SolveIt.m @@ -25,6 +25,10 @@ aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,Co deltG(t+Balance,:)=0; %% dxdy=aa\yy; +%% KLU +%spy(aa) +%dxdy = klu(aa,'\',full(yy)); +%% dX=dxdy(1:ContrlCount); dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum); dL=Lz+deltG'*dX; diff --git a/pf.m b/pf.m index 575e701..9cd1f77 100644 --- a/pf.m +++ b/pf.m @@ -6,7 +6,7 @@ function [kmax,Precision,Uangle,U,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Li % 程序编者： % 编制时间：2010.12 %************************************************************************** -clc; +%clc; tic; %% 读取数据文件 [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,... Transforx,Transfork0,Branchi,Branchb); [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量 -disp('迭代次数i 最大不平衡量'); +%disp('迭代次数i 最大不平衡量'); %% 循环体计算 for i = 0:kmax [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c); %形成雅克比矩阵 % disp('第一次雅克比'); m = max(abs(PQ)); m=full(m); - fprintf(' %u %.8f \n',i,m); + %fprintf(' %u %.8f \n',i,m); if m > Precision %判断不平衡量是否满足精度要求 [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量 else - disp(['收敛，迭代次数为',num2str(i),'次']); + %disp(['收敛，迭代次数为',num2str(i),'次']); break %若满足精度要求，则计算收敛 end end