Skip to content

Commit 60e4bb7

Browse files
GorazdBrdzman
GorazdB
authored and
rdzman
committed
Tap and Phases calculation for transformers.
> > Four files added, one of them is an example file. > Call taps_and_phases_analysis to execute. > > Files sensit_to_taps_and_phases and shiftJac_taps_phases must be > included in the same folder or in path.
1 parent a3bc668 commit 60e4bb7

4 files changed

Lines changed: 357 additions & 0 deletions

File tree

‎lib/example_tap_phase_calc.m‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
%An example of running the tap changer and phase shifter analysis on the
2+
%IEEE 300 bus system with two tap changers and two phase shifters
3+
4+
5+
%% Define case data
6+
file = 'case300';
7+
8+
%% 1. Define the transformers
9+
%tap_changers_data = [line of insertion, control node, control voltage]
10+
tap_changers_data = [
11+
1,9001,.95
12+
2,9005, 1
13+
];
14+
15+
%phase_shifters_data = [line of insertion, control power]
16+
phase_shifters_data = [
17+
100 , 0
18+
110 , 1
19+
];
20+
21+
%% Run the analysis
22+
[results, stevilo_iteracij, success] = taps_and_phases_analysis(file, tap_changers_data, phase_shifters_data);
23+

‎lib/sensit_to_taps_and_phases.m‎

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
function [ sensitivity_of_H_to_U_taps_and_phasess ] = sensit_to_taps_and_phases( results, tap_changers, phase_shifters, Jacobian , tip_vozlisc)
2+
%sensit_to_taps_and_phases Calculates the first orded sensitivities of
3+
%regulated active powers and nodal voltages to taps and phases.
4+
% [sensitivity_of_H_to_U_taps_and_phasess] = sensit_to_taps_and_phases(results, tap_changers, phase_shifters, Jacobian , node_types)
5+
%
6+
% Inputs :
7+
% results: results from a runpf(), interanly idnexed (use ext2int())
8+
% tap_changers: tap changers data as in taps_and_phases_analysis()
9+
% phase_shifters: phase shifter data as in taps_and_phases_analysis()
10+
% Jacobian: power missmatch Jacobian matrix with the appropriate
11+
% ordering, rows go P1, Q1, P2, Q2...., columns go V1, delta1, V2,
12+
% delta2.....
13+
% tip_vozlisc: node types, 2 for PQ, 1 for PV, same as nubmer of
14+
% equations for each node, used for indexing
15+
%
16+
% Output :
17+
% sensitivity_of_H_to_U_taps_and_phasess : first order sensitivities
18+
% of p.u. changes in nodal votlages and line active powers to changes
19+
% in calculated taps and phases of transformers.
20+
%
21+
% by Gorazd Bone, Faculty of Electrical Engineering, Ljubljana
22+
23+
ind_node_eq = tip_vozlisc * 0;
24+
for k=1:size(tip_vozlisc,1)
25+
ind_node_eq(k) = sum(tip_vozlisc(1:k-1)) + 1;
26+
end
27+
28+
stevilo_tap_ch = size(tap_changers,1);
29+
stevilo_pha_sh = size(phase_shifters,1);
30+
31+
%% define named indices into bus, gen, branch matrices
32+
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
33+
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; %#ok<*NASGU,*ASGLU>
34+
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
35+
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
36+
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
37+
38+
%% Define Partial derivative matrices
39+
% G is the LF problem power injection missmatch form (sum of nodal powers)
40+
GpoU = zeros(sum(tip_vozlisc) , stevilo_tap_ch + stevilo_pha_sh); % derivative of LF problem w.r.t. taps & phases
41+
GpoX = Jacobian; % derivative of LF problem w.r.t. voltages and angles
42+
HpoU = zeros(stevilo_tap_ch + stevilo_pha_sh); %derivative of criterion function w.r.t. taps & phases
43+
HpoX = zeros(stevilo_tap_ch + stevilo_pha_sh , sum(tip_vozlisc));%derivative of criterion function w.r.t. voltages and angles
44+
45+
for k=1:stevilo_tap_ch % partial derivatives for taps
46+
veja = tap_changers(k,1);%tap branch
47+
reg_voz = tap_changers(k,2); %ragulated node
48+
voz1 = results.branch(veja,F_BUS);
49+
voz2 = results.branch(veja,T_BUS);
50+
indeks_P1 = ind_node_eq(voz1); %active power and volt. amp index for f node
51+
indeks_P2 = ind_node_eq(voz2); %active power and volt. amp index for t node
52+
indeks_reg_voz = ind_node_eq(reg_voz);%volt. amp index for reg. node
53+
54+
Uvozi = results.bus(voz1,VM);
55+
Uvozj = results.bus(voz2,VM);
56+
di = results.bus(voz1,VA)/180*pi;
57+
dj = results.bus(voz2,VA)/180*pi;
58+
Z = results.branch(veja,BR_R) + 1i*results.branch(veja,BR_X);
59+
Bsh = results.branch(veja,BR_B);
60+
G = real(1/Z);
61+
B = imag(1/Z);
62+
tap = results.branch(veja,TAP);
63+
ph = results.branch(veja,SHIFT)/180*pi;
64+
65+
if tip_vozlisc(reg_voz) == 1 % tap changer regulating PV node
66+
error_string = strcat(2 , 'Tap changer' , int2str(k) , 'regulating the voltage of a PV node');
67+
errordlg(error_string );
68+
return
69+
end
70+
71+
dPijdtap = (Uvozi*(-2*G*Uvozi+G*tap*Uvozj*cos(di-dj-ph)+B*tap*Uvozj*sin(di-dj-ph)))/tap^3;
72+
dPjidtap = (Uvozi*Uvozj*(G*cos(di-dj-ph)-B*sin(di-dj-ph)))/tap^2;
73+
dQijdtap = (Uvozi*((2*B+Bsh)*Uvozi-B*tap*Uvozj*cos(di-dj-ph)+G*tap*Uvozj*sin(di-dj-ph)))/tap^3;
74+
dQjidtap = -((Uvozi*Uvozj*(B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap^2);
75+
76+
%% dG/dU
77+
GpoU(indeks_P1 , k) = - dPijdtap; %dGfrom / dtap
78+
if tip_vozlisc(voz1) == 2 %PQ type of node
79+
GpoU(indeks_P1 + 1 , k) = - dQijdtap;
80+
end
81+
GpoU(indeks_P2 , k) = - dPjidtap;%dGto / dtap
82+
if tip_vozlisc(voz2) == 2 %PQ type of node
83+
GpoU(indeks_P2 + 1 , k) = - dQjidtap;
84+
end
85+
86+
%% dH/dU
87+
HpoU(k,k) = 0;%odvod vozliscne napetosti po tap-u je 0
88+
if stevilo_pha_sh
89+
if find(phase_shifters(:,1) == veja)
90+
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
91+
end
92+
end
93+
94+
%% dH/dX
95+
HpoX(k,indeks_reg_voz) = 1; %dRegNodeVolt / dVoltAmp = unity
96+
end
97+
98+
for k=1:stevilo_pha_sh % partial derivatives for phases
99+
veja = phase_shifters(k,1);%tap branch
100+
voz1 = results.branch(veja,F_BUS);
101+
voz2 = results.branch(veja,T_BUS);
102+
indeks_P1 = ind_node_eq(voz1); %active power and volt. amp index for f node
103+
indeks_P2 = ind_node_eq(voz2); %active power and volt. amp index for t node
104+
105+
Uvozi = results.bus(voz1,VM);
106+
Uvozj = results.bus(voz2,VM);
107+
di = results.bus(voz1,VA)/180*pi;
108+
dj = results.bus(voz2,VA)/180*pi;
109+
Z = results.branch(veja,BR_R) + 1i*results.branch(veja,BR_X);
110+
Bsh = results.branch(veja,BR_B);
111+
G = real(1/Z);
112+
B = imag(1/Z);
113+
tap = results.branch(veja,TAP);
114+
ph = results.branch(veja,SHIFT)/180*pi;
115+
116+
dPijdph = (Uvozi*Uvozj*(B*cos(di-dj-ph)-G*sin(di-dj-ph)))/tap;
117+
dPjidph = -((Uvozi*Uvozj*(B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap);
118+
dPijUi = (2*G*Uvozi-tap*Uvozj*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap^2;
119+
dPijUj = -((Uvozi*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap);
120+
dPijdi = (Uvozi*Uvozj*(-B*cos(di-dj-ph)+G*sin(di-dj-ph)))/tap;
121+
dPijdj = (Uvozi*Uvozj*(B*cos(di-dj-ph)-G*sin(di-dj-ph)))/tap;
122+
dQijdph = (Uvozi*Uvozj*(G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap;
123+
dQjidph = (Uvozi*Uvozj*(-G*cos(di-dj-ph)+B*sin(di-dj-ph)))/tap;
124+
125+
%% dG/dU
126+
GpoU(indeks_P1, k + stevilo_tap_ch) = - dPijdph;%dGfrom / dphase
127+
if tip_vozlisc(voz1) == 2 %PQ type of node
128+
GpoU(indeks_P1 + 1, k + stevilo_tap_ch) = - dQijdph;%dGfrom / dphase
129+
end
130+
GpoU(indeks_P2, k + stevilo_tap_ch) = - dPjidph;%dGto / dphase
131+
if tip_vozlisc(voz2) == 2 %PQ type of node
132+
GpoU(indeks_P2 + 1, k + stevilo_tap_ch) = - dQjidph;%dGto / dphase
133+
end
134+
135+
%% dH/dU
136+
HpoU(k + stevilo_tap_ch, k + stevilo_tap_ch) = dPijdph;%dLinePower / dphase
137+
138+
%% dH/dX
139+
if tip_vozlisc(voz1) == 2 %from bus is PQ type
140+
HpoX(k + stevilo_tap_ch, indeks_P1) = dPijUi;%dLinePower / dVoltAmp
141+
HpoX(k + stevilo_tap_ch, indeks_P1+1) = dPijdi;%dLinePower / dVoltAng
142+
else %from bus is PV type
143+
HpoX(k + stevilo_tap_ch, indeks_P1) = dPijdi;%dLinePower / dVoltAng
144+
end
145+
146+
if tip_vozlisc(voz2) == 2 %same thing as abocve for to bus
147+
HpoX(k + stevilo_tap_ch, indeks_P2) = dPijUj;
148+
HpoX(k + stevilo_tap_ch, indeks_P2+1) = dPijdj;
149+
else
150+
HpoX(k + stevilo_tap_ch, indeks_P2) = dPijdj;
151+
end
152+
153+
end
154+
155+
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
156+
en

‎lib/shiftJac_taps_phases.m‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
function [ J ] = shiftJac_taps_phases( results )
2+
%shiftJac_taps_phases Constructs the Jacobian matrix using Matpower
3+
%supplied makeJac and reorderes it for taps and phases analysis.
4+
% [J ] = shiftJac_taps_phases( results )
5+
% Inputs :
6+
% results: results from a runpf(), interanly idnexed (use ext2int)
7+
%
8+
% Output :
9+
% J : Jacobian matrix for power mismatch load flow, ordereing
10+
% appropriate fo sensit_to_taps_and_phases()
11+
%
12+
% by Gorazd Bone, Faculty of Electrical Engineering, Ljubljana
13+
14+
%% define named indices into bus, gen, branch matrices
15+
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
16+
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
17+
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
18+
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
19+
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; %#ok<*ASGLU,*NASGU>
20+
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
21+
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
22+
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
23+
24+
pq = results.bus(:,BUS_TYPE)==PQ;
25+
pv = results.bus(:,BUS_TYPE)==PV;
26+
ref = find(pq + pv == 0);
27+
28+
tip_vozlisc = pq * 2 + pv;
29+
ind_node_eq = tip_vozlisc * 0;
30+
for k=1:size(tip_vozlisc,1)
31+
ind_node_eq(k) = sum(tip_vozlisc(1:k-1)) + 1;
32+
end
33+
loc_P_pv = ind_node_eq(pv==1);
34+
loc_P_pq = ind_node_eq(pq==1);
35+
loc_Q_pq = loc_P_pq+1;
36+
37+
38+
per_rows = sparse([loc_P_pv;loc_P_pq;loc_Q_pq] , 1:sum(pq * 2 + pv),1);
39+
per_cols = sparse(1:sum(pq * 2 + pv) , [loc_P_pv;loc_Q_pq;loc_P_pq] , 1);
40+
%
41+
% MatPower Jacobi has the following oredering for the missmatch:
42+
%
43+
% P_vect_of_pvnodes
44+
% P_vect_of_pqnodes
45+
% Q_vect_of_pqnodes
46+
%
47+
% and for the unknowns:
48+
%
49+
% angle_vect_of_pvnodes
50+
% angle_vect_of_pqnodes
51+
% voltage_vect_of_pqnodes
52+
%
53+
% the used Jacobi has the ordering node-by-node, for the missmatch
54+
%
55+
% P_of_node1
56+
% Q_of_node1 (if PQ)
57+
% P_of_node2
58+
% Q_of_node2 (if PQ)
59+
%
60+
% and for the unknowns:
61+
%
62+
% voltage_of_node1 (if PQ)
63+
% anle_of_node1
64+
% voltage_of_node2 (if PQ)
65+
% anle_of_node2 (if PQ)
66+
%
67+
J = - per_rows * makeJac(results.baseMVA, results.bus, results.branch, results.gen) * per_cols;
68+

‎lib/taps_and_phases_analysis.m‎

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
function [results, stevilo_iteracij, success] = taps_and_phases_analysis(casedata, tap_changers_data, phase_shifters_data)
2+
%taps_and_phases_analysis() caltulates taps and phases by iterating pwoer flows
3+
% [RESULTS, ITERATION_COUNT,SUCCESS] = taps_and_phases_analysis(casedata, tap_changers_data, phase_shifters_data)
4+
%
5+
% Runs a Newtonian routine that calcualtes the taps and phases from
6+
% specified nodal voltages and active power flow criteria. Taps and
7+
% phases not specified in the data specifically are fixed. Initial taps
8+
% and phases are read from data.
9+
%
10+
% Inputs (all are mandatory), if there are no transformers use []:
11+
% casedata : A string containing the name of the file with the case
12+
% data (e.g. input = 'case300')
13+
% tap_changers_data : Data of all the tap changers connected into the
14+
% system in the form of a matrix, the first column contains branches
15+
% the second contains regulated nodes and the third the voltage
16+
% magnitude in p.u.. For every transformer it is a line in the form
17+
% of tap_changers_data = [branch , controlled node , votlage in p.u.]
18+
%
19+
% phase_shifters_data : Data of all the phase shifters connected into
20+
% the system. The first column contains branches and the second
21+
% contains the active power flow into the line at the first node of
22+
% the branch of the device's insertion in p.u.. Every line is in the
23+
% form of phase_shifters_data = [branch , active power in p.u.]
24+
%
25+
% Outputs (all are optional):
26+
% RESULTS : results struct
27+
% ITERATION_COUNT : number of successive load flow calcualtions ran
28+
% which is also the iteration count with this method
29+
% SUCCESS : 0 if failed to converge, 1 if converged
30+
31+
% Example:
32+
% results = taps_and_phases_analysis('case300',[1,9001,.95],[100,0]);
33+
% By Gorazd Bone, Faculty of Electrical Engineerinc, Ljubljana
34+
35+
36+
st_tap_ch = size(tap_changers_data,1); %number of tap changers
37+
st_pha_sh = size(phase_shifters_data,1); %number of phase shifters
38+
39+
% U %the control inputs are defined from initial conditions
40+
H = zeros(st_tap_ch + st_pha_sh,1); %the criteria vector
41+
itandp = zeros(st_tap_ch + st_pha_sh,1); %initial taps and phases
42+
43+
%% define named indices into bus, gen, branch matrices
44+
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
45+
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
46+
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
47+
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
48+
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
49+
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
50+
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
51+
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; %#ok<*NASGU,*ASGLU>
52+
53+
54+
%% Read data
55+
mpc_ext = loadcase(casedata);
56+
mpc_int = ext2int(mpc_ext);
57+
58+
%% Index reordering for tap changer control bus
59+
ind_voz = [(1:size(mpc_ext.bus(:,BUS_I),1))' , mpc_ext.bus(:,BUS_I)];
60+
for k=1: st_tap_ch
61+
tap_changers_data(k,2) = ind_voz(ind_voz(:,2) == tap_changers_data(k,2), 1);
62+
end
63+
64+
%% Initial state determination
65+
mpopt = mpoption('out.all', 0, 'verbose', 0);
66+
results = ext2int(runpf(mpc_int,mpopt));
67+
results.branch(results.branch(:,TAP)==0,TAP) = 1; % 0 initial tap is regarded as 1
68+
Sbazni = mpc_int.baseMVA; %razmerje kako so podane moci v vektorju results.branch
69+
if size(tap_changers_data,1>0)
70+
itandp(1 : st_tap_ch) = results.branch(tap_changers_data(:,1),TAP);
71+
H(1 : st_tap_ch) = results.bus(tap_changers_data(:,2),VM) - tap_changers_data(:,3);
72+
end
73+
if size(phase_shifters_data,1>0)
74+
itandp(st_tap_ch + 1 : st_tap_ch + st_pha_sh) = results.branch(phase_shifters_data(:,1),SHIFT)/180*pi;
75+
H(st_tap_ch + 1 : st_tap_ch + st_pha_sh) = results.branch(phase_shifters_data(:,1),PF)/Sbazni - phase_shifters_data(:,2);
76+
end
77+
U = itandp; % the input vector
78+
79+
%% Read node types
80+
pqtip = results.bus(:,BUS_TYPE)==PQ;
81+
pvtip = results.bus(:,BUS_TYPE)==PV;
82+
tip_vozlisc = pqtip*2 + pvtip; % nodal type, 2 for PQ 1 for PV, corresponds to the nubmer of equations
83+
84+
%% Iterate taps and phases
85+
stevilo_iteracij = 0; %iteration count
86+
while (norm(H,2)>1e-5) && (stevilo_iteracij<40)
87+
J = shiftJac_taps_phases(results); % jacobian matrix
88+
obcutljivosti = sensit_to_taps_and_phases( results, tap_changers_data, phase_shifters_data, J , tip_vozlisc);% system sensitivity calculation
89+
U = U - obcutljivosti\H;%step by Newton's method
90+
91+
if st_tap_ch
92+
results.branch(tap_changers_data(:,1),TAP) = U(1:st_tap_ch); %apply new taps
93+
end
94+
if st_pha_sh
95+
results.branch(phase_shifters_data(:,1),SHIFT) = U(st_tap_ch + 1 : st_tap_ch + st_pha_sh)*180/pi;%apply new phases
96+
end
97+
98+
results = ext2int(runpf(results,mpopt)); % rerun load flow with new values
99+
if size(tap_changers_data,1>0)
100+
H(1 : st_tap_ch) = results.bus(tap_changers_data(:,2),VM) - tap_changers_data(:,3); % criteria for tap changers
101+
end
102+
if size(phase_shifters_data,1>0)
103+
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
104+
end
105+
stevilo_iteracij = stevilo_iteracij + 1;
106+
end
107+
108+
results = int2ext(results);
109+
110+
success = (~(norm(H,2)>1e-5));

0 commit comments

Comments
 (0)