From acab2b5d17dba936188d59d8aaf06f8c93974f47 Mon Sep 17 00:00:00 2001 From: "facat@lab.com" Date: Tue, 4 Nov 2014 17:21:06 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8A=A0=E6=B5=8B=E8=AF=95=E7=A8=8B=E5=BA=8F?= =?UTF-8?q?=E4=BB=A3=E7=A0=81=EF=BC=8C=E5=B9=B6=E5=8A=A0=E5=85=A5=E8=AF=84?= =?UTF-8?q?=E4=BB=B7=E6=8C=87=E6=A0=87=E4=BB=A3=E7=A0=81=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat@lab.com --- FormLw.m | 4 +-- FormLz.m | 4 +-- OPF.m | 92 +++++++++++++++++++++++++++++++++++--------------------- test.m | 22 ++++++++++++++ 4 files changed, 84 insertions(+), 38 deletions(-) create mode 100644 test.m diff --git a/FormLw.m b/FormLw.m index a75b7d0..c411549 100644 --- a/FormLw.m +++ b/FormLw.m @@ -26,8 +26,8 @@ function Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P) % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF; % PDU=1.8*ones(length(Loadi),1); % QDU=1.8*ones(length(Loadi),1); -PDU=1.05*rPD3P; -QDU=1.05*rQD3P; +PDU=1.1*rPD3P; +QDU=1.1*rQD3P; t1=([PDU',QDU'])'; t2=Mat_G+Init_U'-t1; Lw=t2; diff --git a/FormLz.m b/FormLz.m index 1796386..5627360 100644 --- a/FormLz.m +++ b/FormLz.m @@ -23,8 +23,8 @@ function Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P) % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF; % PDL=0*ones(length(Loadi),1); % QDL=0*ones(length(Loadi),1); -PDL=0.95*rPD3P; -QDL=0.95*rQD3P; +PDL=0.9*rPD3P; +QDL=0.9*rQD3P; t1=([PDL',QDL'])'; t2=Mat_G-Init_L'-t1; Lz=t2; diff --git a/OPF.m b/OPF.m index 21dda4a..47741df 100644 --- a/OPF.m +++ b/OPF.m @@ -1,16 +1,18 @@ + +function [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF() tic clc clear -lineZ=readLineZ('feeder123\lineParameter.txt'); +lineZ=readLineZ('feeder13\lineParameter.txt'); [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ... phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ... - cap]=dataRead(lineZ,'feeder123\data.txt'); -phaseASpotLoadP(phaseASpotLoadP==0)=0.002; -phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002; -phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002; -phaseASpotLoadQ(phaseASpotLoadQ==0)=0.002; -phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002; -phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002; + cap]=dataRead(lineZ,'feeder13\data1.txt'); +% phaseASpotLoadP(phaseASpotLoadP==0)=0.002; +% phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002; +% phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002; +% phaseASpotLoadQ(phaseASpotLoadQ==0)=0.002; +% phaseBSpotLoadQ(phaseBSpotLoadQ==0)=0.002; +% phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002; @@ -73,9 +75,9 @@ while(k<=kmax && maxD> EPS) SB=VoltpB.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1)); SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1)); %先不要电容 -% SA=0; -% SB=0; -% SC=0; + % SA=0; + % SB=0; + % SC=0; iterPD=PD+real(SA+SB+SC)/3; iterQD=QD+imag(SA+SB+SC)/3; iterPhaseASpotLoadP=phaseASpotLoadP+real(SA); @@ -151,12 +153,12 @@ PQi3P(3:3:end)=(PQi-1)*3+3; Loadi=PQi3P; PD3P=sparse(Busnum*3,1); QD3P=sparse(Busnum*3,1); -PD3P(1:3:end)=phaseASpotLoadP*0.5; -PD3P(2:3:end)=phaseBSpotLoadP*0.5; -PD3P(3:3:end)=phaseCSpotLoadP*0.5; -QD3P(1:3:end)=phaseASpotLoadQ*0.5; -QD3P(2:3:end)=phaseBSpotLoadQ*0.5; -QD3P(3:3:end)=phaseCSpotLoadQ*0.5; +PD3P(1:3:end)=phaseASpotLoadP*0.9; +PD3P(2:3:end)=phaseBSpotLoadP*0.9; +PD3P(3:3:end)=phaseCSpotLoadP*0.9; +QD3P(1:3:end)=phaseASpotLoadQ*0.9; +QD3P(2:3:end)=phaseBSpotLoadQ*0.9; +QD3P(3:3:end)=phaseCSpotLoadQ*0.9; PD3P=PD3P(Loadi); QD3P=QD3P(Loadi); QGi3P=zeros(length(QGi)*3,1); @@ -192,10 +194,10 @@ rVoltABCV(1:3:end)=abs(VoltpABC(1,:)); rVoltABCV(2:3:end)=abs(VoltpABC(2,:)); rVoltABCV(3:3:end)=abs(VoltpABC(3,:)); %三相电压相角 -rVoltABCA=zeros(busNum*3,1); -rVoltABCA(1:3:end)=angle(VoltpABC(1,:)); -rVoltABCA(2:3:end)=angle(VoltpABC(2,:)); -rVoltABCA(3:3:end)=angle(VoltpABC(3,:)); +% rVoltABCA=zeros(busNum*3,1); +% rVoltABCA(1:3:end)=angle(VoltpABC(1,:)); +% rVoltABCA(2:3:end)=angle(VoltpABC(2,:)); +% rVoltABCA(3:3:end)=angle(VoltpABC(3,:)); %三相负荷 rPD3P=zeros(busNum*3,1); rPD3P(1:3:end)=phaseASpotLoadP; @@ -216,9 +218,12 @@ mPD3P=rPD3P.*PD3PSigma; QD3PSigma=(1+normrnd(0,sigma,length(rQD3P),1)); mQD3P=rPD3P.*QD3PSigma; %量测方差 -wVolt=1./(abs(mVoltABCV*sigma).^2)*1; -wPD=1./(abs(mPD3P*sigma).^2)*1; -wQD=1./(abs(mQD3P*sigma).^2)*1; +wVolt=1./(abs(mVoltABCV*sigma).^2); +wPD=1./(abs(mPD3P*sigma).^2); +wQD=1./(abs(mQD3P*sigma).^2); +wVolt(1:3:end)=0; +wPD(1:3:end)=0; +wQD(1:3:end)=0; % wVolt=1./abs(VoltSigma).^2; % wPD=1./abs(PD3PSigma).^2; % wQD=1./abs(QD3PSigma).^2; @@ -244,26 +249,26 @@ while(abs(Gap)>Precision) break; end Init_u=Gap/2/RestraintCount*CenterA; - AngleIJMat=0; +% AngleIJMat=0; %% 开始计算OPF %% 形成等式约束的雅克比 deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi); %% 形成不等式约束的雅克比 deltG=func_deltG(Busnum,Loadi); %% - L_1Z=diag(Init_Z./Init_L); - U_1W=diag(Init_W./Init_U); +% L_1Z=diag(Init_Z./Init_L); +% U_1W=diag(Init_W./Init_U); %% 形成海森阵 deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount); -% deltdeltF=0; +% deltdeltF=0; %% 形成ddHy ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Yangle,Loadi,ContrlCount); %% 开始构建ddg -% ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD); + % ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD); ddg=0; %% 开始构建deltF deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,Busnum,Loadi); -% deltF=0; +% deltF=0; %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); @@ -273,9 +278,9 @@ while(abs(Gap)>Precision) Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P); Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); -% YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); + % YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 -% fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); + % 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,Lx,Balance,Busnum,Loadi); %%取各分量 [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); @@ -284,7 +289,26 @@ while(abs(Gap)>Precision) fprintf('%f\n',full(Gap)); KK=KK+1; end -(rVoltABCV-Volt)./rVoltABCV*100 +(rVoltABCV-Volt)./rVoltABCV*100; +%% 计算统计量 +%目标函数均值 +JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2); +%估计误差统计 +%量测量数量 +mCount=Busnum*3+length(Loadi)*3*2; +%估计量质量 +%有PD，QD的 +SEMeasurement=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((rPD3P-PD3P)./mPD3P./sigma).^2)+sum(((rQD3P-QD3P)./mQD3P./sigma).^2); +SEMeasurement=(SEMeasurement/mCount)^.5; +%没有PD，QD的 +SEMeasurementNOPQ=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2); +SEMeasurementNOPQ=(SEMeasurementNOPQ/mCount)^.5; +%量测量质量 +SEMeasurementM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2)+sum(((rPD3P-mPD3P)./mPD3P./sigma).^2)+sum(((rQD3P-mQD3P)./mQD3P./sigma).^2); +SEMeasurementM=(SEMeasurementM/mCount)^.5; +%没有PD，QD的 +SEMeasurementNOPQM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2); +SEMeasurementNOPQM=(SEMeasurementNOPQM/mCount)^.5; fprintf('迭代次数%d\n',KK); toc - +end diff --git a/test.m b/test.m new file mode 100644 index 0000000..380621f --- /dev/null +++ b/test.m @@ -0,0 +1,22 @@ +close all +clear +clc +JMeasurementSum=0; +SEMeasurementSum=0; +SEMeasurementNOPQSum=0; +SEMeasurementMSum=0; +SEMeasurementNOPQMSum=0; +N=1000; +for I=1:N + [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF(); + JMeasurementSum=JMeasurement+JMeasurementSum; + SEMeasurementSum=SEMeasurement+SEMeasurementSum; + SEMeasurementNOPQSum=SEMeasurementNOPQ+SEMeasurementNOPQSum; + SEMeasurementMSum=SEMeasurementMSum+SEMeasurementM; + SEMeasurementNOPQMSum=SEMeasurementNOPQMSum+SEMeasurementNOPQM; +end +JMeasurementSum=JMeasurementSum/N; +SEMeasurementSum=SEMeasurementSum/N; +SEMeasurementNOPQSum=SEMeasurementNOPQSum/N; +SEMeasurementMSum=SEMeasurementMSum/N; +SEMeasurementNOPQMSum=SEMeasurementNOPQMSum/N; \ No newline at end of file