clear all
clc

% ----------------------------------------------------------------- %

% Sample time

st = 1e-3;

% ----------------------------------------------------------------- %

% Network Parameters
fn = 50;
wn = 2*pi*fn;
baseMVA = 100e6;
basekV = 400e3;
HR = 1e6;
HRbase = HR/baseMVA;

% ----------------------------------------------------------------- %

% 22 kV Generator Parameters 

% Stator
Vg = 22e3;
Sgob = 355e6;
Pgen = 320e6;
Qgmin = -42e6;
Qgmax = 131e6;
Qpmin = -46e6;
Qpmax = 148e6;
Qkmin = -155e6;
Qkmax = 150e6;
% Field
Vfflg = 218;
Ifflg = 2430;
Vfnlg = 69;
Ifnlg = 1070;
Polesg = 18;
Inertiag = 4;
% Reactances pu
Rsg = 0; % Stator resistance (assumption)
Xdg = 1.5; % Direct axis reactance
Xqg = 0.88; % Quadrature axis reactance
Xlg = 0.15; % Leakage reactance (assumption)
X1dg = 0.42; % Transient reactance
X2dg = 0.24; % Subtransient reactance
X0g = 0.10; % Zero sequence reactance 
X2g = 0.26; % Negative sequence reactance
T2d0g = 0.05; %Subtransient time constant direct axis (Assumption)
T2q0g = 0.09; %Subtransient time constant quadrature axis (Assumption)

% Base Values Calculation
% Stator
Sgbase = Sgob;
Vgbase = Vg*(sqrt(2)/sqrt(3));
Igbase = (Sgbase*sqrt(2))/(Vg*sqrt(3));
Zgbase = Vgbase/Igbase;
% Lgbase = Zgbase/wn; 
% Rotor Values Refered to stator
Ifgbase = Ifnlg*(Xdg - Xlg); 
Vfgbase = Sgbase/Ifgbase;
Zfgbase = Vfgbase/Ifgbase;
% Lfgbase = Zfgbase/wn;
% Ratio field - stator
ratiog = (2/3)*(Ifgbase/Igbase);
%Rotor Internal parameters 
Rfsg = (3/2)*(Vfnlg/Ifnlg)*(ratiog^2);
% Llfsg = Xlg*Lgbase;
Rfsgpu = Rfsg/Zgbase;
% Llfsgpu = Xlg;

syms lfdg l2qg

% Synchronous Machine pu Standard missing parameters calculation
Xfdg = solve(X1dg == Xlg + ((Xdg - Xlg)*lfdg)/((Xdg - Xlg) + lfdg), lfdg); %Field reactance direct axis
X2qg = double(solve(X2g == 0.5*(X2dg + l2qg), l2qg)); % Subtransient reactance quadrature axis
T1d0g = double((1/wn)*(((Xdg - Xlg) + Xfdg)/(Rfsgpu))); %Transient time constant direct axis

% Generator Base values 
factorGen = baseMVA/Sgob;
Vgbase1 = 1;
Sgbase1 = Sgob;
Pgenbase = Pgen/Sgob;
Qgminbase = Qgmin/Sgob;
Qgmaxbase = Qgmax/Sgob;
Xdgbase = Xdg*factorGen; 
Xqgbase = Xqg*factorGen; 
Xlgbase = Xlg*factorGen; 
X1dgbase = X1dg*factorGen; 
X2dgbase = X2dg*factorGen; 
X2qgbase = X2qg*factorGen;

% ----------------------------------------------------------------- %

% 10.5 kV Condensator Parameters

% Stator
Vc = 10.5e3;
Scob = 355e6;
Pcon = 0.5e6;
Qcmin = -55.28e6;
Qcmax = 110.59e6;
% Field
Vfflc = 218;
Ifflc = 2430;
Vfnlc = Vfnlg/1.6;
Ifnlc = Ifnlg/1.6;
Polesc = Polesg/3;
Inertiac = Inertiag/3;
% Reactances pu
Rsc = 0 ; % Stator resistance (assumption)
Xdc = 2.2; % Direct axis reactance
Xqc = 1.15; % Quadrature axis reactance
Xlc = 0.15; % Leakage reactance (assumption)
X1dc = 0.55; % Transient reactance
X2dc = 0.40; % Subtransient reactance
X0c = 0.20; % Zero sequence reactance
X2c = 0.40; % Negative sequence reactance
T2d0c = 0.035; %Subtransient time constant direct axis (Assumption)
T2q0c = 0.035; %Subtransient time constant quadrature axis (Assumption)

% Base Values Calculation
% Stator
Scbase = Scob;
Vcbase = (Vc*sqrt(2))/sqrt(3);
Icbase = (Scbase*sqrt(2))/(Vc*sqrt(3));
Zcbase = Vcbase/Icbase;
% Lcbase = Zcbase/wn; 
% Rotor Values Refered to stator
Ifcbase = Ifnlc*(Xdc - Xlc); 
Vfcbase = Scbase/Ifcbase;
Zfcbase = Vfcbase/Ifcbase;
% Lfcbase = Zfcbase/wn;
% Ratio field - stator
ratioc = (2/3)*(Ifcbase/Icbase);
%Rotor Internal parameters 
Rfsc = (3/2)*(Vfnlc/Ifnlc)*(ratioc^2);
% Llfsc = Xlc*Lcbase;
Rfscpu = Rfsc/Zcbase;
% Llfsgcu = Xlc;

syms lfdc l2qc 

% Synchronous Machine pu Standard parameters calculation
Xfdc = solve(X1dc == Xlc + ((Xdc - Xlc)*lfdc)/((Xdc - Xlc) + lfdc), lfdc); %Field reactance direct axis
X2qc = double(solve(X2c == 0.5*(X2dc + l2qc), l2qc)); % Subtransient reactance quadrature axis
T1d0c = double((1/wn)*(((Xdc - Xlc) + Xfdc)/(Rfscpu))); %Transient time constant direct axis

% Condenser Base values 
factorCon = baseMVA/Scob;
Vcbase1 = 1;
Scbase1 = Scob;
Pconbase = Pcon;
Qcminbase = Qcmin/Scob;
Qcmaxbase = Qcmax/Scob;
Xdcbase = Xdc*factorCon; 
Xqcbase = Xqc*factorCon; 
Xlcbase = Xlc*factorCon; 
X1dcbase = X1dc*factorCon; 
X2dcbase = X2dc*factorCon; 
X2qcbase = X2qc*factorCon;

% ----------------------------------------------------------------- %

% Two winding transformer

%Operational parameters
St12 = 360e6;
Up12 = 420e3;
Us12 = 22e3;
% T1, T2 SC and OC test parameters
Po12 = 233e3;
Io12 = 0.4;
uk12 = 0.14;
Pk12 = 750e3;
% Internal parameters calculation
r12 = (Pk12)/(St12);
x12 = sqrt((uk12)^2 - (r12)^2);
% Transformer parameters
T12r = r12/2;
T12x = x12/2;
T12b = Po12/St12;

%Windings Transformer base values
factorT2w = baseMVA/St12;
factorkV = (Up12/basekV)^2;
T12xbase = T12x*factorT2w;

%Two windings Transformer base values
T2Wxbase = uk12*factorT2w*factorkV;

% ----------------------------------------------------------------- %

% Three winding transformers

%Operational parameters
St40x1 = 350e6;
St40x3 = 100e6;
Up40x = 420e3;
Us40x = 120e3;
Ut40x = 10.5e3;
factorT3w1 = baseMVA/St40x1;
factorT3w2 = baseMVA/St40x3;

%Transformer T401 SC and OC test parameters
% No load losses
T401Po = 84.05e3;
% Impedance base St40x1, SC losses base St40x1
T401uk12 = 13.92/100;
T401Pk12 = 559.21e3;
% Impedance base St40x1, SC losses base St40x3
T401uk13 = 36.65/100;
T401Pk13 = 141.75e3;
% Impedance base St40x1, SC losses base St40x3
T401uk23 = 19.64/100;
T401Pk23 = 131.96e3;

% T401 Internal parameters calculation
T401b = T401Po/St40x1;
T401r12 = (T401Pk12)/(St40x1);
T401x12 = sqrt((T401uk12)^2 - (T401r12)^2);
T401r13 = (T401Pk13)/(St40x3);
T401x13 = sqrt((T401uk13)^2 - (T401r13)^2);
T401r23 = (T401Pk23)/(St40x3);
T401x23 = sqrt((T401uk23)^2 - (T401r23)^2);

% T401 Transformer parameters
T401r1 = (T401r12 + T401r13 - T401r23)/2;
T401r2 = (T401r12 - T401r13 + T401r23)/2;
T401r3 = (-T401r12 + T401r13 + T401r23)/2;
T401x1 = (T401x12 + T401x13 - T401x23)/2;
T401x2 = (T401x12 - T401x13 + T401x23)/2;
T401x3 = (-T401x12 + T401x13 + T401x23)/2;

%T401 windings Transformer base values
T401z1 = (T401uk12 + T401uk13 - T401uk23)/2;
T401z2 = (T401uk12 - T401uk13 + T401uk23)/2;
T401z3 = (-T401uk12 + T401uk13 + T401uk23)/2;

T401z12 = (T401z1*T401z2 + T401z1*T401z3 + T401z2*T401z3)/T401z3;
T401z23 = (T401z1*T401z2 + T401z1*T401z3 + T401z2*T401z3)/T401z1;
T401z13 = (T401z1*T401z2 + T401z1*T401z3 + T401z2*T401z3)/T401z2;

T401x12base = T401z12*factorT3w1*factorkV;
T401x23base = T401z23*factorT3w1;
T401x13base = T401z13*factorT3w1*factorkV;

%Transformer T402 parameters SC and OC test parameters
% No load losses
T402Po = 88.34e3;
% Impedance base St40x1, SC losses base St40x1
T402uk12 = 13.85/100;
T402Pk12 = 542.90e3;
% Impedance base St40x1, SC losses base St40x3
T402uk13 = 36.42/100;
T402Pk13 = 138.04e3;
% Impedance base St40x1, SC losses base St40x3
T402uk23 = 19.56/100;
T402Pk23 = 131.07e3;

% T402 Internal parameters calculation
T402b = T402Po/St40x1;
T402r12 = (T402Pk12)/(St40x1);
T402x12 = sqrt((T402uk12)^2 - (T402r12)^2);
T402r13 = (T402Pk13)/(St40x3);
T402x13 = sqrt((T402uk13)^2 - (T402r13)^2);
T402r23 = (T402Pk23)/(St40x3);
T402x23 = sqrt((T402uk23)^2 - (T402r23)^2);

% T402 Transformer parameters
T402r1 = (T402r12 + T402r13 - T402r23)/2;
T402r2 = (T402r12 - T402r13 + T402r23)/2;
T402r3 = (-T402r12 + T402r13 + T402r23)/2;
T402x1 = (T402x12 + T402x13 - T402x23)/2;
T402x2 = (T402x12 - T402x13 + T402x23)/2;
T402x3 = (-T402x12 + T402x13 + T402x23)/2;

%T402 windings Transformer base values
T402z1 = (T402uk12 + T402uk13 - T402uk23)/2;
T402z2 = (T402uk12 - T402uk13 + T402uk23)/2;
T402z3 = (-T402uk12 + T402uk13 + T402uk23)/2;

T402z12 = (T402z1*T402z2 + T402z1*T402z3 + T402z2*T402z3)/T402z3;
T402z23 = (T402z1*T402z2 + T402z1*T402z3 + T402z2*T402z3)/T402z1;
T402z13 = (T402z1*T402z2 + T402z1*T402z3 + T402z2*T402z3)/T402z2;

T402x12base = T402z12*factorT3w1*factorkV;
T402x23base = T402z23*factorT3w1;
T402x13base = T402z13*factorT3w1*factorkV;

%Transformer T403 parameters
% No load losses
T403Po = 94.48e3;
% Impedance base St40x1, SC losses base St40x1
T403uk12 = 14.06/100;
T403Pk12 = 588e3;
% Impedance base St40x1, SC losses base St40x3
T403uk13 = 36.50/100;
T403Pk13 = 132e3;
% Impedance base St40x1, SC losses base St40x3
T403uk23 = 19.32/100;
T403Pk23 = 133e3;

% Internal parameters calculation
T403b = T403Po/St40x1;
T403r12 = (T403Pk12)/(St40x1);
T403x12 = sqrt((T403uk12)^2 - (T403r12)^2);
T403r13 = (T403Pk13)/(St40x3);
T403x13 = sqrt((T403uk13)^2 - (T403r13)^2);
T403r23 = (T403Pk23)/(St40x3);
T403x23 = sqrt((T403uk23)^2 - (T403r23)^2);

% T403 Transformer parameters
T403r1 = (T403r12 + T403r13 - T403r23)/2;
T403r2 = (T403r12 - T403r13 + T403r23)/2;
T403r3 = (-T403r12 + T403r13 + T403r23)/2;
T403x1 = (T403x12 + T403x13 - T403x23)/2;
T403x2 = (T403x12 - T403x13 + T403x23)/2;
T403x3 = (-T403x12 + T403x13 + T403x23)/2;

%T403 windings Transformer base values
T403z1 = (T403uk12 + T403uk13 - T403uk23)/2;
T403z2 = (T403uk12 - T403uk13 + T403uk23)/2;
T403z3 = (-T403uk12 + T403uk13 + T403uk23)/2;

T403z12 = (T403z1*T403z2 + T403z1*T403z3 + T403z2*T403z3)/T403z3;
T403z23 = (T403z1*T403z2 + T403z1*T403z3 + T403z2*T403z3)/T403z1;
T403z13 = (T403z1*T403z2 + T403z1*T403z3 + T403z2*T403z3)/T403z2;

T403x12base = T403z12*factorT3w1*factorkV;
T403x23base = T403z23*factorT3w1;
T403x13base = T403z13*factorT3w1*factorkV;

% ----------------------------------------------------------------- %

% Transmission Line Parameters

VTL = 400e3;
ZTL = (VTL^2)/baseMVA;
r1TL = 1.37;
x1TL = 16.94;
r1TLbase = r1TL/ZTL;
x1TLbase = x1TL/ZTL;

% ----------------------------------------------------------------- %

% Network equivalent
MVAsc = 15118*10^6;
XR = 15;
VTS = 400e3;
ZNE = (VTS^2)/baseMVA;

XTS = (VTS^2)/MVAsc;
RTS = XTS/XR;

XTSbase = XTS/ZNE;  
RTSbase = RTS/ZNE;

% ----------------------------------------------------------------- %

% Reactor
%VSR Inductance Calculation
MVAR = (1e6)*[60.00; 61.40; 62.80; 64.30; 65.80; 67.30; 68.90; 70.60; 72.30; 74.00; 75.90; 77.80; 79.50; 81.30; 83.10; 82.60; 84.80; 87.00; 89.30; 91.80; 94.20; 96.80; 99.50; 102.30; 105.20; 108.20; 111.30; 114.20; 117.10; 120.00]; 

l = zeros(1,length(MVAR));
for i = 1:length(MVAR) 
    l(i) = (Up40x^2)/(2*fn*pi*MVAR(i));
end

Inductance = transpose(l);

%VSR Tap position
Tap_Position = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30];

%VSR Losses
Rlosses = (0.2/100)*(Up40x^2)/MVAR(1);

% ----------------------------------------------------------------- %

Sn = 2*Sgob + 2*Scob;
agen = Sgob/Sn;
acon = Scob/Sn;

TQ = 5;
Ts = 50;
TLim = 10;

Vkra0 = 1;

a11 = (T2Wxbase + agen*x1TLbase)/Vkra0;
a12 = (x1TLbase*agen)/Vkra0;
a21 = a12;
a22 = a11;

a33 = (T401x23base*T401x13base)/((T401x13base - acon*T401x23base)*Vkra0);
a44 = (T402x23base*T402x13base)/((T402x13base - acon*T402x23base)*Vkra0);
 
bk1 = -acon*T401x23base/(T401x13base - acon*T401x23base);
bk2 = -acon*T402x23base/(T401x13base - acon*T401x23base);

bd1 = T401x13base/(T401x13base - acon*T401x23base);
bd2 = T402x13base/(T401x13base - acon*T401x23base);

A = [a11 a12 0 0; a21 a22 0 0; 0 0 a33 0; 0 0 0 a44];
B1 = [1; 1; bk1; bk2];
B2 = [0; 0; bd1; bd2];


C1 = A\B1;
C2 = A\B2;

KPG = Vkra0/(T2Wxbase + x1TLbase);
KIG = (Vkra0/Ts)*((1/XTSbase) + (1/(T2Wxbase + x1TLbase)));

KPC1 = Vkra0*(T401x23base*T401x13base)/(T401x13base - T401x23base);
KIC1 = (Vkra0/Ts)*((1/XTSbase) + (T401x23base*T401x13base)/(T401x13base - T401x23base));

KPC2 = Vkra0*(T402x23base*T402x13base)/(T402x13base - T402x23base);
KIC2 = (Vkra0/Ts)*((1/XTSbase) + (T402x23base*T402x13base)/(T402x13base - T402x23base));




