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.

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);
+

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

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;
+

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));