293 lines
10 KiB
Matlab
293 lines
10 KiB
Matlab
clc
|
|
clear
|
|
close all
|
|
arrayA=zeros(21,10);
|
|
for I=1:1
|
|
close
|
|
[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:\算例\feeder33\feeder33.txt');
|
|
sigma=0.01;
|
|
RealPD=PD;
|
|
RealQD=QD;
|
|
rVolt=Volt;
|
|
Loadi=find(PD~=0);
|
|
PD0=sparse(Busnum,1);
|
|
QD0=sparse(Busnum,1);
|
|
PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
|
|
QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
|
|
mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
|
|
|
|
%加载保存的变量
|
|
% PD0=load('PD0');
|
|
% PD0=PD0.PD0;
|
|
% QD0=load('QD0');
|
|
% QD0=QD0.QD0;
|
|
% mVolt=load('mVolt');
|
|
% mVolt=mVolt.mVolt;
|
|
|
|
% mVolt(19)=rVolt(19)*(1+sigma*4);
|
|
%% 画Case A的图
|
|
figure('Color',[1 1 1]);
|
|
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
|
|
subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]);
|
|
CaseAREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case A
|
|
% CaseAREV=CaseAREV(2:end)*100;
|
|
plot(1:length(CaseAREV),CaseAREV,'k.:','Marker','diamond');
|
|
box off;
|
|
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
|
|
xlabel('节点号');
|
|
ylabel('误差%');
|
|
subplot(4,1,2);
|
|
CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
|
|
% CaseAREA=CaseAREA(2:end)*100;
|
|
plot(1:length(CaseAREA),CaseAREA,'k:','Marker','diamond');
|
|
box off;
|
|
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
|
|
xlabel('节点号');
|
|
ylabel('误差%');
|
|
subplot(4,1,3);
|
|
CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
|
|
% CaseAREP=CaseAREP(2:end)*100;
|
|
plot(1:length(CaseAREP),CaseAREP,'k:','Marker','diamond');
|
|
box off;
|
|
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
|
|
xlabel('节点号');
|
|
ylabel('误差%');
|
|
subplot(4,1,4);
|
|
CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
|
|
% CaseAREQ=CaseAREQ(2:end)*100;
|
|
plot(1:length(CaseAREQ),CaseAREQ,'k:','Marker','diamond');
|
|
box off;
|
|
set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
|
|
xlabel('节点号');
|
|
ylabel('误差%');
|
|
%计算Case A的误差
|
|
CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
|
|
objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
|
|
|
|
% arrayA(1:19,I)=Vbi;
|
|
% arrayA(21,I)=CaseAE*1000;
|
|
|
|
end
|
|
%% Case B
|
|
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);%
|
|
subplot(4,1,1);
|
|
hold on;
|
|
CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
|
|
% CaseBREV=CaseBREV(2:end)*100;
|
|
plot(1:length(CaseBREV),CaseBREV,'b.:','Marker','square');
|
|
subplot(4,1,2);
|
|
hold on;
|
|
CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
|
|
% CaseBREA=CaseBREA(2:end)*100;
|
|
plot(1:length(CaseBREA),CaseBREA,'b:','Marker','square');
|
|
subplot(4,1,3);
|
|
hold on;
|
|
CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
|
|
% CaseBREP=CaseBREP(2:end)*100;
|
|
plot(1:length(CaseBREP),CaseBREP,'b:','Marker','square');
|
|
subplot(4,1,4);
|
|
hold on;
|
|
CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
|
|
% CaseBREQ=CaseBREQ(2:end)*100;
|
|
plot(1:length(CaseBREQ),CaseBREQ,'b:','Marker','square');
|
|
CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
|
|
objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
|
|
%% Case C
|
|
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:17],PD0,QD0,mVolt,sigma);%
|
|
subplot(4,1,1);
|
|
hold on;
|
|
CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
|
|
% CaseCREV=CaseCREV(2:end)*100;
|
|
plot(1:length(CaseCREV),CaseCREV,'r.:','Marker','o');
|
|
subplot(4,1,2);
|
|
hold on;
|
|
CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
|
|
% CaseCREA=CaseCREA(2:end)*100;
|
|
plot(1:length(CaseCREA),CaseCREA,'r:','Marker','o');
|
|
subplot(4,1,3);
|
|
hold on;
|
|
CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
|
|
% CaseCREP=CaseCREP(2:end)*100;
|
|
plot(1:length(CaseCREP),CaseCREP,'r:','Marker','o');
|
|
subplot(4,1,4);
|
|
hold on;
|
|
CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
|
|
% CaseCREQ=CaseCREQ(2:end)*100;
|
|
plot(1:length(CaseCREQ),CaseCREQ,'r:','Marker','o');
|
|
% 画legend
|
|
subplot(4,1,1);
|
|
% title('Voltage');
|
|
ld=legend('算例A','算例B','算例C');
|
|
set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
|
|
subplot(4,1,2);
|
|
% title('Voltage Angle');
|
|
ld=legend('算例A','算例B','算例C');
|
|
set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]);
|
|
subplot(4,1,3);
|
|
% title('Active load power');
|
|
ld=legend('算例A','算例B','算例C');
|
|
set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
|
|
subplot(4,1,4);
|
|
% title('Reactive load power');
|
|
ld=legend('算例A','算例B','算例C');
|
|
set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]);
|
|
CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
|
|
objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
|
|
% fprintf('目标函数值 %.2f\n',full(obj));
|
|
fprintf('Case A Case B Case C 的误差\n')
|
|
fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE);
|
|
fprintf('三个Case目标值\n')
|
|
fprintf('%f\t%f\t%f \n',objA,objB,objC)
|
|
%% 画测量值
|
|
% subplot(4,1,1);
|
|
% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')
|
|
% subplot(4,1,3);
|
|
% plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram');
|
|
% subplot(4,1,4);
|
|
% plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram');
|
|
%% 直方图
|
|
% 电压
|
|
figure('Name','电压直方图')
|
|
split_number=20;
|
|
%Case A
|
|
subplot(1,3,1)
|
|
y=CaseAREV;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
xlabel('Error');
|
|
ylabel('Number of buses');
|
|
title('算例A');
|
|
% ylim([0 4])
|
|
%Case B
|
|
subplot(1,3,2)
|
|
y=CaseBREV;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
xlabel('Error');
|
|
ylabel('Number of buses');
|
|
title('算例D');
|
|
% ylim([0 4])
|
|
%Case C
|
|
subplot(1,3,3)
|
|
y=CaseCREV;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
xlabel('Error');
|
|
ylabel('Number of buses');
|
|
title('算例C');
|
|
% ylim([0 4])
|
|
figure('Name','相角直方图')
|
|
% 相角
|
|
split_number=20;
|
|
%Case A
|
|
subplot(2,2,1)
|
|
y=CaseAREA;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case B
|
|
subplot(2,2,2)
|
|
y=CaseBREA;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case C
|
|
subplot(2,2,3)
|
|
y=CaseCREA;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
%PD
|
|
figure('Name','有功直方图')
|
|
split_number=20;
|
|
%Case A
|
|
subplot(2,2,1)
|
|
y=CaseAREP;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case B
|
|
subplot(2,2,2)
|
|
y=CaseBREP;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case C
|
|
subplot(2,2,3)
|
|
y=CaseCREP;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
%QD
|
|
figure('Name','无功直方图')
|
|
split_number=20;
|
|
%Case A
|
|
subplot(2,2,1)
|
|
y=CaseAREQ;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case B
|
|
subplot(2,2,2)
|
|
y=CaseBREQ;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
% ylim([0 4])
|
|
%Case C
|
|
subplot(2,2,3)
|
|
y=CaseCREQ;
|
|
ymin=min(y);
|
|
ymax=max(y);
|
|
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
|
|
yy=hist(y,x); %计算各个区间的个数
|
|
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
|
|
bar(x,yy) %画出概率密度分布图
|
|
%画收敛曲线
|
|
fz=find(abs(plotGapA)==0);
|
|
fz=fz(1);
|
|
figure('Name','互补曲线')
|
|
plot(1:fz-1,plotGapA(1:fz-1));
|