V2G and G2V MATLAB Simulation

%% V2G_G2V_Parameters.m  — Run BEFORE opening the Simulink model
% Projectsatbangalore | IEEE 2026 V2G/G2V Project

%% ── BATTERY PARAMETERS ──────────────────────────────────
clear; clc;
Vbat_nom   = 350;        % Nominal battery voltage (V)
Cbat_Ah    = 60;         % Battery capacity (Ah)
SOC_init   = 0.50;       % Initial SOC (50%)
SOC_min    = 0.20;       % Minimum SOC before forced G2V
SOC_max    = 0.90;       % Maximum SOC before V2G enabled
R0         = 0.01;       % Internal resistance — ohm (ESR)
R1         = 0.005;      % R-C branch resistance
C1         = 1500;       % R-C branch capacitance (F)

%% ── DC-DC CONVERTER ─────────────────────────────────────
Vbus       = 700;         % DC bus voltage target (V)
L_dc       = 1.5e-3;     % Boost inductor (H)
Cin        = 470e-6;     % Input capacitor (F)
Cout       = 1000e-6;    % Output/bus capacitor (F)
fsw_dc     = 20e3;       % DC-DC switching frequency (Hz)
Kp_dcdc    = 0.8;        % PI proportional gain
Ki_dcdc    = 50;         % PI integral gain

%% ── H-BRIDGE VSI ────────────────────────────────────────
Vgrid_rms  = 230;         % Grid phase voltage RMS (V)
f_grid     = 50;          % Grid frequency (Hz)
fsw_inv    = 10e3;        % VSI switching frequency (Hz)
Ts         = 1/fsw_inv/20; % Simulation timestep (5 µs)

%% ── LC / LCL FILTER ─────────────────────────────────────
L1_filt    = 2e-3;        % Inverter-side inductance (H)
L2_filt    = 0.5e-3;     % Grid-side inductance (H)
Cf_filt    = 10e-6;      % Filter capacitance (F)
Rd_filt    = 2.5;         % Damping resistance (Ω)
f_res      = 1/(2*pi) * sqrt((L1_filt+L2_filt)/(L1_filt*L2_filt*Cf_filt));
fprintf('LCL Resonant frequency: %.1f Hz\n', f_res);

%% ── DROOP / FREQUENCY CONTROL ───────────────────────────
R_droop    = 0.04;        % Droop coefficient (4%)
P_rated    = 7400;        % EV charger rated power (W, 7.4 kW)
H_inertia  = 4;           % Virtual inertia constant (s)
deadband_f = 0.1;         % Frequency deadband (Hz)
Tv         = 0.02;        % Virtual inertia filter time const (s)

disp('✔ V2G/G2V parameters loaded into workspace. Open V2G_Model.slx');

%% SOC_StateLogic.m — Embedded MATLAB Function block in Simulink
% Inputs: SOC (0-1), f_grid (Hz)
% Outputs: mode (0=standby,1=G2V,2=V2G), Iref (A)

function [mode, Iref] = SOC_StateLogic(SOC, f_grid)
    % Persistent variable for mode memory
    persistent current_mode;
    if isempty(current_mode); current_mode = 0; end

    I_CC  = 50;   % CC phase current = 50A (~0.83C for 60Ah bat)
    I_max = 32;   % Max grid current RMS (A) — 7.4kW/230V
    f0    = 50.0; % Nominal frequency
    db    = 0.1;  % Deadband (Hz)
    R     = 0.04; % Droop (4%)

    % ── Priority 1: Battery protection ──────────────────
    if SOC < 0.20
        mode  = 1;      % Force G2V (critical charge)
        Iref  = I_CC;   % CC phase

    % ── Priority 2: Full battery, shed to grid ──────────
    elseif SOC > 0.95
        mode  = 2;      % V2G: small discharge to avoid overcharge
        Iref  = -5;    % Trickle discharge (A, negative = discharge)

    % ── Priority 3: Frequency regulation ────────────────
    elseif (f_grid < f0 - db) && (SOC > 0.30)
        % Frequency low → inject power (V2G)
        delta_f = f0 - f_grid;
        delta_P = (1/R) * (delta_f / f0) * 7400; % Watts
        mode  = 2;
        Iref  = -min(delta_P/230, I_max); % A (negative)

    elseif (f_grid > f0 + db) && (SOC < 0.85)
        % Frequency high → absorb power (G2V fast charge)
        delta_f = f_grid - f0;
        delta_P = (1/R) * (delta_f / f0) * 7400;
        mode  = 1;
        Iref  = min(delta_P/230, I_max); % A (positive)

    else
        % Normal G2V smart charging (off-peak tariff)
        mode  = 1;
        Iref  = I_CC * 0.5; % 50% rate for off-peak
    end
    current_mode = mode;
end

%% THD_FFT_Analysis.m — Run after Simulink simulation completes
% Analyses grid current harmonic content and computes THD

%% Load simulation output (from To Workspace block "Igrid_out")
t     = Igrid_out.time;
Igrid = Igrid_out.signals.values;

%% Steady-state window (skip first 3 cycles = 60ms)
f0     = 50;  fs = 1/Ts;
t_start = 0.06;  t_end = t_start + 10/f0; % 10 cycles
idx    = (t >= t_start) & (t <= t_end);
i_ss   = Igrid(idx);
N      = length(i_ss);

%% FFT
Y      = fft(i_ss, N);
f_axis = (0:N-1) * fs / N;  % Frequency axis
Ymag   = 2/N * abs(Y(1:N/2));
f_axis = f_axis(1:N/2);

%% Find fundamental and harmonics
[~, idx50] = min(abs(f_axis - f0));
I1_mag     = Ymag(idx50);       % Fundamental amplitude
harm_ord   = [3 5 7 9 11 13]; % Odd harmonics to check
THD_sq     = 0;

fprintf('\n── Harmonic Analysis ────────────────────────────────\n');
fprintf('Order  Frequency(Hz)  Magnitude(A)  %% of Fundamental\n');
for h = harm_ord
    [~, idxh] = min(abs(f_axis - h*f0));
    Ih = Ymag(idxh);
    pct = 100*Ih/I1_mag;
    THD_sq = THD_sq + (Ih/I1_mag)^2;
    fprintf('  %2d       %4d          %6.3f A       %5.2f%%\n', h, h*f0, Ih, pct);
end
THD = 100 * sqrt(THD_sq);
fprintf('────────────────────────────────────────────────────\n');
fprintf('Total THD = %.2f%%  (IEEE 519 limit: <5%%)\n', THD);

%% Plot
figure('Name','V2G/G2V Harmonic Spectrum');
subplot(2,1,1); plot(t(idx), i_ss, 'b', 'LineWidth', 1.2);
xlabel('Time (s)'); ylabel('Current (A)');
title('Grid Current Waveform (steady state)'); grid on;

subplot(2,1,2); stem(f_axis(1:300), Ymag(1:300), 'filled', 'Color', [.85 .4 0]);
xlabel('Frequency (Hz)'); ylabel('Magnitude (A)');
title(['Harmonic Spectrum — THD = ', num2str(THD,'%.2f'), '%']); grid on;
xlim([0 1500]);

%% V2G_Droop_Simulation.m  
% Sweeps grid frequency from 49 to 51 Hz, plots V2G/G2V power response
% Validates droop characteristic P vs Δf

V2G_G2V_Parameters;   % Load base parameters

f_range  = 49.0 : 0.05 : 51.0;   % Grid frequency sweep
P_droop  = zeros(1, length(f_range));
P_inert  = zeros(1, length(f_range));

for k = 1:length(f_range)
    f_now    = f_range(k);
    delta_f  = f_now - 50;
    
    % Droop response (primary)
    if abs(delta_f) < deadband_f
        dP = 0;
    else
        dP = -(1/R_droop) * (delta_f/50) * P_rated;
        dP = max(min(dP, P_rated), -P_rated);   % Clamp
    end
    P_droop(k) = dP;
    
    % Virtual inertia (simplified, df/dt estimated)
    if k > 1
        dfdt = (f_range(k) - f_range(k-1)) / 0.05;
        P_inert(k) = P_droop(k) - 2*H_inertia * dfdt;
    else
        P_inert(k) = P_droop(k);
    end
end

%% Plot droop characteristic
figure('Name','V2G Droop Characteristic');
plot(f_range, P_droop/1000, 'r-', 'LineWidth', 2); hold on;
plot(f_range, P_inert/1000, 'b--','LineWidth', 2);
yline( P_rated/1000,'k:','P\_max G2V');
yline(-P_rated/1000,'k:','P\_max V2G');
xline(49.9,'--'); xline(50.1,'--');
fill([49.9 50.1 50.1 49.9], [-8 -8 8 8], 'y','FaceAlpha',.15);
legend('Droop only','Droop + Virtual Inertia');
xlabel('Grid Frequency (Hz)'); ylabel('EV Power (kW)');
title('V2G/G2V Power vs Grid Frequency — Droop Characteristic (R=4%)');
grid on; grid minor;
% +kW = G2V (charging), -kW = V2G (discharging to grid)

%% PLL_Design.m — SRF-PLL design and Bode plot verification
% Synchronous Reference Frame PLL (SRF-PLL) for V2G grid sync

omega0 = 2*pi*50;   % Nominal grid angular frequency (rad/s)
% PLL PI controller gains (tuned for bandwidth ~30 Hz)
Kp_pll = 132;       % Proportional gain  (rad/s per rad)
Ki_pll = 3600;      % Integral gain      (rad/s² per rad)

% Transfer function of PLL open loop: Kp + Ki/s × 1/s (VCO integrator)
s     = tf('s');
C_pll = Kp_pll + Ki_pll/s;       % PI controller
VCO   = 1/s;                     % VCO = integrator
G_ol  = C_pll * VCO;             % Open-loop TF

% Closed loop
G_cl  = feedback(G_ol, 1);
fprintf('PLL Closed-Loop Poles:\n'); disp(pole(G_cl));

% Check bandwidth and phase margin
[Gm, Pm, Wcg, Wcp] = margin(G_ol);
fprintf('Phase Margin: %.1f deg at %.1f Hz\n', Pm, Wcp/(2*pi));
fprintf('Gain  Margin: %.1f dB\n', 20*log10(Gm));

% Step response: grid phase jump of 30°
theta_jump = pi/6;    % 30 degrees in radians
t_sim = 0:0.0001:0.1;
figure('Name','PLL Step Response');
step(theta_jump * G_cl, t_sim);
title('SRF-PLL Phase Tracking — 30° Phase Jump');
xlabel('Time (s)'); ylabel('Phase angle (rad)'); grid on;
% Expected: settle within 2 cycles (~40 ms) with <5% overshoot

% Bode plot
figure('Name','PLL Open-Loop Bode');
bode(G_ol); grid on;
title('PLL Open-Loop Bode — verify Pm > 45°, BW ~ 200 rad/s');