% ------------------------------------------------------------------------
% Matlab script for the technical note Inductor selection for DAB TN###
% Copyright : Imperix Switzerland ltd 2024
% Engineers : B. Sauter, J. Orsinger
% ------------------------------------------------------------------------

function dab_inductor_sizing_phaseShift()
% ------------------------
% User parameters
% ------------------------
% Full rated power (positive power implies power flow from V1->V2)
P = 1e3;                    %[W] rated active power

% DAB converter specs
V1 = 720;                   %[V] input (primary) dc voltage
V2 = 40;                    %[V] output (secondary) dc voltage
fsw = 100e3;                %[Hz] switching frequency (N.B. actual switching frequency of a single switch is fsw/2)
Ts = 1 / fsw;               %[s] switching period
f_clk = 250e6;              %[Hz] clock frequency of the modulator in FPGA 

% Transformer specs
Llk = 10e-6;                % [H] leakage inductance, referred to primary
n = 15;                     % [-] turns ration (N_prim/N_sec)

% Semiconductor specs
global Coss_pri Ron_pri Coss_sec Ron_sec
Coss_pri = 92e-12/10000;          % [F] equivalent output capacitance of MOSFET in primary full-bridge
Ron_pri = 80e-3;            % [Ohm] equivalent drain-source on-state resistance of MOSFET in primary full-bridge
Coss_sec = 260e-12;         % [F] equivalent output capacitance of MOSFET in secondary full-bridge
Ron_sec = 21e-3;            % [Ohm] equivalent drain-source on-state resistance of MOSFET in secondary full-bridge

% Minimum power point where converter must be operating in ZVS
P_min = 0.5 * P;

% Maximum allowable change in power when minimum phase shift adjustment is made [W] - max allowed is a percentage of full power
resolution_percentage = .025;
max_allowed_deltaP = resolution_percentage * P;

% ------------------------
% Lower bound on L for ZVS
% ------------------------
% Minimum L to achieve ZVS (use P_min because it will be a stricter minimum L)
% Selected inductor should be minimized while still meeting ZVS and resolution requirements at light load
% Coss primary is used because that is where ZVS conditions need achieved
Iout_min = P_min / V2;
Lmin = Coss_pri * (2 * Ts * V1 / ((1/n)*Iout_min*Ts + 8*Coss_pri*V1))^2;
disp(['Lmin (ZVS):', num2str(Lmin)])
a = -2;
b = 1;
c = -P_min * fsw * Lmin / n / V1 / V2;
% Solve for phase shift and controller resolution
phase_shift_min = 1 / f_clk / Ts;
phi = min(roots([a,b,c]));
max_deltaP = (n * V1 * V2 * phi * (1 - 2*phi) / fsw / Lmin) - (n * V1 * V2 * (phi - phase_shift_min) * (1 - 2*(phi - phase_shift_min)) / fsw / Lmin);

% whether or not the resolution limits the minimum L value
resolution_limited = false;

% if min L to achieve ZVS doesn't provide enough resolution iterate up until it does
while max_deltaP > max_allowed_deltaP
    resolution_limited = true;
    Lmin = Lmin * 1.0001;
    a = -2;
    b = 1;
    c = -P_min * fsw * Lmin / n / V1 / V2;
    phi = min(real(roots([a,b,c])));
    max_deltaP = (n * V1 * V2 * phi * (1 - 2*phi) / fsw / Lmin) - (n * V1 * V2 * (phi - phase_shift_min) * (1 - 2*(phi - phase_shift_min)) / fsw / Lmin);
end

if resolution_limited
    limiter = 'resolution requirement';
else
    limiter = 'ZVS condition at minimum power';
end

% ------------------------------------
% Upper bound on L for power transfer
% ------------------------------------

% Maximum L for max power transfer
phi_max = 0.25;
Lmax = (n * V1 * V2 * 0.25 * (1 - 2*phi_max)) / (P * fsw);
min_deltaP = (n * V1 * V2 * phi_max * (1 - 2*phi_max) / fsw / Lmax) - (n * V1 * V2 * (phi_max - phase_shift_min) * (1 - 2*(phi_max - phase_shift_min)) / fsw / Lmax);

% Feasibility check on lower/upper bounds
if Lmax < Lmin
    warning('on', 'all')
    warning(['Maximum L required to operate at full power is too small, limited by ', limiter, ' - loss model assumes ZVS and will, therefore, be inaccurate']);
    
    disp(['Maximum L = ', num2str(Lmax/1e-6), ' uH for phase shift = ', num2str(phi_max), ' @ full power: ', num2str(P), 'W but will not achieve ', limiter])
    disp(['Resolution @ full power (minimum delta P): ', num2str(min_deltaP), char(10)]);
else
    disp(['Minimum L = ', num2str(Lmin/1e-6), ' uH for phase shift = ', num2str(phi), ' @ min power: ', num2str(P_min), 'W', ' limited by ', limiter])
    disp('(Based on ZVS and resolution @ min rated power)')
    disp(['Resolution for minimum L (minimum delta P): ', num2str(max_deltaP) 'W', char(10)])
    
    disp(['Maximum L = ', num2str(Lmax/1e-6), ' uH for phase shift = ', num2str(phi_max), ' @ full power: ', num2str(P), 'W', ' limited by power transfer requirement'])
    disp(['Resolution for maximum L (minimum delta P): ', num2str(min_deltaP) 'W' , char(10)])
    disp(['Operating outside of specified conditions will yield a different operation than is modeled here', char(10)])
    
    % --------------------------------------------------------
    % Plot power losses vs resolution for inductance selection
    % --------------------------------------------------------
    Ploss_vs_resolution(P, Lmin, Lmax, fsw, n, V1, V2, Ts, phase_shift_min);
    Ploss_vs_resolution(P_min, Lmin, Lmax, fsw, n, V1, V2, Ts, phase_shift_min);
    
end

% ----------------------------
% Manual inductance selection
% ----------------------------
Ltot = input('Based on the provided information and your available inductance values, please select a total inductance value: ');
if Ltot < Lmin || Ltot > Lmax
    disp('WARNING: select a value between the reccomended bounds to achieve desired operating conditions')
end

% -----------------------------------------------------------
% Computation and display of resulting design characteristics
% -----------------------------------------------------------
% Calculate operating point at full power
a = -2;
b = 1;
c = -P * fsw * Ltot / n / V1 / V2;
% Solve for phase shift and controller resolution
phi = min(roots([a,b,c]));
resolution = (n * V1 * V2 * phi * (1 - 2*phi) / fsw / Ltot) - (n * V1 * V2 * (phi - phase_shift_min) * (1 - 2*(phi - phase_shift_min)) / fsw / Ltot);

% phase shift in seconds
t_phi = phi * Ts;

% Print physical inductor sizing
Lext = Ltot - Llk;
disp(['Selected inductance is ', num2str(10^6 * Ltot), 'uH, due to leakage Llk = ', num2str(Llk), ' required external inductance is ', num2str(10^6 * Lext), 'uH', char(10)])
disp(['At full power, required phase shift is ', num2str(phi), ' = ' , num2str(t_phi), ' seconds and resolution is ', num2str(resolution)])

% find max power for this inductor
P_max_inductanceChoice = n * V2 * V1 * 0.25 * (1-2*0.25) / fsw / Ltot;
disp(['For inductance choice, the max power that could be operated at is ', num2str(P_max_inductanceChoice), 'W'])

% find full ZVS range
P_min_inductanceChoice = (n * V2 / Ts) * (2 * Ts * V1 * sqrt(Coss_pri) / sqrt(Ltot) - 8 * Coss_pri * V1);
disp(['For inductance choice, the minimum power that could achieve ZVS (resolution requirement not guaranteed) is  ', num2str(P_min_inductanceChoice), 'W'])

% RMS current on the primary and secondary sides
I_pri_RMS = Irms_prim_phaseShift(phi, V1, V2, n, fsw, Ltot);
I_sec_RMS = I_pri_RMS * n;

% RMS current running through two switches at all times (assuming negligible dead time)
pri_Pcond = 2 * Ron_pri * I_pri_RMS^2;
sec_Pcond = 2 * Ron_sec * I_sec_RMS^2;
total_Pcond = pri_Pcond + sec_Pcond;

disp(['[At full power, using input inductance value] Expected primary RMS switch currents are: ', num2str(I_pri_RMS), ' and expected secondary RMS switch currents are: ', num2str(I_sec_RMS)])
disp(['Conduction losses are predicted to be:', char(10), num2str(pri_Pcond), 'W (primary)', char(10), num2str(sec_Pcond), 'W (secondary)', char(10), num2str(total_Pcond), 'W (total)', char(10)])

% Use all components from loss model to predict efficiency
Pin = P;
Pout = P - total_Pcond;
nu = (Pout / Pin);

disp(['With current elements of loss model, the predicted efficiency is ', num2str(nu*100), '% @ full power'])

% ---------------------------------------------------------
% Computation and plot of characteristics over power range
% ---------------------------------------------------------
Ploss_vs_P(P_min, P, Ltot, fsw, n, V1, V2)
end


function Ploss_vs_resolution(power, Lmin, Lmax, fsw, n, V1, V2, Ts, phase_shift_min)
% ------------------------------------------------------------
% Computation and plot of losses and resolution vs inductance
% ------------------------------------------------------------
global Ron_pri Ron_sec

% initialize vectors
phis = [];
Pconds = [];
resolutions = [];
Ls = [];

% iterate through inductor values
for L = Lmin : 1e-6 : Lmax
    Ls(end+1) = L;
    a = -2;
    b = 1;
    c = -power * fsw * L / n / V1 / V2;
    % Solve for phase shift and controller resolution
    phi = min(roots([a,b,c]));
    phis(end+1) = phi;
    resolutions(end+1) = (n * V1 * V2 * phi * (1 - 2*phi) / fsw / L) - (n * V1 * V2 * (phi - phase_shift_min) * (1 - 2*(phi - phase_shift_min)) / fsw / L);
    
    % RMS current on the primary and secondary sides
    I_pri_RMS = Irms_prim_phaseShift(phi, V1, V2, n, fsw, L);
    I_sec_RMS = I_pri_RMS * n;
    
    % RMS current running through two switches at all times (assuming negligible dead time)
    pri_Pcond = 2 * Ron_pri * I_pri_RMS^2;
    sec_Pcond = 2 * Ron_sec * I_sec_RMS^2;
    Pconds(end+1) = pri_Pcond + sec_Pcond;
end

% plot losses and resolution against inductance
figure
hold off
yyaxis left
plot(Ls*1e6, Pconds, 'color', [0 0.4470 0.7410])
hold on
ylabel('Losses [W]')
xlabel(['Total inductance [' char(181) 'H]'])
yyaxis right
plot(Ls*10^6, resolutions, 'color', [0.8500 0.3250 0.0980])
ylabel('Resolution [W]')
legend('Conduction losses', 'Resolution')
title(['Conduction loss/Resolution Tradeoff vs. Inductance @ ' num2str(power) 'W' ]);
grid on
end

function Ploss_vs_P(Pmin, Pmax, Ltot, fsw, n, V1, V2)
% -------------------------------------------------------------
% Computation and plot of losses and efficiency vs power level
% -------------------------------------------------------------
global Ron_pri Ron_sec

powers = [];
Pconds = [];
Iout_rms = [];
effs = [];

for P = Pmin : 1 : Pmax
    powers(end+1) = P;
    
    a = -2;
    b = 1;
    c = -P * fsw * Ltot / n / V1 / V2;
    % Solve for phase shift and controller resolution
    phi = min(roots([a,b,c]));
    
    % RMS current on the primary and secondary sides
    I_pri_RMS = Irms_prim_phaseShift(phi, V1, V2, n, fsw, Ltot);
    I_sec_RMS = I_pri_RMS * n;
    Iout_rms(end+1) = I_sec_RMS;
    
    % RMS current running through two switches at all times (assuming negligible dead time)
    pri_Pcond = 2 * Ron_pri * I_pri_RMS^2;
    sec_Pcond = 2 * Ron_sec * I_sec_RMS^2;
    Pconds(end+1) = pri_Pcond + sec_Pcond;
    
    % Efficiency
    effs(end+1) = 1 - Pconds(end) / P;
    
end

% plot conduction losses and resolution against power level
figure
hold off
yyaxis left;
plot(powers, Pconds, 'color', [0 0.4470 0.7410])
hold on
xlabel('Power [W]')
ylabel('Power losses [W]')
yyaxis right;
plot(powers, effs*100, 'color', [0.8500 0.3250 0.0980])
ylabel('Efficiency [%]')

legend('Conduction losses', 'Efficiency')
title(['Conduction losses and efficiency vs. Power Level @ L_{tot} = ' num2str(Ltot*1e6) char(181) 'H']);

grid on

end

function [Irms] = Irms_prim_phaseShift(phi,V1,V2,n,fsw,L)
% ------------------------------------------------------------
% Computation of primary RMS current for given operating point
% ------------------------------------------------------------

% source: Krismer_2011_03_17_Modeling_and_Optimization_of_Bidirectional_Dual_Active_Bridge_DC-DC_Converter_Topologies.pdf
    % (warning: phi is defined without unit here, whereas Krismer uses [rad] --> phi = phi_Krismer/2/pi)

% RMS current for phase-shift modulation (duty ratios are 0.5)
D1 = 0.5;
D2 = 0.5;

eRMS = (1-4*abs(phi)).*(1-2*abs(phi)+4*phi.^2-3*(D1.*(1-D1)+D2.*(1-D2)));
Irms = 1/2./fsw./L.*sqrt(D1.^2.*V1.^2.*(1-4/3.*D1)+D2.^2.*n^2.*V2.^2.*(1-4/3.*D2)+n*V1.*V2./3.*eRMS);

end