facat
/
distributionnetwork-power2current
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

1.添加了画柱状图的函数 

2.加了结果展示代码

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2014-11-03 11:07:25 +08:00
parent 971d2e3455
commit 0c9c2e2eb7
2 changed files with 54 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
barPlot.m Normal file
View File

@ -0,0 +1,13 @@
function barPlot( data, split_number,xWord,yWord,figTitle)
ymin=min(data);
ymax=max(data);
% split_number=10;
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
yy=hist(data,x); %计算各个区间的个数
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
bar(x,yy) %画出概率密度分布图
xlabel(xWord);
ylabel(yWord);
title(figTitle);
end

50
run.m
View File

@ -1,5 +1,6 @@
%% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。 %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。
% %
close all
clc clc
clear clear
lineZ=readLineZ('feeder13\lineParameter.txt'); lineZ=readLineZ('feeder13\lineParameter.txt');
@ -197,8 +198,8 @@ sigmaI1=normrnd(0,sigma,length(Loadi),1);
I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流 I1measurement=mIf1(Loadi).*(1+sigmaI1);%测量值是等效发电机电流
wI1r=abs( real(mIf1(Loadi).*sigmaI1) ); wI1r=abs( real(mIf1(Loadi).*sigmaI1) );
wI1i=abs( imag(mIf1(Loadi).*sigmaI1) ); wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
[ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 ); % [ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2); % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(rV1)-V1r).^2)+sum((imag(rV1)-V1i).^2);
% fprintf('目标值 %f\n',full(f)); % fprintf('目标值 %f\n',full(f));
%% 算负序的 %% 算负序的
fprintf('负序\n'); fprintf('负序\n');
@ -219,8 +220,8 @@ sigmaI2=normrnd(0,sigma,length(Loadi),1);
I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流 I2measurement=mIf2(Loadi).*(1+sigmaI2);%测量值是等效发电机电流
wI2r=abs( real(mIf2(Loadi).*sigmaI2) ); wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
wI2i=abs( imag(mIf2(Loadi).*sigmaI2) ); wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
[ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 ); % [ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2); % f=sum(([real(I2measurement);imag(I2measurement)]-[I2r;I2i]).^2)+sum((real(rV2)-V2r).^2)+sum((imag(rV2)-V2i).^2);
% fprintf('目标值 %f\n',full(f)); % fprintf('目标值 %f\n',full(f));
%% 算零序 %% 算零序
fprintf('零序\n'); fprintf('零序\n');
@ -241,8 +242,25 @@ sigmaI0=normrnd(0,sigma,length(Loadi),1);
I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流 I0measurement=mIf0(Loadi).*(1+sigmaI0);%测量值是等效发电机电流
wI0r=abs( real(mIf0(Loadi).*sigmaI0) ); wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
wI0i=abs( imag(mIf0(Loadi).*sigmaI0) ); wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
[ V0r,V0i,I0r,I0i ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 ); % matlabpool local 3
f=sum(([real(I0measurement);imag(I0measurement)]-[I0r;I0i]).^2)+sum((real(rV0)-V0r).^2)+sum((imag(rV0)-V0i).^2); tic
for II=1:3
if II==1
[ V1r,V1i,I1r,I1i ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
end
if II==2
[ V2r,V2i,I2r,I2i ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
end
if II==3
[ V0r,V0i,I0r,I0i ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 );
end
end
toc
% matlabpool close
% f=sum(([real(I0measurement);imag(I0measurement)]-[I0r;I0i]).^2)+sum((real(rV0)-V0r).^2)+sum((imag(rV0)-V0i).^2);
% fprintf('目标值 %f\n',full(f)); % fprintf('目标值 %f\n',full(f));
%% 把三序合成三相 %% 把三序合成三相
% 三相电压 % 三相电压
@ -257,7 +275,21 @@ SEIf1=I1r+1j*I1i;
SEIf2=I2r+1j*I2i; SEIf2=I2r+1j*I2i;
SEIpABC=full(Tf2p*conj([SEIf0';SEIf1(1:3)';SEIf2'])); SEIpABC=full(Tf2p*conj([SEIf0';SEIf1(1:3)';SEIf2']));
%看一下差多少 %看一下差多少
(full(abs(VoltpABC))-abs(SEVoltpABC))./abs(SEVoltpABC)*100 VError=(full(abs(VoltpABC))-abs(SEVoltpABC))./abs(VoltpABC)*100;
VError=reshape(VError,size(VError,1)*size(VError,2),1);
barPlot( VError, 10,['相对误差%'],['分布密度'],['电压幅值']);
% 三相负荷 % 三相负荷
SEVoltpABC(:,setxor(1:size(SEVoltpABC,2),Balance)).*conj(-SEIpABC) rThreeLoad=[ phaseASpotLoadP'+1j*phaseASpotLoadQ';
phaseBSpotLoadP'+1j*phaseBSpotLoadQ';
phaseCSpotLoadP'+1j*phaseCSpotLoadQ';
];
rThreeLoad=rThreeLoad(:,setxor(1:size(SEVoltpABC,2),Balance));
SEThreeLoad=SEVoltpABC(:,setxor(1:size(SEVoltpABC,2),Balance)).*conj(-SEIpABC);
phaseLoadPError=real(rThreeLoad-SEThreeLoad)./real(rThreeLoad)*100;
phaseLoadQError=imag(rThreeLoad-SEThreeLoad)./imag(rThreeLoad)*100;
phaseLoadPError=reshape(phaseLoadPError,size(phaseLoadPError,1)*size(phaseLoadPError,2),1);
phaseLoadQError=reshape(phaseLoadQError,size(phaseLoadQError,1)*size(phaseLoadQError,2),1);
figure()
barPlot( phaseLoadPError, 10,['相对误差%'],['分布密度'],['有功负荷误差']);
figure()
barPlot( phaseLoadQError, 10,'相对误差%','分布密度','无功负荷误差');