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