From 0f8edee30bbc03a18310c6b9e0dc780550d0e2c8 Mon Sep 17 00:00:00 2001 From: facat Date: Mon, 21 Jan 2013 21:47:54 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E5=85=A5=E4=BA=86=E5=90=84=E7=A7=8D?= =?UTF-8?q?=E7=BB=9F=E8=AE=A1=E5=80=BC=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- @ForThesis/AddPDQDPGQG.m | 14 +++ @ForThesis/ForThesis.m | 30 +++++ @ForThesis/MaxDeviation.asv | 7 ++ @ForThesis/MaxDeviation.m | 11 ++ @ForThesis/MeanPD.m | 7 ++ @ForThesis/MeanPG.m | 7 ++ @ForThesis/MeanQD.m | 7 ++ @ForThesis/MeanQG.m | 7 ++ @ForThesis/StatDeviation.asv | 15 +++ @ForThesis/StatDeviation.m | 20 +++ FormLw.m | 10 +- FormLz.m | 10 +- OPF.asv | 7 +- OPF.m | 238 ++++++++++++++++++----------------- OPF_Init.m | 4 +- 15 files changed, 266 insertions(+), 128 deletions(-) create mode 100644 @ForThesis/AddPDQDPGQG.m create mode 100644 @ForThesis/ForThesis.m create mode 100644 @ForThesis/MaxDeviation.asv create mode 100644 @ForThesis/MaxDeviation.m create mode 100644 @ForThesis/MeanPD.m create mode 100644 @ForThesis/MeanPG.m create mode 100644 @ForThesis/MeanQD.m create mode 100644 @ForThesis/MeanQG.m create mode 100644 @ForThesis/StatDeviation.asv create mode 100644 @ForThesis/StatDeviation.m diff --git a/@ForThesis/AddPDQDPGQG.m b/@ForThesis/AddPDQDPGQG.m new file mode 100644 index 0000000..f5a179c --- /dev/null +++ b/@ForThesis/AddPDQDPGQG.m @@ -0,0 +1,14 @@ +function this=AddPDQDPGQG( this, PD,QD,PG,QG) +%ADDPD Summary of this function goes here +% Detailed explanation goes here +if this.currentPos>this.sampleNum + error('出入的数太多了'); + return; +end +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; +end + diff --git a/@ForThesis/ForThesis.m b/@ForThesis/ForThesis.m new file mode 100644 index 0000000..8cb1866 --- /dev/null +++ b/@ForThesis/ForThesis.m @@ -0,0 +1,30 @@ +classdef ForThesis + %UNTITLED Summary of this class goes here + % Detailed explanation goes here + + properties + currentPos=0; + sampleNum=0; + LoadNum=0; + PDArray=0; + + QDArray=0; + PGArray=0; + QGArray=0; + + end + + methods + function this=ForThesis(sampleNum,LoadNum) + this.currentPos=0;%记录当前用到的列，不要超过100列 + this.sampleNum=sampleNum; + this.LoadNum=LoadNum; + this.PDArray=zeros(LoadNum,sampleNum); + this.QDArray=zeros(LoadNum,sampleNum); + this.PGArray=zeros(sampleNum,1); + this.QGArray=zeros(sampleNum,1); + end + + end + +end \ No newline at end of file diff --git a/@ForThesis/MaxDeviation.asv b/@ForThesis/MaxDeviation.asv new file mode 100644 index 0000000..bcad661 --- /dev/null +++ b/@ForThesis/MaxDeviation.asv @@ -0,0 +1,7 @@ +function MaxDeviation(this,PG0,QG0,PD0,QD0) +PD0Array=repmat(PD0,1,this.sampleNum); +QD0Array=repmat(QD0,1,this.sampleNum); +PDMaxDev=max(abs(this.PDArray-PD0Array),[],2) +QDMaxDev=max(abs(this.QDArray-QD0Array),[],2) +PG0Array=repmat() +end diff --git a/@ForThesis/MaxDeviation.m b/@ForThesis/MaxDeviation.m new file mode 100644 index 0000000..006d74e --- /dev/null +++ b/@ForThesis/MaxDeviation.m @@ -0,0 +1,11 @@ +function [output_arg]=MaxDeviation(this,PG0,QG0,PD0,QD0) +PD0Array=repmat(PD0,1,this.sampleNum); +QD0Array=repmat(QD0,1,this.sampleNum); +PDMaxDev=max(abs( (this.PDArray-PD0Array)./PD0Array ),[],2); +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]); +end diff --git a/@ForThesis/MeanPD.m b/@ForThesis/MeanPD.m new file mode 100644 index 0000000..fca29ca --- /dev/null +++ b/@ForThesis/MeanPD.m @@ -0,0 +1,7 @@ +function [ output_args ] = MeanPD(this) +%MEANPD Summary of this function goes here +% Detailed explanation goes here +output_args=sum(this.PDArray,2)/this.sampleNum; + +end + diff --git a/@ForThesis/MeanPG.m b/@ForThesis/MeanPG.m new file mode 100644 index 0000000..a6e12bd --- /dev/null +++ b/@ForThesis/MeanPG.m @@ -0,0 +1,7 @@ +function [ output_args ] = MeanPG( this ) +%MEANPG Summary of this function goes here +% Detailed explanation goes here +output_args=sum(this.PGArray,1)/this.sampleNum; + +end + diff --git a/@ForThesis/MeanQD.m b/@ForThesis/MeanQD.m new file mode 100644 index 0000000..e976e9c --- /dev/null +++ b/@ForThesis/MeanQD.m @@ -0,0 +1,7 @@ +function [ output_args ] = MeanQD(this) +%MEANQD Summary of this function goes here +% Detailed explanation goes here +output_args=sum(this.QDArray,2)/this.sampleNum; + +end + diff --git a/@ForThesis/MeanQG.m b/@ForThesis/MeanQG.m new file mode 100644 index 0000000..97aa772 --- /dev/null +++ b/@ForThesis/MeanQG.m @@ -0,0 +1,7 @@ +function [ output_args ] = MeanQG( this ) +%MEAQG Summary of this function goes here +% Detailed explanation goes here +output_args=sum(this.QGArray,1)/this.sampleNum; + +end + diff --git a/@ForThesis/StatDeviation.asv b/@ForThesis/StatDeviation.asv new file mode 100644 index 0000000..7687bc2 --- /dev/null +++ b/@ForThesis/StatDeviation.asv @@ -0,0 +1,15 @@ +function [ output_args ] = StatDeviation( this,PG0,QG0,PD0,QD0 )%统计误差 +%STATDEVIATION Summary of this function goes here +% Detailed explanation goes here + +PD0Array=repmat(PD0,1,this.sampleNum); +QD0Array=repmat(QD0,1,this.sampleNum); +PDMaxDev=max(abs(this.PDArray-PD0Array),[],2); +QDMaxDev=max(abs(this.QDArray-QD0Array),[],2); +PG0Array=repmat(PG0,this.sampleNum,1); +QG0Array=repmat(QG0,this.sampleNum,1); +PGMaxDev=max(abs(PG0Array)); +QGMaxDev=max(abs(QG0Array)); + +end + diff --git a/@ForThesis/StatDeviation.m b/@ForThesis/StatDeviation.m new file mode 100644 index 0000000..abd08ed --- /dev/null +++ b/@ForThesis/StatDeviation.m @@ -0,0 +1,20 @@ +function [ output_args ] = StatDeviation( this,PG0,QG0,PD0,QD0 )%统计误差 +%STATDEVIATION Summary of this function goes here +% Detailed explanation goes here + +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; +QGDev=(QG0Array-this.QGArray)/0.01; +wholeMat=[PDDev;QDDev;PGDev';QGDev']; +t1=wholeMat.^2; +t2=sum(t1,1); +t3=t2/size(t1,1); +t4=t3.^2; +output_args=sum(t4)/length(t4); +end + diff --git a/FormLw.m b/FormLw.m index d4811a9..056fd3c 100644 --- a/FormLw.m +++ b/FormLw.m @@ -8,14 +8,14 @@ VoltU=(1.1)*ones(1,Busnum); %VoltU=10*ones(1,Busnum); PDU=PD0(Loadi); % PDU=noDataTransCapacity; -PDU(PDU>0)=1.200*PDU(PDU>0); -PDU(PDU<0)=0.800*PDU(PDU<0); +PDU(PDU>0)=1.300*PDU(PDU>0); +PDU(PDU<0)=0.700*PDU(PDU<0); PDU(PDU==0)=0.400; %PDU=10*ones(length(Loadi),1); QDU=QD0(Loadi); -QDU(QDU>0)=1.200*QDU(QDU>0); -QDU(QDU<0)=0.800*QDU(QDU<0); -QDU(QDU==0)=0.200; +QDU(QDU>0)=1.300*QDU(QDU>0); +QDU(QDU<0)=0.700*QDU(QDU<0); +QDU(QDU==0)=0.400; % PF=0.85; % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF; t1=([PU',QU',PDU',QDU',VoltU])'; diff --git a/FormLz.m b/FormLz.m index 4229485..8bdd741 100644 --- a/FormLz.m +++ b/FormLz.m @@ -7,14 +7,14 @@ QL=-5*ones(length(PVQL(:,1)),1); VoltL=(0.9)*ones(1,Busnum); %VoltL=-10*ones(1,Busnum); PDL=PD0(Loadi); -PDL(PDL>0)=0.800*PDL(PDL>0); -PDL(PDL<0)=1.200*PDL(PDL<0); +PDL(PDL>0)=0.700*PDL(PDL>0); +PDL(PDL<0)=1.300*PDL(PDL<0); PDL(PDL==0)=-0.400; %PDL=-10*ones(length(Loadi),1); QDL=QD0(Loadi); -QDL(QDL>0)=0.800*QDL(QDL>0); -QDL(QDL<0)=1.200*QDL(QDL<0); -QDL(QDL==0)=-0.200; +QDL(QDL>0)=0.700*QDL(QDL>0); +QDL(QDL<0)=1.300*QDL(QDL<0); +QDL(QDL==0)=-0.400; % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF; t1=([PL',QL',PDL',QDL',VoltL])'; t2=Mat_G-Init_L'-t1; diff --git a/OPF.asv b/OPF.asv index 90ab4c9..551113e 100644 --- a/OPF.asv +++ b/OPF.asv @@ -51,7 +51,8 @@ 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(fileName); +%[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; @@ -83,7 +84,7 @@ while(abs(Gap)>Precision) 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); + 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); %% 开始解方程 @@ -97,7 +98,7 @@ while(abs(Gap)>Precision) 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); +% DrawGap(plotGap); %% %Volt=full(Volt'); %PD=full(PD); diff --git a/OPF.m b/OPF.m index 79e8e75..0fec178 100644 --- a/OPF.m +++ b/OPF.m @@ -1,120 +1,132 @@ 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\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视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'; -%% 初值-即测量值 -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\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-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; + 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=KK+1; + 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/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,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); + 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)-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); + 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(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi)); +thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi)) +toc diff --git a/OPF_Init.m b/OPF_Init.m index cc1284f..8b555d5 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -17,8 +17,8 @@ tPL=sparse(GenL(:,2));% tQL=sparse(PVQL(:,1));% 无功下限 PG(PGi)=(tPU+tPL)/2; QG(PVi)=(tQU+tQL)/2; -wPG=100*ones(size(PGi,1),1); -wQG=100*ones(size(PVi,1),1); +wPG=0*ones(size(PGi,1),1); +wQG=0*ones(size(PVi,1),1); %randInt=randperm(size(Loadi,1)); %randPDind=randInt(1:10); randPDind=0;