clc
clear
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:\算例\17\17.csv');
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))';
%% 画Case A的图
figV=figure();
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
subplot(4,1,1);
plot(1:Busnum,Volt-rVolt,'b.:','Marker','diamond');
subplot(4,1,2);
plot(1:Busnum,UAngel-rUAngel,'b:','Marker','diamond');
subplot(4,1,3);
plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'b:','Marker','diamond');
subplot(4,1,4);
plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'b:','Marker','diamond');
%% Case B
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,5,11,12,13,16],[18,19]),PD0,QD0,mVolt,sigma);%
% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,6,7],[]),PD0,QD0,mVolt,sigma);%
subplot(4,1,1);
hold on;
plot(1:Busnum,Volt-rVolt,'g.:','Marker','square');
subplot(4,1,2);
hold on;
plot(1:Busnum,UAngel-rUAngel,'g:','Marker','square');
subplot(4,1,3);
hold on;
plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'g:','Marker','square');
subplot(4,1,4);
hold on;
plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'g:','Marker','square');
%% Case C
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([1:17],PD0,QD0,mVolt,sigma);%
subplot(4,1,1);
hold on;
plot(1:Busnum,Volt-rVolt,'r.:','Marker','o');
subplot(4,1,2);
hold on;
plot(1:Busnum,UAngel-rUAngel,'r:','Marker','o');
subplot(4,1,3);
hold on;
plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'r:','Marker','o');
subplot(4,1,4);
hold on;
plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'r:','Marker','o');
% 画legend
subplot(4,1,1);
title('Voltage');
legend('Case A','Case B','Case C');
subplot(4,1,2);
title('Voltage Angle');
legend('Case A','Case B','Case C');
subplot(4,1,3);
title('Active load power');
legend('Case A','Case B','Case C');
subplot(4,1,4);
title('Reactive load  power');
legend('Case A','Case B','Case C');
obj=sum(Vbi)+sum(PDbi)+sum(QDbi);
fprintf('目标函数值 %.2f\n',full(obj));