diff --git a/0PD0.mat b/0PD0.mat new file mode 100644 index 0000000..c78db61 Binary files /dev/null and b/0PD0.mat differ diff --git a/0QD0.mat b/0QD0.mat new file mode 100644 index 0000000..19bc650 Binary files /dev/null and b/0QD0.mat differ diff --git a/12PD0.mat b/12PD0.mat new file mode 100644 index 0000000..85525f4 Binary files /dev/null and b/12PD0.mat differ diff --git a/12QD0.mat b/12QD0.mat new file mode 100644 index 0000000..fe326ed Binary files /dev/null and b/12QD0.mat differ diff --git a/16PD0.mat b/16PD0.mat new file mode 100644 index 0000000..37a7220 Binary files /dev/null and b/16PD0.mat differ diff --git a/16QD0.mat b/16QD0.mat new file mode 100644 index 0000000..d03a88e Binary files /dev/null and b/16QD0.mat differ diff --git a/20PD0.mat b/20PD0.mat new file mode 100644 index 0000000..066fff1 Binary files /dev/null and b/20PD0.mat differ diff --git a/20QD0.mat b/20QD0.mat new file mode 100644 index 0000000..1d7dcd4 Binary files /dev/null and b/20QD0.mat differ diff --git a/4PD0.mat b/4PD0.mat new file mode 100644 index 0000000..7181fd5 Binary files /dev/null and b/4PD0.mat differ diff --git a/4QD0.mat b/4QD0.mat new file mode 100644 index 0000000..9edccdc Binary files /dev/null and b/4QD0.mat differ diff --git a/8PD0.mat b/8PD0.mat new file mode 100644 index 0000000..486594f Binary files /dev/null and b/8PD0.mat differ diff --git a/8QD0.mat b/8QD0.mat new file mode 100644 index 0000000..f8fc076 Binary files /dev/null and b/8QD0.mat differ diff --git a/@ForThesis/MaxBranchDeviation.m b/@ForThesis/MaxBranchDeviation.m index dfdfca4..fadccac 100644 --- a/@ForThesis/MaxBranchDeviation.m +++ b/@ForThesis/MaxBranchDeviation.m @@ -1,7 +1,34 @@ -function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel) +function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,FileName,PD0,QD0) %% 最大支路功率偏差 % 支路功率包括线路和变压器 -%% +%% 先用加了误差的负荷功率计算潮流值 +[Busnum,Balance,PQstandard,Precision,~,~,~,~,~,kmax,~ ,... + ~,~,~,~,~,~,~,Pointpoweri,PG,QG,~,~,PVi,PVu,~,~,~,~,~,~,~]= openfile2(FileName); +PD=PD0; +QD=QD0; +%% 形成节点导纳矩阵 +[~,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,... + Transforx,Transfork0,Branchi,Branchg,Branchb); +[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum); %求功率不平衡量 +%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); + if m > Precision %判断不平衡量是否满足精度要求 + [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U); %求解修正方程，更新电压变量 + else + disp(['收敛，迭代次数为',num2str(i),'次']); + break %若满足精度要求，则计算收敛 + end +end + +Volt0=U; +UAngel0=Uangle; + [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); diff --git a/@ForThesis/MaxDeviation.m b/@ForThesis/MaxDeviation.m index 006d74e..8bfa6c3 100644 --- a/@ForThesis/MaxDeviation.m +++ b/@ForThesis/MaxDeviation.m @@ -7,5 +7,5 @@ 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]); +output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev;QGMaxDev;]); end diff --git a/@ForThesis/StatBranchDeviation.m b/@ForThesis/StatBranchDeviation.m index 4acc382..6e9ff72 100644 --- a/@ForThesis/StatBranchDeviation.m +++ b/@ForThesis/StatBranchDeviation.m @@ -14,6 +14,6 @@ t22=dispTransloss0(:,3)>1e-5;% 太小 t3=[t1(t11);t2(t22)]; t4=t3.^2; t5=sum(t4).^.5; -output_args=t5 +output_args=t5; end diff --git a/FormLw.m b/FormLw.m index 5088230..8c02646 100644 --- a/FormLw.m +++ b/FormLw.m @@ -6,13 +6,13 @@ VoltU=(1.1)*ones(1,Busnum); %VoltU=10*ones(1,Busnum); PDU=PD0(Loadi); % PDU=noDataTransCapacity; -PDU(PDU>0)=1.300*PDU(PDU>0); -PDU(PDU<0)=0.700*PDU(PDU<0); +PDU(PDU>0)=1.200*PDU(PDU>0); +PDU(PDU<0)=0.800*PDU(PDU<0); PDU(PDU==0)=0.400; %PDU=10*ones(length(Loadi),1); QDU=QD0(Loadi); -QDU(QDU>0)=1.300*QDU(QDU>0); -QDU(QDU<0)=0.700*QDU(QDU<0); +QDU(QDU>0)=1.200*QDU(QDU>0); +QDU(QDU<0)=0.800*QDU(QDU<0); QDU(QDU==0)=0.400; % PF=0.85; % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF; diff --git a/FormLz.m b/FormLz.m index 48a736e..1f58ebc 100644 --- a/FormLz.m +++ b/FormLz.m @@ -4,13 +4,13 @@ KK=999; VoltL=(0.9)*ones(1,Busnum); %VoltL=-10*ones(1,Busnum); PDL=PD0(Loadi); -PDL(PDL>0)=0.700*PDL(PDL>0); -PDL(PDL<0)=1.300*PDL(PDL<0); +PDL(PDL>0)=0.800*PDL(PDL>0); +PDL(PDL<0)=1.200*PDL(PDL<0); PDL(PDL==0)=-0.400; %PDL=-10*ones(length(Loadi),1); QDL=QD0(Loadi); -QDL(QDL>0)=0.700*QDL(QDL>0); -QDL(QDL<0)=1.300*QDL(QDL<0); +QDL(QDL>0)=0.800*QDL(QDL>0); +QDL(QDL<0)=1.200*QDL(QDL<0); QDL(QDL==0)=-0.400; % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF; t1=([PDL',QDL',VoltL])'; diff --git a/Lineloss.asv b/Lineloss.asv index 3132a3d..16e5eaa 100644 --- a/Lineloss.asv +++ b/Lineloss.asv @@ -1,10 +1,10 @@ %% 计算线损 -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有功') -cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle); -cmpVolt=cmpVolt'; +% fprintf('功率为有名值\n'); +% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar') +cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle'); +% cmpVolt=cmpVolt'; y0=1i*Lineb2; %yj0=1i*standardinput(:,7); yij=1./(Liner+1i*Linex); @@ -13,16 +13,26 @@ Sij=cmpVolt(Linei).*conj( cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Line Sji=cmpVolt(Linej).*conj( cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij ); %Sij1==Sij2 deltLineS=Sij+Sji; + +%% 另一种计算方式begin +ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij)); +ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij); +ss=Volt(Linei)'.*Volt(Linej)'.*abs(yij).*cos() +%% 另一种计算方式end dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100]; %full(dispLineloss) -dispLineloss=sortrows(dispLineloss,-3); -full(dispLineloss) +% dispLineloss=sortrows(dispLineloss,-3); +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) +% dispTransloss=sortrows(dispTransloss,-3); +full(dispTransloss); end \ No newline at end of file diff --git a/Lineloss.m b/Lineloss.m index 3af8837..f2c1428 100644 --- a/Lineloss.m +++ b/Lineloss.m @@ -3,8 +3,8 @@ function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Tr %format long % fprintf('功率为有名值\n'); % fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar') -cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle); -cmpVolt=cmpVolt'; +cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle'); +% cmpVolt=cmpVolt'; y0=1i*Lineb2; %yj0=1i*standardinput(:,7); yij=1./(Liner+1i*Linex); @@ -13,6 +13,12 @@ Sij=cmpVolt(Linei).*conj( cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Line Sji=cmpVolt(Linej).*conj( cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij ); %Sij1==Sij2 deltLineS=Sij+Sji; + +%% 另一种计算方式begin +% ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij)); +% ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij); +% ss=abs(yij).*cos(angle(yij)).*( Volt(Linei)'.^2+Volt(Linej)'.^2- 2*Volt(Linei)'.*Volt(Linej)' .*cos(Angle(Linei)' - Angle(Linej)') ); +%% 另一种计算方式end dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100]; %full(dispLineloss) dispLineloss=sortrows(dispLineloss,-3); @@ -24,7 +30,8 @@ Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij ); deltTransS=Sij+Sji; %% 接地支路损耗 % 没有考虑变压器变比 -deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb)))); +deltTransS =deltTransS+cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))); +ss=Volt(Branchi)'.^2.*(Branchg) %% dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100]; dispTransloss=sortrows(dispTransloss,-3); diff --git a/OPF.m b/OPF.m index 76f63d0..f41efec 100644 --- a/OPF.m +++ b/OPF.m @@ -4,15 +4,13 @@ clear %% 存在问题 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123 %% -thesis=ForThesis(4,8); -for II=1:4 +thesis=ForThesis(1,45); +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('E:\算例\东际911_2751267_2012-09-05\newFIle20.txt'); + pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle20.txt'); %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); - %pf('c:/file31.txt'); - - + %pf('c:/file31.txt'); %% 计算功率因数 Loadi=QD~=0 | PD~=0; PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2)) @@ -53,17 +51,16 @@ for II=1:4 QD_NON_ZERO=QD(PD==0 & QD~=0); QD_NON_ZERO_IND=find(PD==0 & QD~=0); %% - Precision=Precision/10; - + 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)); + PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.5,length(Loadi),1)); + QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.5,length(Loadi),1)); + load('0PD0.mat'); + load('0QD0.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)); %% 读变压器容量 - %[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; @@ -87,7 +84,6 @@ for II=1:4 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); @@ -136,12 +132,12 @@ 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)) +thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi)) +thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(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.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle0.txt',PD0,QD0) thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel) toc diff --git a/OPF_Init.asv b/OPF_Init.asv index 61a11ab..c0cfc2d 100644 --- a/OPF_Init.asv +++ b/OPF_Init.asv @@ -1,7 +1,7 @@ -function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD) +function [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) Loadi=find(QD~=0 | PD~=0); %Loadi=[1:Busnum]'; -RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1+1; %约束条件数,只放开一个QD +RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD t_Bal_volt=Volt(Balance); Volt=sparse(1*ones(1,Busnum)); Volt(Balance)=t_Bal_volt; @@ -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; -QG(PVi)=(tQU+tQL)/2; -wPG=1*ones(size(PGi,1),1); -wQG=1*zeros(size(PVi,1),1); +wPG=0; +wQG=0; %randInt=randperm(size(Loadi,1)); %randPDind=randInt(1:10); randPDind=0; -wD=1*zeros(size(Loadi,1),1); +wPD=1/.5^2*zeros(size(Loadi,1),1); +wPD(1:2:end)=0; +wQD=1/.5^2*ones(size(Loadi,1),1); +wQD(1:2:end)=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 \ No newline at end of file diff --git a/OPF_Init.m b/OPF_Init.m index b9fdfc6..6dfd73d 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -20,10 +20,10 @@ wQG=0; %randInt=randperm(size(Loadi,1)); %randPDind=randInt(1:10); randPDind=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; +wPD=1/.5^2*ones(size(Loadi,1),1); +% wPD(1:2:end)=0; +wQD=1/.5^2*ones(size(Loadi,1),1); +% wQD(1:2:end)=0; %wD(randPDind)=0;%一些负荷不约束 %wD(7)=0; % wD(11)=0;