From c6f4e608cfafb7dd619f352716d6c25d8a5c9190 Mon Sep 17 00:00:00 2001 From: aldalahmeh <109239143+aldalahmeh@users.noreply.github.com> Date: Fri, 29 Sep 2023 11:39:36 +0300 Subject: [PATCH] Fix problem in state estimation It seems that there is a bug in computing the S matrix that affects the compulation of the derivative matrix H. As, such a new function has been developed to compute S. Moreover, the code is hardwired to solve the problem 6.7 in the book "Computational Methods for Electric Power Systems" by Mariesa Crow, 1st edition. In particular, it only accepts the measurements in that problem. So, the code is expanded to accommodate any sets of measurements. A test file "test_se_Ex_6_17.m" is provided that solves example 6.17 in the above 2nd edit of the book. The file shows the difference in computing H between the original, which is edited to allow more measurements, and the modified code. The modified code produced results that coincide with the solution of example 6.17 in the book. --- se/case3bus_Ex6_17.m | 57 +++++++++ se/cmptSmat.m | 23 ++++ se/doSEerr.m | 294 +++++++++++++++++++++++++++++++++++++++++++ se/doSEmod.m | 276 ++++++++++++++++++++++++++++++++++++++++ se/run_se_err.m | 85 +++++++++++++ se/run_se_mod.m | 85 +++++++++++++ se/test_se_Ex_6_17.m | 71 +++++++++++ 7 files changed, 891 insertions(+) create mode 100644 se/case3bus_Ex6_17.m create mode 100644 se/cmptSmat.m create mode 100644 se/doSEerr.m create mode 100644 se/doSEmod.m create mode 100644 se/run_se_err.m create mode 100644 se/run_se_mod.m create mode 100644 se/test_se_Ex_6_17.m diff --git a/se/case3bus_Ex6_17.m b/se/case3bus_Ex6_17.m new file mode 100644 index 0000000..4e41c3e --- /dev/null +++ b/se/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/se/cmptSmat.m b/se/cmptSmat.m new file mode 100644 index 0000000..a8c3ee8 --- /dev/null +++ b/se/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/se/doSEerr.m b/se/doSEerr.m new file mode 100644 index 0000000..c5924a3 --- /dev/null +++ b/se/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/se/doSEmod.m b/se/doSEmod.m new file mode 100644 index 0000000..7f21fa7 --- /dev/null +++ b/se/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/se/run_se_err.m b/se/run_se_err.m new file mode 100644 index 0000000..b8c4ef0 --- /dev/null +++ b/se/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/se/run_se_mod.m b/se/run_se_mod.m new file mode 100644 index 0000000..08011de --- /dev/null +++ b/se/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/se/test_se_Ex_6_17.m b/se/test_se_Ex_6_17.m new file mode 100644 index 0000000..c55cf9f --- /dev/null +++ b/se/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);