Skip to content

Transformer taps and phases #16

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
GorazdB wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Transformer taps and phases #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7d622b6
Tap and Phases calculation for transformers.
GorazdB Jun 14, 2017
f75d11b
Tap and Phases calculation for transformers.
GorazdB Jun 14, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 23 additions & 0 deletions lib/example_tap_phase_calc.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
%An example of running the tap changer and phase shifter analysis on the
%IEEE 300 bus system with two tap changers and two phase shifters


%% Define case data
file = 'case300';

%% 1. Define the transformers
%tap_changers_data = [line of insertion, control node, control voltage]
tap_changers_data = [
1,9001,.95
2,9005, 1
];

%phase_shifters_data = [line of insertion, control power]
phase_shifters_data = [
100 , 0
110 , 1
];

%% Run the analysis
[results, stevilo_iteracij, success] = taps_and_phases_analysis(file, tap_changers_data, phase_shifters_data);

155 changes: 155 additions & 0 deletions lib/sensit_to_taps_and_phases.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
function [ sensitivity_of_H_to_U_taps_and_phasess ] = sensit_to_taps_and_phases( results, tap_changers, phase_shifters, Jacobian , tip_vozlisc)
%sensit_to_taps_and_phases Calculates the first orded sensitivities of
%regulated active powers and nodal voltages to taps and phases.
% [sensitivity_of_H_to_U_taps_and_phasess] = sensit_to_taps_and_phases(results, tap_changers, phase_shifters, Jacobian , node_types)
%
% Inputs :
% results: results from a runpf(), interanly idnexed (use ext2int())
% tap_changers: tap changers data as in taps_and_phases_analysis()
% phase_shifters: phase shifter data as in taps_and_phases_analysis()
% Jacobian: power missmatch Jacobian matrix with the appropriate
% ordering, rows go P1, Q1, P2, Q2...., columns go V1, delta1, V2,
% delta2.....
% tip_vozlisc: node types, 2 for PQ, 1 for PV, same as nubmer of
% equations for each node, used for indexing
%
% Output :
% sensitivity_of_H_to_U_taps_and_phasess : first order sensitivities
% of p.u. changes in nodal votlages and line active powers to changes
% in calculated taps and phases of transformers.
%
% by Gorazd Bone, Faculty of Electrical Engineering, Ljubljana

ind_node_eq = tip_vozlisc * 0;
for k=1:size(tip_vozlisc,1)
ind_node_eq(k) = sum(tip_vozlisc(1:k-1)) + 1;
end

stevilo_tap_ch = size(tap_changers,1);
stevilo_pha_sh = size(phase_shifters,1);

%% 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; %#ok<*NASGU,*ASGLU>
[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, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% Define Partial derivative matrices
% G is the LF problem power injection missmatch form (sum of nodal powers)
GpoU = zeros(sum(tip_vozlisc) , stevilo_tap_ch + stevilo_pha_sh); % derivative of LF problem w.r.t. taps & phases
GpoX = Jacobian; % derivative of LF problem w.r.t. voltages and angles
HpoU = zeros(stevilo_tap_ch + stevilo_pha_sh); %derivative of criterion function w.r.t. taps & phases
HpoX = zeros(stevilo_tap_ch + stevilo_pha_sh , sum(tip_vozlisc));%derivative of criterion function w.r.t. voltages and angles

for k=1:stevilo_tap_ch % partial derivatives for taps
veja = tap_changers(k,1);%tap branch
reg_voz = tap_changers(k,2); %ragulated node
voz1 = results.branch(veja,F_BUS);
voz2 = results.branch(veja,T_BUS);
indeks_P1 = ind_node_eq(voz1); %active power and volt. amp index for f node
indeks_P2 = ind_node_eq(voz2); %active power and volt. amp index for t node
indeks_reg_voz = ind_node_eq(reg_voz);%volt. amp index for reg. node

Uvozi = results.bus(voz1,VM);
Uvozj = results.bus(voz2,VM);
di = results.bus(voz1,VA)/180*pi;
dj = results.bus(voz2,VA)/180*pi;
Z = results.branch(veja,BR_R) + 1i*results.branch(veja,BR_X);
Bsh = results.branch(veja,BR_B);
G = real(1/Z);
B = imag(1/Z);
tap = results.branch(veja,TAP);
ph = results.branch(veja,SHIFT)/180*pi;

if tip_vozlisc(reg_voz) == 1 % tap changer regulating PV node
error_string = strcat(2 , 'Tap changer' , int2str(k) , 'regulating the voltage of a PV node');
errordlg(error_string );
return
end

dPijdtap = (Uvozi*(-2*G*Uvozi+G*tap*Uvozj*cos(di-dj-ph)+B*tap*Uvozj*sin(di-dj-ph)))/tap^3;
dPjidtap = (Uvozi*Uvozj*(G*cos(di-dj-ph)-B*sin(di-dj-ph)))/tap^2;
dQijdtap = (Uvozi*((2*B+Bsh)*Uvozi-B*tap*Uvozj*cos(di-dj-ph)+G*tap*Uvozj*sin(di-dj-ph)))/tap^3;
dQjidtap = -((Uvozi*Uvozj*(B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap^2);

%% dG/dU
GpoU(indeks_P1 , k) = - dPijdtap; %dGfrom / dtap
if tip_vozlisc(voz1) == 2 %PQ type of node
GpoU(indeks_P1 + 1 , k) = - dQijdtap;
end
GpoU(indeks_P2 , k) = - dPjidtap;%dGto / dtap
if tip_vozlisc(voz2) == 2 %PQ type of node
GpoU(indeks_P2 + 1 , k) = - dQjidtap;
end

%% dH/dU
HpoU(k,k) = 0;%odvod vozliscne napetosti po tap-u je 0
if stevilo_pha_sh
if find(phase_shifters(:,1) == veja)
HpoU(find(phase_shifters(:,1) == veja) + stevilo_tap_ch , k) = dPijdtap;%dLinePower / dtap - in case a phase shifter and tap changer are in the same line
end
end

%% dH/dX
HpoX(k,indeks_reg_voz) = 1; %dRegNodeVolt / dVoltAmp = unity
end

for k=1:stevilo_pha_sh % partial derivatives for phases
veja = phase_shifters(k,1);%tap branch
voz1 = results.branch(veja,F_BUS);
voz2 = results.branch(veja,T_BUS);
indeks_P1 = ind_node_eq(voz1); %active power and volt. amp index for f node
indeks_P2 = ind_node_eq(voz2); %active power and volt. amp index for t node

Uvozi = results.bus(voz1,VM);
Uvozj = results.bus(voz2,VM);
di = results.bus(voz1,VA)/180*pi;
dj = results.bus(voz2,VA)/180*pi;
Z = results.branch(veja,BR_R) + 1i*results.branch(veja,BR_X);
Bsh = results.branch(veja,BR_B);
G = real(1/Z);
B = imag(1/Z);
tap = results.branch(veja,TAP);
ph = results.branch(veja,SHIFT)/180*pi;

dPijdph = (Uvozi*Uvozj*(B*cos(di-dj-ph)-G*sin(di-dj-ph)))/tap;
dPjidph = -((Uvozi*Uvozj*(B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap);
dPijUi = (2*G*Uvozi-tap*Uvozj*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap^2;
dPijUj = -((Uvozi*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap);
dPijdi = (Uvozi*Uvozj*(-B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap;
dPijdj = (Uvozi*Uvozj*(B*cos(di-dj-ph)-G*sin(di-dj-ph)))/tap;
dQijdph = (Uvozi*Uvozj*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap;
dQjidph = (Uvozi*Uvozj*(-G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap;

%% dG/dU
GpoU(indeks_P1, k + stevilo_tap_ch) = - dPijdph;%dGfrom / dphase
if tip_vozlisc(voz1) == 2 %PQ type of node
GpoU(indeks_P1 + 1, k + stevilo_tap_ch) = - dQijdph;%dGfrom / dphase
end
GpoU(indeks_P2, k + stevilo_tap_ch) = - dPjidph;%dGto / dphase
if tip_vozlisc(voz2) == 2 %PQ type of node
GpoU(indeks_P2 + 1, k + stevilo_tap_ch) = - dQjidph;%dGto / dphase
end

%% dH/dU
HpoU(k + stevilo_tap_ch, k + stevilo_tap_ch) = dPijdph;%dLinePower / dphase

%% dH/dX
if tip_vozlisc(voz1) == 2 %from bus is PQ type
HpoX(k + stevilo_tap_ch, indeks_P1) = dPijUi;%dLinePower / dVoltAmp
HpoX(k + stevilo_tap_ch, indeks_P1+1) = dPijdi;%dLinePower / dVoltAng
else %from bus is PV type
HpoX(k + stevilo_tap_ch, indeks_P1) = dPijdi;%dLinePower / dVoltAng
end

if tip_vozlisc(voz2) == 2 %same thing as abocve for to bus
HpoX(k + stevilo_tap_ch, indeks_P2) = dPijUj;
HpoX(k + stevilo_tap_ch, indeks_P2+1) = dPijdj;
else
HpoX(k + stevilo_tap_ch, indeks_P2) = dPijdj;
end

end

sensitivity_of_H_to_U_taps_and_phasess = HpoU - HpoX * ( (GpoX) \ (GpoU) ); %total sensitivity of line active powers and nodal voltages w.r.t. all taps and phases
68 changes: 68 additions & 0 deletions lib/shiftJac_taps_phases.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
function [ J ] = shiftJac_taps_phases( results )
%shiftJac_taps_phases Constructs the Jacobian matrix using Matpower
%supplied makeJac and reorderes it for taps and phases analysis.
% [J ] = shiftJac_taps_phases( results )
% Inputs :
% results: results from a runpf(), interanly idnexed (use ext2int)
%
% Output :
% J : Jacobian matrix for power mismatch load flow, ordereing
% appropriate fo sensit_to_taps_and_phases()
%
% by Gorazd Bone, Faculty of Electrical Engineering, Ljubljana

%% 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, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; %#ok<*ASGLU,*NASGU>
[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;

pq = results.bus(:,BUS_TYPE)==PQ;
pv = results.bus(:,BUS_TYPE)==PV;
ref = find(pq + pv == 0);

tip_vozlisc = pq * 2 + pv;
ind_node_eq = tip_vozlisc * 0;
for k=1:size(tip_vozlisc,1)
ind_node_eq(k) = sum(tip_vozlisc(1:k-1)) + 1;
end
loc_P_pv = ind_node_eq(pv==1);
loc_P_pq = ind_node_eq(pq==1);
loc_Q_pq = loc_P_pq+1;


per_rows = sparse([loc_P_pv;loc_P_pq;loc_Q_pq] , 1:sum(pq * 2 + pv),1);
per_cols = sparse(1:sum(pq * 2 + pv) , [loc_P_pv;loc_Q_pq;loc_P_pq] , 1);
%
% MatPower Jacobi has the following oredering for the missmatch:
%
% P_vect_of_pvnodes
% P_vect_of_pqnodes
% Q_vect_of_pqnodes
%
% and for the unknowns:
%
% angle_vect_of_pvnodes
% angle_vect_of_pqnodes
% voltage_vect_of_pqnodes
%
% the used Jacobi has the ordering node-by-node, for the missmatch
%
% P_of_node1
% Q_of_node1 (if PQ)
% P_of_node2
% Q_of_node2 (if PQ)
%
% and for the unknowns:
%
% voltage_of_node1 (if PQ)
% anle_of_node1
% voltage_of_node2 (if PQ)
% anle_of_node2 (if PQ)
%
J = - per_rows * makeJac(results.baseMVA, results.bus, results.branch, results.gen) * per_cols;

110 changes: 110 additions & 0 deletions lib/taps_and_phases_analysis.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
function [results, stevilo_iteracij, success] = taps_and_phases_analysis(casedata, tap_changers_data, phase_shifters_data)
%taps_and_phases_analysis() caltulates taps and phases by iterating pwoer flows
% [RESULTS, ITERATION_COUNT,SUCCESS] = taps_and_phases_analysis(casedata, tap_changers_data, phase_shifters_data)
%
% Runs a Newtonian routine that calcualtes the taps and phases from
% specified nodal voltages and active power flow criteria. Taps and
% phases not specified in the data specifically are fixed. Initial taps
% and phases are read from data.
%
% Inputs (all are mandatory), if there are no transformers use []:
% casedata : A string containing the name of the file with the case
% data (e.g. input = 'case300')
% tap_changers_data : Data of all the tap changers connected into the
% system in the form of a matrix, the first column contains branches
% the second contains regulated nodes and the third the voltage
% magnitude in p.u.. For every transformer it is a line in the form
% of tap_changers_data = [branch , controlled node , votlage in p.u.]
%
% phase_shifters_data : Data of all the phase shifters connected into
% the system. The first column contains branches and the second
% contains the active power flow into the line at the first node of
% the branch of the device's insertion in p.u.. Every line is in the
% form of phase_shifters_data = [branch , active power in p.u.]
%
% Outputs (all are optional):
% RESULTS : results struct
% ITERATION_COUNT : number of successive load flow calcualtions ran
% which is also the iteration count with this method
% SUCCESS : 0 if failed to converge, 1 if converged

% Example:
% results = taps_and_phases_analysis('case300',[1,9001,.95],[100,0]);
% By Gorazd Bone, Faculty of Electrical Engineerinc, Ljubljana


st_tap_ch = size(tap_changers_data,1); %number of tap changers
st_pha_sh = size(phase_shifters_data,1); %number of phase shifters

% U %the control inputs are defined from initial conditions
H = zeros(st_tap_ch + st_pha_sh,1); %the criteria vector
itandp = zeros(st_tap_ch + st_pha_sh,1); %initial taps and phases

%% 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, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[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; %#ok<*NASGU,*ASGLU>


%% Read data
mpc_ext = loadcase(casedata);
mpc_int = ext2int(mpc_ext);

%% Index reordering for tap changer control bus
ind_voz = [(1:size(mpc_ext.bus(:,BUS_I),1))' , mpc_ext.bus(:,BUS_I)];
for k=1: st_tap_ch
tap_changers_data(k,2) = ind_voz(ind_voz(:,2) == tap_changers_data(k,2), 1);
end

%% Initial state determination
mpopt = mpoption('out.all', 0, 'verbose', 0);
results = ext2int(runpf(mpc_int,mpopt));
results.branch(results.branch(:,TAP)==0,TAP) = 1; % 0 initial tap is regarded as 1
Sbazni = mpc_int.baseMVA; %razmerje kako so podane moci v vektorju results.branch
if size(tap_changers_data,1>0)
itandp(1 : st_tap_ch) = results.branch(tap_changers_data(:,1),TAP);
H(1 : st_tap_ch) = results.bus(tap_changers_data(:,2),VM) - tap_changers_data(:,3);
end
if size(phase_shifters_data,1>0)
itandp(st_tap_ch + 1 : st_tap_ch + st_pha_sh) = results.branch(phase_shifters_data(:,1),SHIFT)/180*pi;
H(st_tap_ch + 1 : st_tap_ch + st_pha_sh) = results.branch(phase_shifters_data(:,1),PF)/Sbazni - phase_shifters_data(:,2);
end
U = itandp; % the input vector

%% Read node types
pqtip = results.bus(:,BUS_TYPE)==PQ;
pvtip = results.bus(:,BUS_TYPE)==PV;
tip_vozlisc = pqtip*2 + pvtip; % nodal type, 2 for PQ 1 for PV, corresponds to the nubmer of equations

%% Iterate taps and phases
stevilo_iteracij = 0; %iteration count
while (norm(H,2)>1e-5) && (stevilo_iteracij<40)
J = shiftJac_taps_phases(results); % jacobian matrix
obcutljivosti = sensit_to_taps_and_phases( results, tap_changers_data, phase_shifters_data, J , tip_vozlisc);% system sensitivity calculation
U = U - obcutljivosti\H;%step by Newton's method

if st_tap_ch
results.branch(tap_changers_data(:,1),TAP) = U(1:st_tap_ch); %apply new taps
end
if st_pha_sh
results.branch(phase_shifters_data(:,1),SHIFT) = U(st_tap_ch + 1 : st_tap_ch + st_pha_sh)*180/pi;%apply new phases
end

results = ext2int(runpf(results,mpopt)); % rerun load flow with new values
if size(tap_changers_data,1>0)
H(1 : st_tap_ch) = results.bus(tap_changers_data(:,2),VM) - tap_changers_data(:,3); % criteria for tap changers
end
if size(phase_shifters_data,1>0)
H(st_tap_ch + 1 : st_tap_ch + st_pha_sh) = results.branch(phase_shifters_data(:,1),PF)/Sbazni - phase_shifters_data(:,2);% criteria for phase shifters
end
stevilo_iteracij = stevilo_iteracij + 1;
end

results = int2ext(results);

success = (~(norm(H,2)>1e-5));