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