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