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