pes2014-debug2-laplace/OPF.m

clc
clear
close all
arrayA=zeros(21,10);
sumCaseA_SE=0;
sumCaseB_SE=0;
sumCaseC_SE=0;
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:\<EFBFBD><EFBFBD><EFBFBD><EFBFBD>\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))';

%<EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ı<EFBFBD><EFBFBD><EFBFBD>
% PD0=load('PD0');
% PD0=PD0.PD0;
% QD0=load('QD0');
% QD0=QD0.QD0;
% mVolt=load('mVolt');
% mVolt=mVolt.mVolt;

mVolt(2)=rVolt(2)*(1+sigma*8);
%% <EFBFBD><EFBFBD>Case A<EFBFBD><EFBFBD>ͼ
% figure('Color',[1 1 1]);
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%ȫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
end
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;
%<EFBFBD><EFBFBD>ʵֵ<EFBFBD><EFBFBD>
plot(1:length(CaseAREV),abs(CaseAREV),'k.:','Marker','diamond');
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','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('<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>');
ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>%');
subplot(4,1,2);
CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
% CaseAREA=CaseAREA(2:end)*100;
%<EFBFBD><EFBFBD>ʵֵ<EFBFBD><EFBFBD>
plot(1:length(CaseAREA),abs(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('<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>');
ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>%');
subplot(4,1,3);
CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
% CaseAREP=CaseAREP(2:end)*100;
%<EFBFBD><EFBFBD>ʵֵ<EFBFBD><EFBFBD>
plot(1:length(CaseAREP),abs(CaseAREP),'k:','Marker','diamond');
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','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('<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>');
ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>%');
subplot(4,1,4);
CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
% CaseAREQ=CaseAREQ(2:end)*100;
%<EFBFBD><EFBFBD>ʵֵ<EFBFBD><EFBFBD>
plot(1:length(CaseAREQ),abs(CaseAREQ),'k:','Marker','diamond');
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','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('<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>');
ylabel('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>%');
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Case A<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
notZeros=find(PD0~=0);
CaseA_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseA_SE=sumCaseA_SE+CaseA_SE;
% arrayA(1:19,I)=Vbi;
% arrayA(21,I)=CaseAE*1000;
%% 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);%
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),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),abs(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),abs(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),abs(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),abs(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));
noMeasurei=[2,3,5,20,24,27,28,10,11,12,13];
Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]);
CaseB_SE=sum(((PD0(Measurei)-PD(Measurei))./(PD0(Measurei)*sigma)).^2)+sum(((QD0(Measurei)-QD(Measurei))./(QD0(Measurei)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseB_SE=CaseB_SE+sum(((PD0(noMeasurei)-PD(noMeasurei))./(PD0(noMeasurei)*20*sigma)).^2)+sum(((QD0(noMeasurei)-QD(noMeasurei))./(QD0(noMeasurei)*20*sigma)).^2);
CaseB_SE=(CaseB_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseB_SE=sumCaseB_SE+CaseB_SE;
%% Case C
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],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),abs(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),abs(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),abs(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),abs(CaseCREQ),'r:','Marker','o');
% <EFBFBD><EFBFBD>legend
subplot(4,1,1);
% title('Voltage');
ld=legend('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>A','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>C');
set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,2);
% title('Voltage Angle');
ld=legend('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>A','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>C');
set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,3);
% title('Active load power');
ld=legend('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>A','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>C');
set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
subplot(4,1,4);
% title('Reactive load  power');
ld=legend('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>A','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>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));

CaseC_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseC_SE=sumCaseC_SE+CaseC_SE;

% fprintf('Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD>ֵ %.2f\n',full(obj));
fprintf('Case A    Case B    Case C <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n')
fprintf('%f    %f    %f    \n',CaseAE,CaseBE,CaseCE);
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>CaseĿ<EFBFBD><EFBFBD>ֵ\n')
fprintf('%f\t%f\t%f \n',objA,objB,objC)
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% 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');
%% ֱ<EFBFBD><EFBFBD>ͼ
% <EFBFBD><EFBFBD>ѹ
figure('Name','<EFBFBD><EFBFBD>ѹֱ<EFBFBD><EFBFBD>ͼ')
split_number=20;
%Case A
subplot(1,3,1)
y=CaseAREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
xlabel('Error');
ylabel('Number of buses');
title('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>A');
% ylim([0 4])
%Case B
subplot(1,3,2)
y=CaseBREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
xlabel('Error');
ylabel('Number of buses');
title('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>D');
% ylim([0 4])
%Case C
subplot(1,3,3)
y=CaseCREV;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
xlabel('Error');
ylabel('Number of buses');
title('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>C');
% ylim([0 4])
figure('Name','<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD>ͼ')
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREA;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
%PD
figure('Name','<EFBFBD>й<EFBFBD>ֱ<EFBFBD><EFBFBD>ͼ')
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREP;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
%QD
figure('Name','<EFBFBD>޹<EFBFBD>ֱ<EFBFBD><EFBFBD>ͼ')
split_number=20;
%Case A
subplot(2,2,1)
y=CaseAREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case B
subplot(2,2,2)
y=CaseBREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
% ylim([0 4])
%Case C
subplot(2,2,3)
y=CaseCREQ;
ymin=min(y);
ymax=max(y);
x=linspace(ymin,ymax,split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֳ<EFBFBD>split_number<EFBFBD><EFBFBD><EFBFBD>ȷֵ<EFBFBD>(19<EFBFBD>ȷ<EFBFBD>),Ȼ<EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
yy=hist(y,x); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĸ<EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD>ʽΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺͼ<EFBFBD>:
bar(x,yy) %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܶȷֲ<EFBFBD>ͼ
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
fz=find(abs(plotGapA)==0);
fz=fz(1);
figure('Name','<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>')
plot(1:fz-1,plotGapA(1:fz-1));
figure('Name','<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>');
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% maxDismatchPQ = [0.3123e-10  0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10];
% bar(maxDismatchPQ);
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
calTime=[87 61; 211 87; 108 59;];
bar(calTime,'stacked');