0001 function [V, converged, iterNum, z, z_est, error_sqrsum] = doSE(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0016 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0017 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
0018 RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;
0019 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
0020 GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
0021
0022
0023 tol = 1e-5;
0024 max_it = 100;
0025 verbose = 0;
0026
0027
0028 j = sqrt(-1);
0029 converged = 0;
0030 i = 0;
0031 V = V0;
0032 Va = angle(V);
0033 Vm = abs(V);
0034
0035 nb = size(Ybus, 1);
0036 f = branch(:, F_BUS);
0037 t = branch(:, T_BUS);
0038
0039
0040 nonref = [pv;pq];
0041
0042
0043 z = [
0044 measure.PF
0045 measure.PT
0046 measure.PG
0047 measure.Va
0048 measure.QF
0049 measure.QT
0050 measure.QG
0051 measure.Vm
0052 ];
0053
0054
0055 idx_zPF = idx.idx_zPF;
0056 idx_zPT = idx.idx_zPT;
0057 idx_zPG = idx.idx_zPG;
0058 idx_zVa = idx.idx_zVa;
0059 idx_zQF = idx.idx_zQF;
0060 idx_zQT = idx.idx_zQT;
0061 idx_zQG = idx.idx_zQG;
0062 idx_zVm = idx.idx_zVm;
0063
0064
0065 sigma_vector = [
0066 sigma.sigma_PF*ones(size(idx_zPF, 1), 1)
0067 sigma.sigma_PT*ones(size(idx_zPT, 1), 1)
0068 sigma.sigma_PG*ones(size(idx_zPG, 1), 1)
0069 sigma.sigma_Va*ones(size(idx_zVa, 1), 1)
0070 sigma.sigma_QF*ones(size(idx_zQF, 1), 1)
0071 sigma.sigma_QT*ones(size(idx_zQT, 1), 1)
0072 sigma.sigma_QG*ones(size(idx_zQG, 1), 1)
0073 sigma.sigma_Vm*ones(size(idx_zVm, 1), 1)
0074 ];
0075 sigma_square = sigma_vector.^2;
0076 R_inv = diag(1./sigma_square);
0077
0078
0079 while (~converged & i < max_it)
0080
0081 i = i + 1;
0082
0083
0084 Sfe = V(f) .* conj(Yf * V);
0085 Ste = V(t) .* conj(Yt * V);
0086
0087 gbus = gen(:, GEN_BUS);
0088 Sgbus = V(gbus) .* conj(Ybus(gbus, :) * V);
0089 Sgen = Sgbus * baseMVA + (bus(gbus, PD) + j*bus(gbus, QD));
0090 Sgen = Sgen/baseMVA;
0091 z_est = [
0092 real(Sfe(idx_zPF));
0093 real(Ste(idx_zPT));
0094 real(Sgen(idx_zPG));
0095 angle(V(idx_zVa));
0096 imag(Sfe(idx_zQF));
0097 imag(Ste(idx_zQT));
0098 imag(Sgen(idx_zQG));
0099 abs(V(idx_zVm));
0100 ];
0101
0102
0103 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0104 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0105
0106 genbus_row = gbus;
0107
0108
0109 dPF_dVa = real(dSf_dVa);
0110 dQF_dVa = imag(dSf_dVa);
0111 dPF_dVm = real(dSf_dVm);
0112 dQF_dVm = imag(dSf_dVm);
0113 dPT_dVa = real(dSt_dVa);
0114 dQT_dVa = imag(dSt_dVa);
0115 dPT_dVm = real(dSt_dVm);
0116 dQT_dVm = imag(dSt_dVm);
0117
0118 dPG_dVa = real(dSbus_dVa(genbus_row, :));
0119 dQG_dVa = imag(dSbus_dVa(genbus_row, :));
0120 dPG_dVm = real(dSbus_dVm(genbus_row, :));
0121 dQG_dVm = imag(dSbus_dVm(genbus_row, :));
0122
0123 dVa_dVa = eye(nb);
0124 dVa_dVm = zeros(nb, nb);
0125
0126 dVm_dVa = zeros(nb, nb);
0127 dVm_dVm = eye(nb);
0128 H = [
0129 dPF_dVa(idx_zPF, nonref) dPF_dVm(idx_zPF, nonref);
0130 dPT_dVa(idx_zPT, nonref) dPT_dVm(idx_zPT, nonref);
0131 dPG_dVa(idx_zPG, nonref) dPG_dVm(idx_zPG, nonref);
0132 dVa_dVa(idx_zVa, nonref) dVa_dVm(idx_zVa, nonref);
0133 dQF_dVa(idx_zQF, nonref) dQF_dVm(idx_zQF, nonref);
0134 dQT_dVa(idx_zQT, nonref) dQT_dVm(idx_zQT, nonref);
0135 dQG_dVa(idx_zQG, nonref) dQG_dVm(idx_zQG, nonref);
0136 dVm_dVa(idx_zVm, nonref) dVm_dVm(idx_zVm, nonref);
0137 ];
0138
0139
0140 J = H'*R_inv*H;
0141 F = H'*R_inv*(z-z_est);
0142 if ~isobservable(H, pv, pq)
0143 error('doSE: system is not observable');
0144 end
0145 dx = (J \ F);
0146
0147
0148 normF = norm(F, inf);
0149 if verbose > 1
0150 fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF);
0151 end
0152 if normF < tol
0153 converged = 1;
0154 end
0155
0156
0157 Va(nonref) = Va(nonref) + dx(1:size(nonref, 1));
0158 Vm(nonref) = Vm(nonref) + dx(size(nonref, 1)+1:2*size(nonref, 1));
0159 V = Vm .* exp(j * Va);
0160 Vm = abs(V);
0161 Va = angle(V);
0162 end
0163
0164 iterNum = i;
0165
0166
0167 error_sqrsum = sum((z - z_est).^2./sigma_square);