clc
clear
close all
arrayA=zeros(21,10);
sumCaseA_SE=0;
sumCaseB_SE=0;
sumCaseC_SE=0;
VoltAAE=0;
VAngleAAE=0;
for badDataNode=1:1
    loopN=1;
    maxDVolt=0;
    maxDVAngle=0;
    while 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+unifrnd(-RealPD(Loadi)*sigma,RealPD(Loadi)*sigma,length(Loadi),1));
        QD0(Loadi)=RealQD(Loadi).*(1+unifrnd(-RealQD(Loadi)*sigma,RealQD(Loadi)*sigma,length(Loadi),1));
        mVolt=rVolt.*(1+unifrnd(-rVolt*sigma,rVolt*sigma,1,length(rVolt)));
 
    
        
%         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(3)=rVolt(3)*(1-sigma*6);
        %% 画Case A的图
        % figure('Color',[1 1 1]);
        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
        if isConverge==0
            continue;
        end
        VoltAAE=VoltAAE+sum(abs((Volt-rVolt)./rVolt));
        VAngleAAE=VAngleAAE+sum(abs((UAngel(2:end)-rUAngel(2:end))./rUAngel(2:end)));
        loopN=loopN+1;
        if loopN>=500
            break;
        end
        
    end
    
end
VoltAAE=VoltAAE/(loopN*length(Volt))*100;
VAngleAAE=VAngleAAE/(loopN*length(UAngel(2:end)))*100;
%% PLOTING

% % 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);%Relative Error of Voltage in Case A
% % CaseAREV=CaseAREV(2:end)*100;
% %真实值的
% plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
% %测量值的
% % 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('节点号');
% ylabel('误差%');
% subplot(4,1,2);
% % CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
% CaseAREA=(UAngel-rUAngel);%Relative Error of Angle in Case A
% CaseAREA(1)=0;
% %真实值的
% 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=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A
% CaseAREP(1)=0;
% %真实值的
% plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond');
% %测量值的
% % 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('节点号');
% ylabel('误差%');
% subplot(4,1,4);
% % CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
% CaseAREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case A
% CaseAREQ(1)=0;
% %真实值的
% plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond');
% %测量值的
% % 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('节点号');
% 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));
% notZeros=find(PD0~=0);
% CaseA_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^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=(Volt-rVolt);%Relative Error of Voltage in Case B
% 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=(UAngel-rUAngel);%Relative Error of Angle in Case B
% 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=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case B
% RealPD(1)=0;
% 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=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case B
% CaseBREQ(1)=0;
% 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));
% 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(((RealPD(Measurei)-PD(Measurei))./(RealPD(Measurei)*1)).^2)+sum(((RealQD(Measurei)-QD(Measurei))./(RealQD(Measurei)*1)).^2)+sum(((Volt(Measurei)-rVolt(Measurei))./(rVolt(Measurei)*1)).^2);
% CaseB_SE=CaseB_SE+sum(((RealPD(noMeasurei)-PD(noMeasurei))./(RealPD(noMeasurei)*1)).^2)+sum(((RealQD(noMeasurei)-QD(noMeasurei))./(RealQD(noMeasurei)*1)).^2)+sum(((Volt(noMeasurei)-rVolt(noMeasurei))./(rVolt(noMeasurei)*1)).^2);
% CaseB_SE=(CaseB_SE/(length(notZeros)+length(noMeasurei) +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=(Volt-rVolt);
% 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=(UAngel-rUAngel);
% 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=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C
% CaseCREP(1);
% 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=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C
% CaseCREQ(1)=0;
% 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));
% 
% CaseC_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2);
% CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
% sumCaseC_SE=sumCaseC_SE+CaseC_SE;
% 
% % 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) %画出概率密度分布图
% voltBarCaseAX=x;
% voltBarCaseAY=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) %画出概率密度分布图
% voltBarCaseBX=x;
% voltBarCaseBY=yy;
% xlabel('Error');
% ylabel('Number of buses');
% title('算例B');
% % 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) %画出概率密度分布图
% voltBarCaseCX=x;
% voltBarCaseCY=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) %画出概率密度分布图
% angelBarCaseAX=x;
% angelBarCaseAY=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) %画出概率密度分布图
% angelBarCaseBX=x;
% angelBarCaseBY=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); %计算各个区间的个数
% angelBarCaseCX=x;
% angelBarCaseCY=yy;
% % 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) %画出概率密度分布图
% PDBarCaseAX=x;
% PDBarCaseAY=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) %画出概率密度分布图
% PDBarCaseBX=x;
% PDBarCaseBY=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) %画出概率密度分布图
% PDBarCaseCX=x;
% PDBarCaseCY=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) %画出概率密度分布图
% QDBarCaseAX=x;
% QDBarCaseAY=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) %画出概率密度分布图
% QDBarCaseBX=x;
% QDBarCaseBY=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) %画出概率密度分布图
% QDBarCaseCX=x;
% QDBarCaseCY=yy;
% %画收敛曲线
% fz=find(abs(plotGapA)==0);
% % fz=fz(1);
% figure('Name','互补曲线')
% plot(1:fz-1,plotGapA(1:fz-1));
% figure('Name','最大不平衡量');
% %% 最大不平衡量
% % maxDismatchPQ = [0.3123e-10  0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10];
% % bar(maxDismatchPQ);
% %% 计算时间
% calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;];
% bar(calTime,'stacked');
% % figure();
% % DeviationFigure(2:33,[CaseAREV(2:end);CaseBREV(2:end);CaseCREV(2:end)],[CaseAREA(2:end);CaseBREA(2:end);CaseCREA(2:end)],[CaseAREP(2:end),CaseBREP(2:end),CaseCREP(2:end)],[CaseAREQ(2:end),CaseBREQ(2:end),CaseCREQ(2:end)]);%,[CaseAREA;CaseBREA;CaseCREA],[CaseAREV;CaseBREV;CaseCREV],[CaseAREV;CaseBREV;CaseCREV]);
% % VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY);
% % AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY);
% % PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY);
% % QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY);
% % MaxErrorFigure()
% % MaxBoundErrorFigure();
% % DrawLoadProfile();

%% PLOTING