facat
/
ipmthreephasese
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

1.加了统计功能 

2.用了新算例
This commit is contained in:
dmy@lab 2015-04-02 22:42:17 +08:00
parent f255d7285d
commit 7652f913ea
7 changed files with 114 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CalPGQG.m
View File

@ -1,5 +1,5 @@
function [ output_args ] = CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ ) function [ output_args ] = CalPGQG(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ )
%%利用Fortiscue方法计算得到了相电压和相角后，反推负荷功率，检查是否和给定的功率一致。 %计算PG和QG
busNum=size(phaseABCY,1)/3; busNum=size(phaseABCY,1)/3;
VoltpA=VoltpABC(1,:); VoltpA=VoltpABC(1,:);
VoltpB=VoltpABC(2,:); VoltpB=VoltpABC(2,:);

53
OPF.m
View File

@ -1,5 +1,5 @@
function [JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF() function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=OPF()
tic tic
clc clc
clear clear
@ -137,7 +137,7 @@ fprintf('
busNum=length(nodeNum); busNum=length(nodeNum);
Busnum=busNum; Busnum=busNum;
PQi=setxor(nodeNum,Balance); % PQi=setxor(nodeNum,Balance);
QGi=[Balance]; QGi=[Balance];
% fprintf('开始牛顿法迭代\n'); % fprintf('开始牛顿法迭代\n');
[r,c,GB]=find(phaseABCY); [r,c,GB]=find(phaseABCY);
@ -146,10 +146,10 @@ Yangle=angle(GB);
Vp3=sparse(ones(busNum*3,1));%给电压赋初值 Vp3=sparse(ones(busNum*3,1));%给电压赋初值
Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi); Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi);
Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi); Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi);
PQi3P=zeros(length(PQi)*3,1); PQi3P=zeros(length(Loadi)*3,1);
PQi3P(1:3:end)=(PQi-1)*3+1; PQi3P(1:3:end)=(Loadi-1)*3+1;
PQi3P(2:3:end)=(PQi-1)*3+2; PQi3P(2:3:end)=(Loadi-1)*3+2;
PQi3P(3:3:end)=(PQi-1)*3+3; PQi3P(3:3:end)=(Loadi-1)*3+3;
Loadi=PQi3P; Loadi=PQi3P;
PD3P=sparse(Busnum*3,1); PD3P=sparse(Busnum*3,1);
QD3P=sparse(Busnum*3,1); QD3P=sparse(Busnum*3,1);
@ -194,10 +194,10 @@ rVoltABCV(1:3:end)=abs(VoltpABC(1,:));
rVoltABCV(2:3:end)=abs(VoltpABC(2,:)); rVoltABCV(2:3:end)=abs(VoltpABC(2,:));
rVoltABCV(3:3:end)=abs(VoltpABC(3,:)); rVoltABCV(3:3:end)=abs(VoltpABC(3,:));
%三相电压相角 %三相电压相角
% rVoltABCA=zeros(busNum*3,1); rVoltABCA=zeros(busNum*3,1);
% rVoltABCA(1:3:end)=angle(VoltpABC(1,:)); rVoltABCA(1:3:end)=angle(VoltpABC(1,:));
% rVoltABCA(2:3:end)=angle(VoltpABC(2,:)); rVoltABCA(2:3:end)=angle(VoltpABC(2,:));
% rVoltABCA(3:3:end)=angle(VoltpABC(3,:)); rVoltABCA(3:3:end)=angle(VoltpABC(3,:));
%三相负荷 %三相负荷
rPD3P=zeros(busNum*3,1); rPD3P=zeros(busNum*3,1);
rPD3P(1:3:end)=phaseASpotLoadP; rPD3P(1:3:end)=phaseASpotLoadP;
@ -221,11 +221,11 @@ mQD3P=rPD3P.*QD3PSigma;
wVolt=1./(abs(mVoltABCV*sigma).^2); wVolt=1./(abs(mVoltABCV*sigma).^2);
wPD=1./(abs(mPD3P*sigma).^2); wPD=1./(abs(mPD3P*sigma).^2);
wQD=1./(abs(mQD3P*sigma).^2); wQD=1./(abs(mQD3P*sigma).^2);
wVolt(5)=0; % wVolt(5)=0;
wPD(5)=0; % wPD(5)=0;
wQD(5)=0; % wQD(5)=0;
%% %%
RestraintCount=size(Loadi,1)*2; %约束条件数,放开所有QD RestraintCount=size(Loadi,1)*2; %约束条件数
Init_Z=sparse(ones(1,RestraintCount)); Init_Z=sparse(ones(1,RestraintCount));
Init_W=sparse(-1*ones(1,RestraintCount)); Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=1*sparse(ones(1,RestraintCount)); Init_L=1*sparse(ones(1,RestraintCount));
@ -234,7 +234,7 @@ Init_Y=sparse(1,2*Busnum*3);%
Gap=(Init_L*Init_Z'-Init_U*Init_W'); Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
ContrlCount=size(Loadi,1)*2+Busnum*6; ContrlCount=size(Loadi,1)*2+Busnum*6;
kmax=900; kmax=20;
%% %%
Precision=1e-5; Precision=1e-5;
CenterA=0.1; CenterA=0.1;
@ -292,20 +292,19 @@ end
JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2); JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2);
%估计误差统计 %估计误差统计
%量测量数量 %量测量数量
Busnum=busNum;
mCount=Busnum*3+length(Loadi)*3*2; mCount=Busnum*3+length(Loadi)*3*2;
%估计量质量 %估计量质量
%有PD，QD的 AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(Volt))));
SEMeasurement=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((rPD3P-PD3P)./mPD3P./sigma).^2)+sum(((rQD3P-QD3P)./mQD3P./sigma).^2); AME_VAngle=sum(sum(abs( rVoltABCA-UAngel)));
SEMeasurement=(SEMeasurement/mCount)^.5; AME_PD=sum(sum(abs(rPD3P-PD3P)));
%没有PD，QD的 AME_QD=sum(sum(abs(rQD3P-QD3P)));
SEMeasurementNOPQ=sum(((rVoltABCV-Volt)./mVoltABCV./sigma).^2); %计算与量测值的
SEMeasurementNOPQ=(SEMeasurementNOPQ/mCount)^.5; AME_mVolt=sum(sum(abs( mVoltABCV-rVoltABCV)));
%量测量质量 AME_mPD=sum(sum(abs(rPD3P-mPD3P)));
SEMeasurementM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2)+sum(((rPD3P-mPD3P)./mPD3P./sigma).^2)+sum(((rQD3P-mQD3P)./mQD3P./sigma).^2); AME_mQD=sum(sum(abs(rQD3P-mQD3P)));
SEMeasurementM=(SEMeasurementM/mCount)^.5; %返回收敛信息
%没有PD，QD的 isConverged=1;
SEMeasurementNOPQM=sum(((rVoltABCV-mVoltABCV)./mVoltABCV./sigma).^2);
SEMeasurementNOPQM=(SEMeasurementNOPQM/mCount)^.5;
fprintf('迭代次数%d\n',KK); fprintf('迭代次数%d\n',KK);
toc toc
end end

18
feeder13/data - 副本.txt
View File

@ -1,16 +1,16 @@
650 650
0 0
603 632 645 500 600 632 645 500
602 632 633 500 602 632 633 500
-1 633 634 0 601 633 634 0
603 645 646 300 600 645 646 300
601 650 632 2000 601 650 632 2000
607 684 652 800 602 684 652 800
601 632 671 2000 601 632 671 2000
604 671 684 300 600 671 684 300
601 671 680 1000 602 671 680 1000
-2 671 692 0 600 671 692 0
605 684 611 300 602 684 611 300
606 692 675 500 601 692 675 500
0 0

46
feeder13/data1.txt
View File

@ -1,14 +1,50 @@
650 4.16 650 4.16
0 0
600 632 645 500
602 632 633 500
601 633 634 600
600 645 646 300
601 650 632 2000 601 650 632 2000
602 632 631 1800 602 684 652 800
602 631 630 1700 601 632 671 2000
600 671 684 300
602 671 680 1000
600 671 692 300
602 684 611 300
601 692 675 500
0 0
632 0 0 0 632 0 0 0
0 0
630 80 70 60 50 40 30 634 160 110 120 90 120 90
631 40 10 11 50 30 30 645 123 102 170 125 133 110
632 120 120 120 110 90 90 646 172 150 230 132 211 160
652 128 86 30 10 80 70
675 85 10 68 60 20 12
611 52 42 85 60 70 40
0
600 632 645 500
602 632 633 500
601 633 634 600
600 645 646 300
601 650 632 2000
602 684 652 800
601 632 671 2000
600 671 684 300
602 671 680 1000
600 671 692 300
602 684 611 300
601 692 675 500
0
632 0 0 0
0
634 160 110 120 90 120 90
645 123 102 170 125 133 110
646 172 150 230 132 0 211 160
652 128 86 30 10 80 70
671 385 220 385 220 385 220
675 485 190 68 60 290 212
692 142 140 153 121 170 151
611 156 142 185 60 170 80
0 0

14
feeder13/data1三相不平衡.txt Normal file
View File

@ -0,0 +1,14 @@
650 4.16
0
601 650 632 2000
602 632 631 1800
602 631 630 1700
0
632 0 0 0
0
630 90 70 80 40 30 10
631 40 10 60 30 100 80
632 12 12 120 120 11 11
0

10
func_deltG.m
View File

@ -8,18 +8,18 @@ sizeLoadi=size(Loadi,1);
%% %%
dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg4_dPD=sparse(size(Loadi,1),length(Loadi)); dg4_dPD=sparse(size(Loadi,1),length(Loadi));
dg5_dPD=sparse(size(Loadi,1),Busnum*3); % dg5_dPD=sparse(size(Loadi,1),Busnum*3);
%% %%
dg3_dQD=sparse(length(Loadi),length(Loadi)); dg3_dQD=sparse(length(Loadi),length(Loadi));
dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1)); dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
dg5_dQD=sparse(size(Loadi,1),Busnum*3); % dg5_dQD=sparse(size(Loadi,1),Busnum*3);
%% %%
dg3_dx=sparse(2*Busnum*3,sizeLoadi); dg3_dx=sparse(2*Busnum*3,sizeLoadi);
dg4_dx=sparse(2*Busnum*3,length(Loadi)); dg4_dx=sparse(2*Busnum*3,length(Loadi));
dg5_dx=[sparse(1:Busnum*3,1:Busnum*3,zeros(Busnum*3,1),Busnum*3,Busnum*3); % dg5_dx=[sparse(1:Busnum*3,1:Busnum*3,zeros(Busnum*3,1),Busnum*3,Busnum*3);
sparse(Busnum*3,Busnum*3); % sparse(Busnum*3,Busnum*3);
]; % ];
%% %%
% deltG=[dg3_dPD,dg4_dPD,dg5_dPD; % deltG=[dg3_dPD,dg4_dPD,dg5_dPD;

32
test.m
View File

@ -1,22 +1,26 @@
close all close all
clear clear
clc clc
AME_VoltSum=0;
AME_VAngleSum=0;
JMeasurementSum=0; JMeasurementSum=0;
SEMeasurementSum=0; AME_PDSum=0;
SEMeasurementNOPQSum=0; AME_QDSum=0;
SEMeasurementMSum=0;
SEMeasurementNOPQMSum=0;
N=1000; N=1000;
for I=1:N for I=1:N
[JMeasurement,SEMeasurement,SEMeasurementNOPQ,SEMeasurementM,SEMeasurementNOPQM]=OPF(); [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=OPF();
JMeasurementSum=JMeasurement+JMeasurementSum; if isConverged==0
SEMeasurementSum=SEMeasurement+SEMeasurementSum; I=I-1;
SEMeasurementNOPQSum=SEMeasurementNOPQ+SEMeasurementNOPQSum; continue
SEMeasurementMSum=SEMeasurementMSum+SEMeasurementM; end
SEMeasurementNOPQMSum=SEMeasurementNOPQMSum+SEMeasurementNOPQM; % [JMeasurement]=run();
AME_VoltSum=AME_VoltSum+AME_Volt;
AME_VAngleSum=AME_VAngleSum+AME_VAngle;
AME_PDSum=AME_PDSum+AME_PD;
AME_QDSum=AME_QDSum+AME_QD;
end end
JMeasurementSum=JMeasurementSum/N; JMeasurementSum=JMeasurementSum/N;
SEMeasurementSum=SEMeasurementSum/N; AME_VoltSum=AME_VoltSum/N;
SEMeasurementNOPQSum=SEMeasurementNOPQSum/N; AME_VAngleSum=AME_VAngleSum/N;
SEMeasurementMSum=SEMeasurementMSum/N; AME_PDSum=AME_PDSum/N;
SEMeasurementNOPQMSum=SEMeasurementNOPQMSum/N; AME_QDSum=AME_QDSum/N;