
fprintf('\n%s','=======================================================');
fprintf('\n%s\n','============SE OUTPUT=================================');
fprintf('\t%s%d%s%d%s\n','FOR ',nob,'-BUS & ',nol,'-BRANCH SYSTEM,');
if flag_sim==0
    fprintf('\t%s\n','WITH MEASUREMENTS AVAILABLE FROM THE FIELD DIRECTLY,');
else
    fprintf('\t%s\n','WITH MEASUREMENTS FROM SIMULATION,');
end
if flag_error==2
    fprintf('\t%s\n','USING UNIFORM ERROR DISTRIBUTION, ');
else
    fprintf('\t%s\n','USING GAUSSIAN ERROR DISTRIBUTION, ');
end

if flag_SE==2
    fprintf('\t%s\n','ALL PMU BASED LINEAR SE ');
%     fprintf('\t%s%d%s%d%s\n','FOR ',nob,'-BUS & ',nol,'-LINE SYSTEM');
    if flag_PMUSE==2
        fprintf('\t%s\n','CURRENT PHASOR CONVERTED INTO VOLTAGE PHASOR');
    else
        fprintf('\t%s\n','CONSIDER VOLTAGE & CURRENT PHASOR DIRECTLY');
    end
    
elseif flag_SE==3
    if flag_PMU==2
        if flag_eq==1
            fprintf('\t%s\n','SEQUENTIAL SE WITH EQUALITY CONSTRAINT ');
        else
            fprintf('\t%s\n','SEQUENTIAL SE WITHOUT EQUALITY CONSTRAINT ');
        end
    else
        if flag_eq==1
            fprintf('\t%s\n','HYBRID SE WITH EQUALITY CONSTRAINT ');
        else
            fprintf('\t%s\n','HYBRID SE WITHOUT EQUALITY CONSTRAINT');
        end
    end
%     fprintf('\n\t%s%d%s%d%s\n','FOR ',nob,'-BUS & ',nol,'-LINE SYSTEM');
else
    if flag_eq==1
        fprintf('\t%s\n','CONVENTIONAL SE WITH EQUALITY CONSTRAINT');
    else
        fprintf('\t%s\n','CONVENTIONAL SE WITHOUT EQUALITY CONSTRAINT');
    end
%     fprintf('\n\t%s%d%s%d%s\n','FOR ',nob,'-BUS & ',nol,'-LINE SYSTEM');
end
    


if flag_SE==3 || flag_SE==1
    if flag_SE_method==1
        fprintf('\t%s\n','NORMAL METHOD USING MATLAB \');
%     elseif flag_SE_method==5
%         fprintf('\t%s\n','LU FACTORIZATION BASED METHOD');
    elseif flag_SE_method==2
        fprintf('\t%s\n','PETERS WILKINSON METHOD');
    elseif flag_SE_method==3
        fprintf('\t%s\n','LEAST SQUARE OPTIMIZATION BASED, i.e. PSEUDO INVERSE');
    elseif flag_SE_method==4
        fprintf('\t%s\n','L-CURVE BASED TIKHONOV REGULARIZATION');
    else
        fprintf('\t%s\n','NORMAL METHOD USING MATLAB \');
    end
end
fprintf('%s\n','=======================================================');

if flag_SE==2
    x_est=x_est_PMU;
elseif flag_SE==3
    if flag_PMU==2
        x_est=x_est_II;
    else
        x_est=x_est_HSE;
    end
else
    x_est=x_est_conv;
end
disp('   V_abs V_angle_radian Vreal    Vimag');
disp(x_est);

% Estimated value of measurements related to bus
meas_bus_est=meas_bus;
VV=x_est(:,3)+1i*x_est(:,4);
S_calc=VV.*conj(YBUS*VV);
for i=1:length(meas_bus_est(:,1))
    if meas_bus_est(i,2)==1  %Real Power Injection
        meas_bus_est(i,1)=real(S_calc(meas_bus_est(i,3)));
    elseif meas_bus_est(i,2)==2  %Reactive Power Injection
        meas_bus_est(i,1)=imag(S_calc(meas_bus_est(i,3)));
    elseif meas_bus_est(i,2)==3  %Voltage Magnitude
        meas_bus_est(i,1)=abs(VV(meas_bus_est(i,3)));
    elseif meas_bus_est(i,2)==4  %Voltage Angle in radian
        meas_bus_est(i,1)=angle(VV(meas_bus_est(i,3)));
    end
end

% Estimated value of measurements related to branch
meas_line_est=meas_line;
VM=x_est(:,1);VA=x_est(:,2);
for i=1:length(meas_line_est(:,1))
    k=meas_line_est(i,3);
    if meas_line_est(i,4)~=line_pi(k,1)
        k=k+nol;
    end
    m=line_pi(k,1);
    n=line_pi(k,2);
    g_t=line_pi(k,5);
    b_t=line_pi(k,6);
    gi_t=line_pi(k,3);
    bi_t=line_pi(k,4);
    
    if meas_line_est(i,2)==1  %Real Power flow
        meas_line_est(i,1)=VM(m)^2*(g_t+gi_t)-VM(m)*VM(n)*(g_t*cos(VA(m)-VA(n))+b_t*sin(VA(m)-VA(n)));
    elseif meas_line_est(i,2)==2  %Reactive Power Flow
        meas_line_est(i,1)=-VM(m)^2*(b_t+bi_t)-VM(m)*VM(n)*(g_t*sin(VA(m)-VA(n))-b_t*cos(VA(m)-VA(n)));
    elseif meas_line_est(i,2)==3 %Current Magnitude
        Imnreal=VM(m)*((g_t+gi_t)*cos(VA(m))-(b_t+bi_t)*sin(VA(m)))-VM(n)*(g_t*cos(VA(n))-b_t*sin(VA(n)));
        Imnimag=VM(m)*((g_t+gi_t)*sin(VA(m))+(b_t+bi_t)*cos(VA(m)))-VM(n)*(g_t*sin(VA(n))+b_t*cos(VA(n)));
        meas_line_est(i,1)=abs(Imnreal+1i*Imnimag);
        clear Imnreal Imnimag;
    elseif meas_line_est(i,2)==4  %Current Angle in radian
        Imnreal=VM(m)*((g_t+gi_t)*cos(VA(m))-(b_t+bi_t)*sin(VA(m)))-VM(n)*(g_t*cos(VA(n))-b_t*sin(VA(n)));
        Imnimag=VM(m)*((g_t+gi_t)*sin(VA(m))+(b_t+bi_t)*cos(VA(m)))-VM(n)*(g_t*sin(VA(n))+b_t*cos(VA(n)));
        meas_line_est(i,1)=angle(Imnreal+1i*Imnimag);
        clear Imnreal Imnimag;
    elseif meas_line_est(i,2)==5  %Real Part of the Current
        meas_line_est(i,1)=VM(m)*((g_t+gi_t)*cos(VA(m))-(b_t+bi_t)*sin(VA(m)))-VM(n)*(g_t*cos(VA(n))-b_t*sin(VA(n)));
    elseif meas_line_est(i,2)==6  %Reactive Part of the Current
        meas_line_est(i,1)=VM(m)*((g_t+gi_t)*sin(VA(m))+(b_t+bi_t)*cos(VA(m)))-VM(n)*(g_t*sin(VA(n))+b_t*cos(VA(n)));
    end
end
xlswrite('meas_estimated.xlsx',meas_bus_est,'bus_meas');
xlswrite('meas_estimated.xlsx',meas_line_est,'line_meas');

z_err_bus=abs(meas_bus_est(:,1)-meas_bus(:,1));
z_err_line=abs(meas_line_est(:,1)-meas_line(:,1));
z_err=[z_err_bus;z_err_line];
x_idx=1:length(z_err);
plot(x_idx',z_err);
xlabel('Measurement index');
ylabel('absolute error in measurements');

if flag_sim==1
    plot_state_error;
end
clear VV VM VA S_calc i k m  n g_t b_t gi_t bi_t temp z_err_bus z_err_line;
clear z_err x_idx;
