%% EEET 3016 Power Electronics (2023)
%% Lecture 7: Thyristors
% Dr. Yee Wei Law \langleyeewei.law@unisa.edu.au\rangle
% Table of Contents
% EEET 3016 Power Electronics (2023)
% Lecture 7: Thyristors
%     Transfer function
%     Example 3: RC snubber circuit for limiting dv/dt
%
%% Transfer function V_T(s)/V_S(s)
%

clearvars
syms V_s R_s L_s C_s R L s t
assume([V_s R_s L_s C_s R L], {'positive'});

Z = R_s + 1/(C_s*s);
G = collect(Z/(Z+R+(L_s+L)*s), s) % Transfer function V_T(s)/V_s(s)

%
%% Example 3: RC snubber circuit for limiting dv/dt
%
% Using G(s) = V_T(s)/V_S(s) from before:

syms V_T(s) v_T(t)
V_T = collect(G*(V_s/s), s)
v_T = simplify(ilaplace(V_T, s, t), 10)

% To prevent division by zero:

assumeAlso(C_s*R^2 + 2*C_s*R*R_s + C_s*R_s^2 - 4*L - 4*L_s > 0)
isAlways(C_s*R^2 + 2*C_s*R*R_s + C_s*R_s^2 - 4*L - 4*L_s > 0) % Double check

% Method 1: \dot{v}_T(t) = \mathcal{L}^{-1}(sV_T(s))

syms vdot_T(t)
vdot_T = collect(simplify(ilaplace(V_T*s, s, t)), t)

% Method 2: \dot{v}_T(t) using diff

syms vdot_T_(t)
vdot_T_ = collect(simplify(diff(v_T)), exp(-t/2*(R+R_s)/(L+L_s)))

% Confirm two methods are equivalent:

simplify(vdot_T - vdot_T_)

% Solving \ddot{v}_T(t)=0 yields no solution:

params = [V_s R L    L_s];
values = [200 5 1e-6 1e-6];
vddot_T = subs(diff(vdot_T), params, values);
solve(vddot_T==0, t, 'Real', true)

% Function solve failed to find an explicit solution for a stationary point. Assuming \dot{v}_T changes monotonically, then steady state requires 1️⃣ \lim_{t\to\infty}\dot{v}_T(t) = 0 and 2️⃣ the maximum \dot{v}_T to happen at t=0.
% ⚠ Caution: There are cases where \lim_{t\to\infty}\dot{v}_T(t) \neq 0, as shown:

expand(limit(vdot_T, t, Inf))

% However, let us assume steady state where \lim_{t\to\infty}\dot{v}_T(t) = 0, and thus in the absense of a stationary point, \dot{v}_T(0) is the maximum:

max_vdot_T = simplify(limit(vdot_T, t, 0))

% Since \dot{v}_T(0) depends only on R_s and not C_s, we can now determine R_s from the critical dv/dt:

R_ssol = solve(subs(max_vdot_T,params,values)==100/1e-6, R_s)

% With R_s fixed, we can now determine C_s from the desired \zeta:

zeta = subs((R+R_s)/2*sqrt(C_s/(L+L_s)), [params R_s], [values R_ssol])
C_ssol = double(solve(zeta==1/sqrt(2), C_s))

% Plots of v_T(t), \dot{v}_T(t), \ddot{v}_T(t):

trange = 0:1e-8:10e-6;
v_Tsol = subs(v_T, [params R_s C_s], [values R_ssol C_ssol]);
vdot_Tsol = subs(vdot_T, [params R_s C_s], [values R_ssol C_ssol]);
vddot_Tsol = subs(diff(vdot_T), [params R_s C_s], [values R_ssol C_ssol]);

figure
subplot(2,2,1); plot(trange, double(subs(v_Tsol, t, trange)))
ylabel('v_T')
subplot(2,2,2); plot(trange, double(subs(vdot_Tsol, t, trange)))
ylabel('dv_T/dt')
subplot(2,2,3); plot(trange, double(subs(vddot_Tsol, t, trange)))
ylabel('d^2v_T/dt^2')

% Clearly from the plots above, and output below, the maximum \dot{v}_T(t) exceeds the critical value of 100 V/\mus:

max(double(subs(vdot_Tsol, t, trange)))

%
% C_s and not R_s plays the dominant role in suppressing \dot{v}_T(t), so a sensible approach is to increase C_s until the maximum \dot{v}_T(t) falls below the critical value of 100 V/\mus:

for k=1:0.1:4
    vdot_Tsol = subs(vdot_T, [params R_s C_s], [values R_ssol k*C_ssol]);
    if max(double(subs(vdot_Tsol, t, trange))) <= 100/1e-6
        newzeta = double(subs(zeta, C_s, k*C_ssol));
        fprintf('Revised C_s=%.1fC_s, zeta=%.2fzeta=%.2f\n', ...
            k, newzeta/(1/sqrt(2)), newzeta)
        break
    end
end
trange = 0:1e-8:10e-6;
v_Tsol = subs(v_T, [params R_s C_s], [values R_ssol k*C_ssol]);
vdot_Tsol = subs(vdot_T, [params R_s C_s], [values R_ssol k*C_ssol]);
vddot_Tsol = subs(diff(vdot_T), [params R_s C_s], [values R_ssol k*C_ssol]);

figure
subplot(2,2,1); plot(trange, double(subs(v_Tsol, t, trange)))
ylabel('v_T')
subplot(2,2,2); plot(trange, double(subs(vdot_Tsol, t, trange)))
ylabel('dv_T/dt')
subplot(2,2,3); plot(trange, double(subs(vddot_Tsol, t, trange)))
ylabel('d^2v_T/dt^2')

% With the revised C_s, \zeta becomes 1.73 times higher, and the system becomes overdamped. We are thus opting to meet the safety criterion associated with the critical dv/dt at the expense of overdamping the transient response.
%
%

H = collect(Z/(Z+R), s) % G but with L=L_s=0
V_T = collect((V_s/s)*H, s)

% Method 1: \dot{v}_T(t) = \mathcal{L}^{-1}(sV_T(s))

vdot_T = collect(simplify(ilaplace(V_s*H, s, t)), V_s)

% Method 2: \dot{v}_T(t) using diff

v_T = collect(ilaplace(V_T, s, t), V_s)
vdot_T_ = collect(simplify(diff(v_T)), V_s)

% The results of both methods differ by one impulse term at t=0. Same as before, function solve failed to find an explicit solution for a stationary point:

params = [V_s R];
values = [200 5];
solve(subs(diff(vdot_T), params, values) == 100/1e-6, t, "Real", true)

% \dot{v}_T(0) depends on both R_s and C_s:

max_vdot_T = simplify(limit(vdot_T, t, 0))

% Luckily, R_s can be determined from the constraint on the discharge current:

R_ssol = double(subs(V_s/100, params, values))

% Now, C_s can be determined from \dot{v}_T(0) \leq 100 V/\mus:

C_ssol = double(solve(subs(max_vdot_T, [params R_s], [values R_ssol]) == 100/1e-6, C_s))

% Plots of v_T(t), \dot{v}_T(t), \ddot{v}_T(t) confirm the critical dv/dt is not exceeded:

trange = 0:1e-8:10e-6;
v_Tsol = subs(v_T, [params R_s C_s], [values R_ssol C_ssol]);
vdot_Tsol = subs(vdot_T, [params R_s C_s], [values R_ssol C_ssol]);
vddot_Tsol = subs(diff(vdot_T), [params R_s C_s], [values R_ssol C_ssol]);

figure
subplot(2,2,1); plot(trange, double(subs(v_Tsol, t, trange)))
ylabel('v_T')
subplot(2,2,2); plot(trange, double(subs(vdot_Tsol, t, trange)))
ylabel('dv_T/dt')
subplot(2,2,3); plot(trange, double(subs(vddot_Tsol, t, trange)))
ylabel('d^2v_T/dt^2')

%