0001 function [V, converged, i] = state_est(branch, Ybus, Yf, Yt, Sbus, V0, ref, pv, pq, mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0017 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0018 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0019
0020
0021 if nargin < 10
0022 mpopt = mpoption;
0023 end
0024
0025
0026 tol = mpopt.pf.tol;
0027 max_it = mpopt.pf.nr.max_it;
0028
0029
0030 converged = 0;
0031 i = 0;
0032 nb = length(V0);
0033 nbr = size(Yf, 1);
0034 nref = [pv;pq];
0035 f = branch(:, F_BUS);
0036 t = branch(:, T_BUS);
0037
0038
0039 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V0);
0040 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V0);
0041 H = [
0042 real(dSf_dVa) real(dSf_dVm);
0043 real(dSt_dVa) real(dSt_dVm);
0044 real(dSbus_dVa) real(dSbus_dVm);
0045 speye(nb) sparse(nb,nb);
0046 imag(dSf_dVa) imag(dSf_dVm);
0047 imag(dSt_dVa) imag(dSt_dVm);
0048 imag(dSbus_dVa) imag(dSbus_dVm);
0049 sparse(nb,nb) speye(nb);
0050 ];
0051
0052
0053 z = [
0054 real(Sf);
0055 real(St);
0056 real(Sbus);
0057 angle(V0);
0058 imag(Sf);
0059 imag(St);
0060 imag(Sbus);
0061 abs(V0);
0062 ];
0063
0064
0065 fullscale = 30;
0066 sigma = [
0067 0.02 * abs(Sf) + 0.0052 * fullscale * ones(nbr,1);
0068 0.02 * abs(St) + 0.0052 * fullscale * ones(nbr,1);
0069 0.02 * abs(Sbus) + 0.0052 * fullscale * ones(nb,1);
0070 0.2 * pi/180 * 3*ones(nb,1);
0071 0.02 * abs(Sf) + 0.0052 * fullscale * ones(nbr,1);
0072 0.02 * abs(St) + 0.0052 * fullscale * ones(nbr,1);
0073 0.02 * abs(Sbus) + 0.0052 * fullscale * ones(nb,1);
0074 0.02 * abs(V0) + 0.0052 * 1.1 * ones(nb,1);
0075 ] ./ 3;
0076 ns = length(sigma);
0077 W = sparse(1:ns, 1:ns , sigma .^ 2, ns, ns );
0078 WInv = sparse(1:ns, 1:ns , 1 ./ sigma .^ 2, ns, ns );
0079
0080
0081
0082
0083
0084 err = normrnd( zeros(size(sigma)), sigma );
0085
0086
0087
0088
0089 z = z + err;
0090
0091
0092 V = ones(nb,1);
0093
0094
0095 Sfe = V(f) .* conj(Yf * V);
0096 Ste = V(t) .* conj(Yt * V);
0097 Sbuse = V .* conj(Ybus * V);
0098 z_est = [
0099 real(Sfe);
0100 real(Ste);
0101 real(Sbuse);
0102 angle(V);
0103 imag(Sfe);
0104 imag(Ste);
0105 imag(Sbuse);
0106 abs(V);
0107 ];
0108
0109
0110 delz = z - z_est;
0111 normF = delz' * WInv * delz;
0112 chusqu = err' * WInv * err;
0113
0114
0115 if mpopt.verbose > 1
0116 fprintf('\n it norm( F ) step size');
0117 fprintf('\n---- -------------- --------------');
0118 fprintf('\n%3d %10.3e %10.3e', i, normF, 0);
0119 end
0120 if normF < tol
0121 converged = 1;
0122 if mpopt.verbose > 1
0123 fprintf('\nConverged!\n');
0124 end
0125 end
0126
0127
0128
0129
0130
0131 vv=[(3:nbr), ...
0132 (nbr+1:2*nbr), ...
0133 (2*nbr+2:2*nbr+nb), ...
0134 (2*nbr+nb+2:2*nbr+2*nb), ...
0135 (2*nbr+2*nb+3:3*nbr+2*nb), ...
0136 (3*nbr+2*nb+1:4*nbr+2*nb), ...
0137 (4*nbr+2*nb+2:4*nbr+3*nb), ...
0138 (4*nbr+3*nb+2:4*nbr+4*nb)]';
0139
0140 ww = [ nref; nb+nref ];
0141
0142
0143 one_at_a_time = 1; max_it_bad_data = 50;
0144
0145 ibd = 1;
0146 while (~converged && ibd <= max_it_bad_data)
0147 nm = length(vv);
0148 baddata = 0;
0149
0150
0151 HH = H(vv,ww);
0152 WWInv = WInv(vv,vv);
0153 ddelz = delz(vv);
0154 VVa = angle(V(nref));
0155 VVm = abs(V(nref));
0156
0157
0158
0159
0160
0161
0162
0163 max_it = 100;
0164 i = 0;
0165 while (~converged && i < max_it)
0166
0167 i = i + 1;
0168
0169
0170 F = HH' * WWInv * ddelz;
0171 J = HH' * WWInv * HH;
0172 dx = (J \ F);
0173
0174
0175 VVa = VVa + dx(1:nb-1);
0176 VVm = VVm + dx(nb:2*nb-2);
0177 V(nref) = VVm .* exp(1j * VVa);
0178
0179
0180 Sfe = V(f) .* conj(Yf * V);
0181 Ste = V(t) .* conj(Yt * V);
0182 Sbuse = V .* conj(Ybus * V);
0183 z_est = [
0184 real(Sfe);
0185 real(Ste);
0186 real(Sbuse);
0187 angle(V);
0188 imag(Sfe);
0189 imag(Ste);
0190 imag(Sbuse);
0191 abs(V);
0192 ];
0193
0194
0195 delz = z - z_est;
0196 ddelz = delz(vv);
0197 normF = ddelz' * WWInv * ddelz;
0198
0199
0200 step = dx' * dx;
0201 if mpopt.verbose > 1
0202 fprintf('\n%3d %10.3e %10.3e', i, normF, step);
0203 end
0204 if (step < tol)
0205 converged = 1;
0206 if mpopt.verbose
0207 fprintf('\nState estimator converged in %d iterations.\n', i);
0208 end
0209 end
0210 end
0211 if mpopt.verbose
0212 if ~converged
0213 fprintf('\nState estimator did not converge in %d iterations.\n', i);
0214 end
0215 end
0216
0217
0218 B = zeros(nm,1);
0219 bad_threshold = 6.25;
0220 RR = inv(WWInv) - 0.95 * HH * inv(HH' * WWInv * HH) * HH';
0221
0222
0223
0224 rr = diag(RR);
0225
0226 B = ddelz .^ 2 ./ rr;
0227 [maxB,i_maxB] = max(B);
0228 if one_at_a_time
0229 if maxB >= bad_threshold
0230 rejected = i_maxB;
0231 else
0232 rejected = [];
0233 end
0234 else
0235 rejected = find( B >= bad_threshold );
0236 end
0237 if length(rejected)
0238 baddata = 1;
0239 converged = 0;
0240 if mpopt.verbose
0241 fprintf('\nRejecting %d measurement(s) as bad data:\n', length(rejected));
0242 fprintf('\tindex\t B\n');
0243 fprintf('\t-----\t-------------\n');
0244 fprintf('\t%4d\t%10.2f\n', [ vv(rejected), B(rejected) ]' );
0245 end
0246
0247
0248 k = find( B < bad_threshold );
0249 vv = vv(k);
0250 nm = length(vv);
0251 end
0252
0253 if (baddata == 0)
0254 converged = 1;
0255 if mpopt.verbose
0256 fprintf('\nNo remaining bad data, after discarding data %d time(s).\n', ibd-1);
0257 fprintf('Largest value of B = %.2f\n', maxB);
0258 end
0259 end
0260 ibd = ibd + 1;
0261 end