diff --git a/lib/example_tap_phase_calc.m b/lib/example_tap_phase_calc.m new file mode 100644 index 00000000..fe05180e --- /dev/null +++ b/lib/example_tap_phase_calc.m @@ -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); + diff --git a/lib/sensit_to_taps_and_phases.m b/lib/sensit_to_taps_and_phases.m new file mode 100644 index 00000000..8f30ace4 --- /dev/null +++ b/lib/sensit_to_taps_and_phases.m @@ -0,0 +1,156 @@ +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 +en \ No newline at end of file diff --git a/lib/shiftJac_taps_phases.m b/lib/shiftJac_taps_phases.m new file mode 100644 index 00000000..ef99bbce --- /dev/null +++ b/lib/shiftJac_taps_phases.m @@ -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; + diff --git a/lib/taps_and_phases_analysis.m b/lib/taps_and_phases_analysis.m new file mode 100644 index 00000000..770200a0 --- /dev/null +++ b/lib/taps_and_phases_analysis.m @@ -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)); \ No newline at end of file