diff --git a/case3bus_Ex6_17.m b/case3bus_Ex6_17.m new file mode 100644 index 0000000..4e41c3e --- /dev/null +++ b/case3bus_Ex6_17.m @@ -0,0 +1,57 @@ +function [baseMVA, bus, gen, branch, areas, gencost] = case3bus_Ex6_17 +%CASE3BUS_Ex3_9 Case of 3 bus system. +% From Excecise 3.9 in book 'Computational +% Methods for Electric Power Systems' 2nd Edition by Mariesa Crow page +% 75. +% created by Sami Aldalahmeh on 2/15/2023 +% +% MATPOWER + +%%----- Power Flow Data -----%% +%% system MVA base +baseMVA = 1000; + +%% bus data +% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin +bus = [ + 1 3 0 0 0 0 1 1.02 0 230 1 1.02 0.99; + 2 2 0 0 0 0 1 1.00 0 230 1 1.02 0.99; + 3 1 1200 500 0 0 1 1.00 0 230 1 1.02 0.99; +]; + +%% generator data +% Note: +% 1)It's better of gen to be in number order, otherwise gen and genbid +% should be sorted to make the lp solution output clearly(in number order as well) +% 2)set Pmax to nonzero. set to 999 if no limit +% 3)If change the order of gen, then must change the order in genbid +% accordingly +% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin +gen = [ + 1 0 0 999 -999 1.02 100 1 999 -999; + 2 500 0 999 -999 1.00 100 1 999 -999; +% 3 0 0 999 -999 1.00 baseMVA 1 999 -999; +]; +%gen(:, 9) = 999; % inactive the Pmax constraints + +%% branch data +% fbus tbus r x b rateA rateB rateC ratio angle status +branch = [ + 1 2 0.02 0.3 0.150 0 0 0 0 0 1; + 1 3 0.01 0.1 0.1 0 0 0 0 0 1; + 2 3 0.01 0.1 0.1 0 0 0 0 0 1; +]; + +%%----- OPF Data -----%% +%% area data +areas = [ + 1 1; +]; + +%% generator cost data +% 2 startup shutdown n c(n-1) ... c0 +gencost = [ + 2 0 0 1 1.5 1 0; + 2 0 0 2 1 2 0; +]; +return; diff --git a/cmptSmat.m b/cmptSmat.m new file mode 100644 index 0000000..a8c3ee8 --- /dev/null +++ b/cmptSmat.m @@ -0,0 +1,23 @@ +function [Sfe, Ste] = cmptSmat(V, Ybus, branch) +% Compute the complex power matrices (from-to) "Sfe" and the (to-from) +% "Ste". + + % Number of buses + nb = size(Ybus, 1); + % Overall complex power matrix + Sij = V * V'.*conj(Ybus) - (V .* conj(V) * ones(1,nb)).*conj(Ybus); + + % Find the branch from-to branch number the branch matric + brncInd = branch(:,[1 2]); + + % Convert to linear indicies + sijFromInd = sub2ind( size(Sij), brncInd(:,1), brncInd(:,2) ); + + % Extract from-to elements + Sfe = Sij(sijFromInd); + + % Find linear indicies for the to-from matrix + sijToInd = sub2ind( size(Sij), brncInd(:,2), brncInd(:,1) ); + + % Extract from-to elements + Ste = Sij(sijToInd); \ No newline at end of file diff --git a/doSEerr.m b/doSEerr.m new file mode 100644 index 0000000..c5924a3 --- /dev/null +++ b/doSEerr.m @@ -0,0 +1,294 @@ +function [V, converged, iterNum, z, z_est, error_sqrsum] = doSEext(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma) +%DOSE Do state estimation (extended version). +% created by Rui Bo on 2007/11/12 + +% MATPOWER +% Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) +% by Rui Bo +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MATPOWER/mx-se. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mx-se/ for more info. + +%% Debugging +% load se_org +%% define named indices into bus, gen, branch matrices +[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... + VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; +[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ... + RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch; +[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ... + GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen; + +%% options +tol = 1e-5; % mpopt.pf.tol; +max_it = 100; % mpopt.pf.nr.max_it; +verbose = 0; + +%% initialize +converged = 0; +i = 0; +V = V0; +Va = angle(V); +Vm = abs(V); + +nb = size(Ybus, 1); +f = branch(:, F_BUS); %% list of "from" buses +t = branch(:, T_BUS); %% list of "to" buses + +%% get non reference buses +nonref = [pv;pq]; +% For PV buses the required state is the voltage angle. +% Whereas for PQ buses both the voltage magnitude and angle are required. +nonref_a = [pv;pq]; % angle buses indices +nonref_m = [pq]; % magnitude buses indices + +% !!!!!! +% nonref_a = nonref; +% nonref_m = nonref; + +% Generators buses +gbus = gen(:, GEN_BUS); +gBusInd = 1:length(gbus); + +% Create indencies for state vector. +% The first elements are the voltages' angles then come the magnitudes. +% Note that for the PQ buses, the states are the angles and the magnitudes, +% hence double the number of buses. +Va_ind = 1:(length(pv)+length(pq)); +Vm_ind = length(Va_ind)+1:length(pv)+2*length(pq); + +%% form measurement vector 'z'. NOTE: all are p.u. values +z = [ + measure.PD % **addition to code** + measure.PF + measure.PT + measure.PG + measure.Va + measure.QD % **addition to code** + measure.QF + measure.QT + measure.QG + measure.Vm + ]; + +%% form measurement index vectors +idx_zPD = idx.idx_zPD; % **addition to code** +idx_zQD = idx.idx_zQD; % **addition to code** +idx_zPF = idx.idx_zPF; +idx_zPT = idx.idx_zPT; +idx_zPG = idx.idx_zPG; +idx_zVa = idx.idx_zVa; +idx_zQF = idx.idx_zQF; +idx_zQT = idx.idx_zQT; +idx_zQG = idx.idx_zQG; +idx_zVm = idx.idx_zVm; + +% idx_zPG = gBusInd(ismember(gbus,idx_zPG)); +% idx_zQG = gBusInd(ismember(gbus,idx_zQG)); + +%% get R inverse matrix +% **addition to code** + +if length(sigma.sigma_PF) > 1 + + sigma_vector = [ + sigma.sigma_PD % **addition to code** + sigma.sigma_PF + sigma.sigma_PT + sigma.sigma_PG + sigma.sigma_Va + sigma.sigma_QD % **addition to code** + sigma.sigma_QF + sigma.sigma_QT + sigma.sigma_QG + sigma.sigma_Vm + ]; +else + % % **remove from code** + sigma_vector = [ + sigma.sigma_PD*ones(size(idx_zPD, 1), 1) % **addition to code** + sigma.sigma_PF*ones(size(idx_zPF, 1), 1) + sigma.sigma_PT*ones(size(idx_zPT, 1), 1) + sigma.sigma_PG*ones(size(idx_zPG, 1), 1) + sigma.sigma_Va*ones(size(idx_zVa, 1), 1) + sigma.sigma_QD*ones(size(idx_zQD, 1), 1) % **addition to code** + sigma.sigma_QF*ones(size(idx_zQF, 1), 1) + sigma.sigma_QT*ones(size(idx_zQT, 1), 1) + sigma.sigma_QG*ones(size(idx_zQG, 1), 1) + sigma.sigma_Vm*ones(size(idx_zVm, 1), 1) + ]; % NOTE: zero-valued elements of simga are skipped +end +sigma_square = sigma_vector.^2; +% R_inv = diag(1./sigma_square); +R_inv = inv(diag(1./sigma_square)); + +%% do Newton iterations +while (~converged && i < max_it) + %% update iteration counter + i = i + 1; + + %% --- compute estimated measurement --- + % **remove from code** + Sfe = V(f) .* conj(Yf * V); + Ste = V(t) .* conj(Yt * V); % + + %% Power injection at buses **addition to code** +% [Sfe, Ste] = cmptSmat(V, Ybus, branch); + + + %% Power demand at buses + Sbus = V .* conj(Ybus * V); + + %% compute net injection at generator buses +% gbus = gen(:, GEN_BUS); + Sgbus = V(gbus) .* conj(Ybus(gbus, :) * V); + Sgen = Sgbus * baseMVA + (bus(gbus, PD) + 1j*bus(gbus, QD)); %% inj S + local Sd + Sgen = Sgen/baseMVA; + z_est = [ % NOTE: all are p.u. values + real(Sbus(idx_zPD)); % Demand real power **addition to code** + real(Sfe(idx_zPF)); + real(Ste(idx_zPT)); + real(Sgen(idx_zPG)); + angle(V(idx_zVa)); + imag(Sbus(idx_zQD)); % Demand reactive power **addition to code** + imag(Sfe(idx_zQF)); + imag(Ste(idx_zQT)); + imag(Sgen(idx_zQG)); + abs(V(idx_zVm)); + ]; + + %% --- get H matrix --- + [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V); + [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V); +% genbus_row = findBusRowByIdx(bus, gbus); + genbus_row = gbus; %% rdz, this should be fine if using internal bus numbering + + %% get sub-matrix of H relating to power demand **addition to code** + dPD_dVa = real(dSbus_dVa); + dPD_dVm = real(dSbus_dVm); + + dQD_dVa = imag(dSbus_dVa); + dQD_dVm = imag(dSbus_dVm); + + %% get sub-matrix of H relating to line flow + dPF_dVa = real(dSf_dVa); % from end + dQF_dVa = imag(dSf_dVa); + dPF_dVm = real(dSf_dVm); + dQF_dVm = imag(dSf_dVm); + dPT_dVa = real(dSt_dVa);% to end + dQT_dVa = imag(dSt_dVa); + dPT_dVm = real(dSt_dVm); + dQT_dVm = imag(dSt_dVm); + + %% get sub-matrix of H relating to generator output + dPG_dVa = real(dSbus_dVa(genbus_row, :)); + dQG_dVa = imag(dSbus_dVa(genbus_row, :)); + dPG_dVm = real(dSbus_dVm(genbus_row, :)); + dQG_dVm = imag(dSbus_dVm(genbus_row, :)); + + %% get sub-matrix of H relating to voltage angle + dVa_dVa = eye(nb); + dVa_dVm = zeros(nb, nb); + + %% get sub-matrix of H relating to voltage magnitude + dVm_dVa = zeros(nb, nb); + dVm_dVm = eye(nb); + % **addition to code** + H = [ + dPD_dVa(idx_zPD, nonref_a) dPD_dVm(idx_zPD, nonref_m);% PD deriv. c + dPF_dVa(idx_zPF, nonref_a) dPF_dVm(idx_zPF, nonref_m); + dPT_dVa(idx_zPT, nonref_a) dPT_dVm(idx_zPT, nonref_m); + dPG_dVa(idx_zPG, nonref_a) dPG_dVm(idx_zPG, nonref_m); + dVa_dVa(idx_zVa, nonref_a) dVa_dVm(idx_zVa, nonref_m); + dQD_dVa(idx_zQD, nonref_a) dQD_dVm(idx_zQD, nonref_m);% QD deriv. c + dQF_dVa(idx_zQF, nonref_a) dQF_dVm(idx_zQF, nonref_m); + dQT_dVa(idx_zQT, nonref_a) dQT_dVm(idx_zQT, nonref_m); + dQG_dVa(idx_zQG, nonref_a) dQG_dVm(idx_zQG, nonref_m); + dVm_dVa(idx_zVm, nonref_a) dVm_dVm(idx_zVm, nonref_m); + ]; +% H = [ +% dPD_dVa(idx_zPD, nonref) dPD_dVm(idx_zPD, nonref);% PD deriv. **addition to code** +% dPF_dVa(idx_zPF, nonref) dPF_dVm(idx_zPF, nonref); +% dPT_dVa(idx_zPT, nonref) dPT_dVm(idx_zPT, nonref); +% dPG_dVa(idx_zPG, nonref) dPG_dVm(idx_zPG, nonref); +% dVa_dVa(idx_zVa, nonref) dVa_dVm(idx_zVa, nonref); +% dQD_dVa(idx_zQD, nonref) dQD_dVm(idx_zQD, nonref); +% dQF_dVa(idx_zQF, nonref) dQF_dVm(idx_zQF, nonref); +% dQT_dVa(idx_zQT, nonref) dQT_dVm(idx_zQT, nonref); +% dQG_dVa(idx_zQG, nonref) dQG_dVm(idx_zQG, nonref); +% dVm_dVa(idx_zVm, nonref) dVm_dVm(idx_zVm, nonref); +% ]; + + %% compute update step + J = H'*R_inv*H; % gain matrix + F = H'*R_inv*(z-z_est); % evalute F(x) + if ~isobservable(H, pv, pq) + error('doSE: system is not observable'); + end + dx = (J \ F); + %% Debugging + DEBUG_FLAG = true; + + if DEBUG_FLAG + % Rearange H matrix as ex. 6.17 in book + Hd = H; + Hd([1 end],:)=Hd([end 1],:); + Hd([end end-1],:)=Hd([end-1 end],:); + fprintf('H matrix in interation #%d\n',i) + fprintf('H = \n') + disp(Hd) + + end + + + %% check for convergence + normF = norm(F, inf); + if verbose > 1 + fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF); + end + if normF < tol + converged = 1; + end + + %% update voltage +% Va(nonref) = Va(nonref) + dx(1:size(nonref, 1)); + Va(nonref_a) = Va(nonref_a) + dx(Va_ind); % **addition code** +% Vm(nonref) = Vm(nonref) + dx(size(nonref, 1)+1:2*size(nonref, 1)); + Vm(nonref_m) = Vm(nonref_m) + dx(Vm_ind);% **addition code** + +% y = [Va(nonref_a);Vm(nonref_m)]; % **addition code** + + V = Vm .* exp(1j * Va); % NOTE: angle is in radians in pf solver, but in degree in case data + Vm = abs(V); %% update Vm and Va again in case + Va = angle(V); %% we wrapped around with a negative Vm + +% %% debugging +% re = sum([abs(abs(V_cell{i}) - Vm),abs(angle(V_cell{i}) - Va)]); +% disp(['iter: ', num2str(i), ' Error: ' num2str(re) ]) +% disp('') +end + +iterNum = i; + +%% get weighted sum of squared errors +error_sqrsum = sum((z - z_est).^2./sigma_square); + +%% New function added to code +function BusInd = getInd(ind, pv, pq) +% Return the bus index from the H matrix column index. + +Lpv = length(pv); +Lpq = length(pq); + +for i = 1:length(ind) + if ind(i) <= Lpv + BusInd(i) = pv(ind(i)); + elseif (ind(i) > Lpv) && (ind(i) <= Lpv+Lpq) + BusInd(i) = pq(ind(i) - Lpv); + else + BusInd(i) = pq(ind(i) - Lpv - Lpq); + end +end + diff --git a/doSEmod.m b/doSEmod.m new file mode 100644 index 0000000..7f21fa7 --- /dev/null +++ b/doSEmod.m @@ -0,0 +1,276 @@ +function [V, converged, iterNum, z, z_est, error_sqrsum] = doSEext(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma) +%DOSE Do state estimation (extended version). +% created by Rui Bo on 2007/11/12 + +% MATPOWER +% Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) +% by Rui Bo +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MATPOWER/mx-se. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mx-se/ for more info. + +%% Debugging +% load se_org +%% define named indices into bus, gen, branch matrices +[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... + VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; +[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ... + RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch; +[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ... + GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen; + +%% options +tol = 1e-5; % mpopt.pf.tol; +max_it = 100; % mpopt.pf.nr.max_it; +verbose = 0; + +%% initialize +converged = 0; +i = 0; +V = V0; +Va = angle(V); +Vm = abs(V); + +nb = size(Ybus, 1); +f = branch(:, F_BUS); %% list of "from" buses +t = branch(:, T_BUS); %% list of "to" buses + +%% get non reference buses +nonref = [pv;pq]; +% For PV buses the required state is the voltage angle. +% Whereas for PQ buses both the voltage magnitude and angle are required. +nonref_a = [pv;pq]; % angle buses indices +nonref_m = [pq]; % magnitude buses indices + +% !!!!!! +% nonref_a = nonref; +% nonref_m = nonref; + +% Generators buses +gbus = gen(:, GEN_BUS); +gBusInd = 1:length(gbus); + +% Create indencies for state vector. +% The first elements are the voltages' angles then come the magnitudes. +% Note that for the PQ buses, the states are the angles and the magnitudes, +% hence double the number of buses. +Va_ind = 1:(length(pv)+length(pq)); +Vm_ind = length(Va_ind)+1:length(pv)+2*length(pq); + +%% form measurement vector 'z'. NOTE: all are p.u. values +z = [ + measure.PD % **addition to code** + measure.PF + measure.PT + measure.PG + measure.Va + measure.QD % **addition to code** + measure.QF + measure.QT + measure.QG + measure.Vm + ]; + +%% form measurement index vectors +idx_zPD = idx.idx_zPD; % **addition to code** +idx_zQD = idx.idx_zQD; % **addition to code** +idx_zPF = idx.idx_zPF; +idx_zPT = idx.idx_zPT; +idx_zPG = idx.idx_zPG; +idx_zVa = idx.idx_zVa; +idx_zQF = idx.idx_zQF; +idx_zQT = idx.idx_zQT; +idx_zQG = idx.idx_zQG; +idx_zVm = idx.idx_zVm; + +% idx_zPG = gBusInd(ismember(gbus,idx_zPG)); +% idx_zQG = gBusInd(ismember(gbus,idx_zQG)); + +%% get R inverse matrix +% **addition to code** + +if length(sigma.sigma_PF) > 1 + + sigma_vector = [ + sigma.sigma_PD % **addition to code** + sigma.sigma_PF + sigma.sigma_PT + sigma.sigma_PG + sigma.sigma_Va + sigma.sigma_QD % **addition to code** + sigma.sigma_QF + sigma.sigma_QT + sigma.sigma_QG + sigma.sigma_Vm + ]; +else + % % **remove from code** + sigma_vector = [ + sigma.sigma_PD*ones(size(idx_zPD, 1), 1) % **addition to code** + sigma.sigma_PF*ones(size(idx_zPF, 1), 1) + sigma.sigma_PT*ones(size(idx_zPT, 1), 1) + sigma.sigma_PG*ones(size(idx_zPG, 1), 1) + sigma.sigma_Va*ones(size(idx_zVa, 1), 1) + sigma.sigma_QD*ones(size(idx_zQD, 1), 1) % **addition to code** + sigma.sigma_QF*ones(size(idx_zQF, 1), 1) + sigma.sigma_QT*ones(size(idx_zQT, 1), 1) + sigma.sigma_QG*ones(size(idx_zQG, 1), 1) + sigma.sigma_Vm*ones(size(idx_zVm, 1), 1) + ]; % NOTE: zero-valued elements of simga are skipped +end +sigma_square = sigma_vector.^2; +% R_inv = diag(1./sigma_square); +R_inv = inv(diag(1./sigma_square)); + +%% do Newton iterations +while (~converged && i < max_it) + %% update iteration counter + i = i + 1; + + %% --- compute estimated measurement --- + % **remove from code** +% Sfe = V(f) .* conj(Yf * V); +% Ste = V(t) .* conj(Yt * V); % + + %% Power injection at buses **addition to code** + [Sfe, Ste] = cmptSmat(V, Ybus, branch); + + + %% Power demand at buses + Sbus = V .* conj(Ybus * V); + + %% compute net injection at generator buses +% gbus = gen(:, GEN_BUS); + Sgbus = V(gbus) .* conj(Ybus(gbus, :) * V); + Sgen = Sgbus * baseMVA + (bus(gbus, PD) + 1j*bus(gbus, QD)); %% inj S + local Sd + Sgen = Sgen/baseMVA; + z_est = [ % NOTE: all are p.u. values + real(Sbus(idx_zPD)); % Demand real power **addition to code** + real(Sfe(idx_zPF)); + real(Ste(idx_zPT)); + real(Sgen(idx_zPG)); + angle(V(idx_zVa)); + imag(Sbus(idx_zQD)); % Demand reactive power **addition to code** + imag(Sfe(idx_zQF)); + imag(Ste(idx_zQT)); + imag(Sgen(idx_zQG)); + abs(V(idx_zVm)); + ]; + + %% --- get H matrix --- + [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V); + [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V); +% genbus_row = findBusRowByIdx(bus, gbus); + genbus_row = gbus; %% rdz, this should be fine if using internal bus numbering + + %% get sub-matrix of H relating to power demand **addition to code** + dPD_dVa = real(dSbus_dVa); + dPD_dVm = real(dSbus_dVm); + + dQD_dVa = imag(dSbus_dVa); + dQD_dVm = imag(dSbus_dVm); + + %% get sub-matrix of H relating to line flow + dPF_dVa = real(dSf_dVa); % from end + dQF_dVa = imag(dSf_dVa); + dPF_dVm = real(dSf_dVm); + dQF_dVm = imag(dSf_dVm); + dPT_dVa = real(dSt_dVa);% to end + dQT_dVa = imag(dSt_dVa); + dPT_dVm = real(dSt_dVm); + dQT_dVm = imag(dSt_dVm); + + %% get sub-matrix of H relating to generator output + dPG_dVa = real(dSbus_dVa(genbus_row, :)); + dQG_dVa = imag(dSbus_dVa(genbus_row, :)); + dPG_dVm = real(dSbus_dVm(genbus_row, :)); + dQG_dVm = imag(dSbus_dVm(genbus_row, :)); + + %% get sub-matrix of H relating to voltage angle + dVa_dVa = eye(nb); + dVa_dVm = zeros(nb, nb); + + %% get sub-matrix of H relating to voltage magnitude + dVm_dVa = zeros(nb, nb); + dVm_dVm = eye(nb); + % **addition to code** + H = [ + dPD_dVa(idx_zPD, nonref_a) dPD_dVm(idx_zPD, nonref_m);% PD deriv. c + dPF_dVa(idx_zPF, nonref_a) dPF_dVm(idx_zPF, nonref_m); + dPT_dVa(idx_zPT, nonref_a) dPT_dVm(idx_zPT, nonref_m); + dPG_dVa(idx_zPG, nonref_a) dPG_dVm(idx_zPG, nonref_m); + dVa_dVa(idx_zVa, nonref_a) dVa_dVm(idx_zVa, nonref_m); + dQD_dVa(idx_zQD, nonref_a) dQD_dVm(idx_zQD, nonref_m);% QD deriv. c + dQF_dVa(idx_zQF, nonref_a) dQF_dVm(idx_zQF, nonref_m); + dQT_dVa(idx_zQT, nonref_a) dQT_dVm(idx_zQT, nonref_m); + dQG_dVa(idx_zQG, nonref_a) dQG_dVm(idx_zQG, nonref_m); + dVm_dVa(idx_zVm, nonref_a) dVm_dVm(idx_zVm, nonref_m); + ]; +% H = [ +% dPD_dVa(idx_zPD, nonref) dPD_dVm(idx_zPD, nonref);% PD deriv. **addition to code** +% dPF_dVa(idx_zPF, nonref) dPF_dVm(idx_zPF, nonref); +% dPT_dVa(idx_zPT, nonref) dPT_dVm(idx_zPT, nonref); +% dPG_dVa(idx_zPG, nonref) dPG_dVm(idx_zPG, nonref); +% dVa_dVa(idx_zVa, nonref) dVa_dVm(idx_zVa, nonref); +% dQD_dVa(idx_zQD, nonref) dQD_dVm(idx_zQD, nonref); +% dQF_dVa(idx_zQF, nonref) dQF_dVm(idx_zQF, nonref); +% dQT_dVa(idx_zQT, nonref) dQT_dVm(idx_zQT, nonref); +% dQG_dVa(idx_zQG, nonref) dQG_dVm(idx_zQG, nonref); +% dVm_dVa(idx_zVm, nonref) dVm_dVm(idx_zVm, nonref); +% ]; + + %% compute update step + J = H'*R_inv*H; % gain matrix + F = H'*R_inv*(z-z_est); % evalute F(x) + if ~isobservable(H, pv, pq) + error('doSE: system is not observable'); + end + dx = (J \ F); + %% Debugging + DEBUG_FLAG = true; + + if DEBUG_FLAG + % Rearange H matrix as ex. 6.17 in book + Hd = H; + Hd([1 end],:)=Hd([end 1],:); + Hd([end end-1],:)=Hd([end-1 end],:); + fprintf('H matrix in interation #%d\n',i) + fprintf('H = \n') + disp(Hd) + + end + + + %% check for convergence + normF = norm(F, inf); + if verbose > 1 + fprintf('\niteration [%3d]\t\tnorm of mismatch: %10.3e', i, normF); + end + if normF < tol + converged = 1; + end + + %% update voltage +% Va(nonref) = Va(nonref) + dx(1:size(nonref, 1)); + Va(nonref_a) = Va(nonref_a) + dx(Va_ind); % **addition code** +% Vm(nonref) = Vm(nonref) + dx(size(nonref, 1)+1:2*size(nonref, 1)); + Vm(nonref_m) = Vm(nonref_m) + dx(Vm_ind);% **addition code** + +% y = [Va(nonref_a);Vm(nonref_m)]; % **addition code** + + V = Vm .* exp(1j * Va); % NOTE: angle is in radians in pf solver, but in degree in case data + Vm = abs(V); %% update Vm and Va again in case + Va = angle(V); %% we wrapped around with a negative Vm + +% %% debugging +% re = sum([abs(abs(V_cell{i}) - Vm),abs(angle(V_cell{i}) - Va)]); +% disp(['iter: ', num2str(i), ' Error: ' num2str(re) ]) +% disp('') +end + +iterNum = i; + +%% get weighted sum of squared errors +error_sqrsum = sum((z - z_est).^2./sigma_square); diff --git a/run_se_err.m b/run_se_err.m new file mode 100644 index 0000000..b8c4ef0 --- /dev/null +++ b/run_se_err.m @@ -0,0 +1,85 @@ +function [baseMVA, bus, gen, branch, success, et, z, z_est, error_sqrsum] = ... + run_se_ext(casename, measure, idx, sigma, type_initialguess, V0) +%RUN_SE Run state estimation. +% [INPUT PARAMETERS] +% measure: measurements +% idx: measurement indices +% sigma: measurement variances +% [OUTPUT PARAMETERS] +% z: Measurement Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm (if +% applicable), so it has ordered differently from original measurements +% z_est: Estimated Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm +% (if applicable) +% error_sqrsum: Weighted sum of error squares +% created by Rui Bo on 2007/11/12 + +% MATPOWER +% Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) +% by Rui Bo +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MATPOWER/mx-se. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mx-se/ for more info. + +%% define named indices into data matrices +[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... + MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... + QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; + +%% read data & convert to internal bus numbering +[baseMVA, bus, gen, branch] = loadcase(casename); +[i2e, bus, gen, branch] = ext2int(bus, gen, branch); + +%% get bus index lists of each type of bus +[ref, pv, pq] = bustypes(bus, gen); + +%% build admittance matrices +[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); +Ybus = full(Ybus); +Yf = full(Yf); +Yt = full(Yt); + +%% prepare initial guess +if nargin < 6 + V0 = getV0(bus, gen, type_initialguess); +else + V0 = getV0(bus, gen, type_initialguess, V0); +end + +%% run state estimation +t0 = tic; +[V, success, iterNum, z, z_est, error_sqrsum] = ... + doSEerr(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma); + +% update Pg and Qg using estimated values idx.idx_zPD +if length(idx.idx_zPG) > 0 + % ** addition to code ** +% cnt = length(idx.idx_zPD) + length(idx.idx_zQD)+ length(idx.idx_zPF)... +% + length(idx.idx_zPT); + cnt = length(idx.idx_zPF) + length(idx.idx_zPT); + gen(idx.idx_zPG, PG) = z_est([cnt+1:cnt+length(idx.idx_zPG)]) * baseMVA; +end +if length(idx.idx_zQG) > 0 + % ** addition to code ** +% cnt = length(idx.idx_zPD) + length(idx.idx_zQD)+ length(idx.idx_zPF)... +% + length(idx.idx_zPT) + length(idx.idx_zPG) + length(idx.idx_zVa)... +% + length(idx.idx_zQF) + length(idx.idx_zQT); + cnt = length(idx.idx_zPF) + length(idx.idx_zPT) + length(idx.idx_zPG) + ... + length(idx.idx_zVa) + length(idx.idx_zQF) + length(idx.idx_zQT); + gen(idx.idx_zQG, QG) = z_est([cnt+1:cnt+length(idx.idx_zQG)]) * baseMVA; +end + +%% update data matrices with solution, ie, V +% [bus, gen, branch] = updatepfsoln(baseMVA, bus, gen, branch, Ybus, V, ref, pv, pq); +[bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq); +et = toc(t0); + +%%----- output results ----- +%% convert back to original bus numbering & print results +[bus, gen, branch] = int2ext(i2e, bus, gen, branch); +%% output power flow solution +outputpfsoln(baseMVA, bus, gen, branch, success, et, 1, iterNum); +%% output state estimation solution +outputsesoln(idx, sigma, z, z_est, error_sqrsum); + diff --git a/run_se_mod.m b/run_se_mod.m new file mode 100644 index 0000000..08011de --- /dev/null +++ b/run_se_mod.m @@ -0,0 +1,85 @@ +function [baseMVA, bus, gen, branch, success, et, z, z_est, error_sqrsum] = ... + run_se_ext(casename, measure, idx, sigma, type_initialguess, V0) +%RUN_SE Run state estimation. +% [INPUT PARAMETERS] +% measure: measurements +% idx: measurement indices +% sigma: measurement variances +% [OUTPUT PARAMETERS] +% z: Measurement Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm (if +% applicable), so it has ordered differently from original measurements +% z_est: Estimated Vector. In the order of PF, PT, PG, Va, QF, QT, QG, Vm +% (if applicable) +% error_sqrsum: Weighted sum of error squares +% created by Rui Bo on 2007/11/12 + +% MATPOWER +% Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) +% by Rui Bo +% and Ray Zimmerman, PSERC Cornell +% +% This file is part of MATPOWER/mx-se. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mx-se/ for more info. + +%% define named indices into data matrices +[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... + MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... + QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; + +%% read data & convert to internal bus numbering +[baseMVA, bus, gen, branch] = loadcase(casename); +[i2e, bus, gen, branch] = ext2int(bus, gen, branch); + +%% get bus index lists of each type of bus +[ref, pv, pq] = bustypes(bus, gen); + +%% build admittance matrices +[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); +Ybus = full(Ybus); +Yf = full(Yf); +Yt = full(Yt); + +%% prepare initial guess +if nargin < 6 + V0 = getV0(bus, gen, type_initialguess); +else + V0 = getV0(bus, gen, type_initialguess, V0); +end + +%% run state estimation +t0 = tic; +[V, success, iterNum, z, z_est, error_sqrsum] = ... + doSEmod(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V0, ref, pv, pq, measure, idx, sigma); + +% update Pg and Qg using estimated values idx.idx_zPD +if length(idx.idx_zPG) > 0 + % ** addition to code ** +% cnt = length(idx.idx_zPD) + length(idx.idx_zQD)+ length(idx.idx_zPF)... +% + length(idx.idx_zPT); + cnt = length(idx.idx_zPF) + length(idx.idx_zPT); + gen(idx.idx_zPG, PG) = z_est([cnt+1:cnt+length(idx.idx_zPG)]) * baseMVA; +end +if length(idx.idx_zQG) > 0 + % ** addition to code ** +% cnt = length(idx.idx_zPD) + length(idx.idx_zQD)+ length(idx.idx_zPF)... +% + length(idx.idx_zPT) + length(idx.idx_zPG) + length(idx.idx_zVa)... +% + length(idx.idx_zQF) + length(idx.idx_zQT); + cnt = length(idx.idx_zPF) + length(idx.idx_zPT) + length(idx.idx_zPG) + ... + length(idx.idx_zVa) + length(idx.idx_zQF) + length(idx.idx_zQT); + gen(idx.idx_zQG, QG) = z_est([cnt+1:cnt+length(idx.idx_zQG)]) * baseMVA; +end + +%% update data matrices with solution, ie, V +% [bus, gen, branch] = updatepfsoln(baseMVA, bus, gen, branch, Ybus, V, ref, pv, pq); +[bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq); +et = toc(t0); + +%%----- output results ----- +%% convert back to original bus numbering & print results +[bus, gen, branch] = int2ext(i2e, bus, gen, branch); +%% output power flow solution +outputpfsoln(baseMVA, bus, gen, branch, success, et, 1, iterNum); +%% output state estimation solution +outputsesoln(idx, sigma, z, z_est, error_sqrsum); + diff --git a/test_se_Ex_6_17.m b/test_se_Ex_6_17.m new file mode 100644 index 0000000..c55cf9f --- /dev/null +++ b/test_se_Ex_6_17.m @@ -0,0 +1,71 @@ +function test_se_Ex_6_17 +%TEST_SE Test state estimation. +% created by Rui Bo on 2007/11/12 + +% MATPOWER +% Copyright (c) 2009-2016, Power Systems Engineering Research Center (PSERC) +% by Rui Bo +% +% This file is part of MATPOWER/mx-se. +% Covered by the 3-clause BSD License (see LICENSE file for details). +% See https://github.com/MATPOWER/mx-se/ for more info. + +%%------------------------------------------------------ +% using data in Problem 6.17 in book 'Computational +% Methods for Electric Power Systems' by Mariesa Crow +%%------------------------------------------------------ +%% which measurements are available +idx.idx_zPD = [3]; % P3 -> Power demand at bus 3**addition to code** +idx.idx_zQD = []; % Q3 -> Reactive Power demand **addition to code** +idx.idx_zPF = [2]; % P13 -> 2nd entry in the branch matrix **addition to code** +idx.idx_zPT = []; +idx.idx_zPG = []; +idx.idx_zVa = []; +idx.idx_zQF = []; +idx.idx_zQT = [1]; % Q21 -> 1st entry in the branch matrix +idx.idx_zQG = [2]; % Q2 -> Injected reactive power at bus 2 +idx.idx_zVm = [3]; % V3 -> Voltage magnitude at bus 3 + +%% specify measurements +measure.PD = [-1.181]; % **addition to code** +measure.QD = []; % **addition to code** +measure.PF = [0.668]; +measure.PT = []; +measure.PG = []; +measure.Va = []; +measure.QF = []; +measure.QT = [-0.082]; +measure.QG = [-0.086]; +measure.Vm = [0.975]; + +%% specify measurement variances +sigma.sigma_PD = [0.050] ; % **addition to code** +sigma.sigma_QD = [] ; % **addition to code** +sigma.sigma_PF = [0.050] ; +sigma.sigma_PT = []; +sigma.sigma_PG = []; +sigma.sigma_Va = []; +sigma.sigma_QF = []; +sigma.sigma_QT = [0.075]; +sigma.sigma_QG = [0.075]; +sigma.sigma_Vm = 0.01; + +%% check input data integrity +nbus = 3; +[success, measure, idx, sigma] = checkDataIntegrity(measure, idx, sigma, nbus); +% if ~success +% error('State Estimation input data are not complete or sufficient!'); +% end + +%% run state estimation +casename = 'case3bus_Ex6_17.m'; +type_initialguess = 2; % flat start + +fprintf('State estimation using the original doSE implemetination\n') +[baseMVA, bus, gen, branch, success, et, z, z_est, error_sqrsum] ... + = run_se_err(casename, measure, idx, sigma, type_initialguess); + +fprintf('\nState estimation using the modified doSE implemetination\n') +[baseMVA, bus, gen, branch, success, et, z, z_est, error_sqrsum] ... + = run_se_mod(casename, measure, idx, sigma, type_initialguess); +% = run_se_ext(casename, measure, idx, sigma, type_initialguess);