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