V2G Smart Grid Integration MATLAB Simulation

%% SmartGrid_V2G_Parameters.m
% Run before opening the multi-EV Simulink model
% Projectsatbangalore | IEEE 2026 V2G Smart-Grid Integration

clear; clc;

%% ── FLEET DEFINITION ────────────────────────────────────
N_EV          = 8;                    % Number of EVs in the aggregator
P_rated_EV    = [7.4 7.4 11 11 22 7.4 11 22]*1e3; % W
Cbat_Ah       = [40 50 60 60 80 45 55 75];   % Ah
Vbat_nom      = 350*ones(1,N_EV);      % V
SOC_init      = [0.35 0.55 0.70 0.45 0.80 0.25 0.60 0.90];
SOC_min       = 0.15;
SOC_max       = 0.95;
departure_hr  = [2 4 6 3 8 1.5 5 7]; % hours from t=0

%% ── GRID / FEEDER ───────────────────────────────────────
Vgrid_ll      = 400;                  % V (line-line RMS)
f0            = 50;                   % Hz
S_base        = 500e3;                % VA (feeder rating)
R_line        = 0.02;                 % pu
X_line        = 0.08;                 % pu

%% ── RENEWABLE PLANTS ────────────────────────────────────
P_pv_rated    = 150e3;                % W
P_wind_rated  = 80e3;                 % W
% Simple irradiance / wind profiles (1-s resolution, 2 h)
t_prof        = 0:1:7200;
irr           = 800 + 200*sin(2*pi*t_prof/3600);  % W/m²
wind_spd      = 8 + 3*sin(2*pi*t_prof/1800);   % m/s

%% ── AGGREGATOR / CONTROL ────────────────────────────────
R_droop       = 0.04;                 % 4 % droop
deadband_f    = 0.1;                  % Hz
Ts_agg        = 5;                    % Aggregator sample time (s)
Ts_power      = 5e-6;                 % Power-stage sample time
comm_delay_ms = 120;                  % Typical ISO 15118 latency (ms)

%% ── LCL FILTER (per charger) ────────────────────────────
L1_filt       = 1.5e-3;
L2_filt       = 0.4e-3;
Cf_filt       = 15e-6;
Rd_filt       = 2.0;

disp('✔ Smart-Grid V2G parameters loaded. Open SmartGrid_V2G_Model.slx');

%% Aggregator_Dispatch.m — Embedded MATLAB Function for Simulink
% Inputs : SOC (1xN), f_grid, P_renewable, P_load, departure_flag (1xN)
% Outputs: Iref (1xN), mode (1xN)   mode: 0=standby, 1=G2V, 2=V2G

function [Iref, mode] = Aggregator_Dispatch(SOC, f_grid, P_ren, P_load, dep_flag)
    N      = length(SOC);
    Iref   = zeros(1,N);
    mode   = zeros(1,N);
    f0     = 50;
    R      = 0.04;
    db     = 0.1;
    P_rated = [7.4 7.4 11 11 22 7.4 11 22]*1e3;
    V_nom  = 230;

    % ── 1. Compute required fleet power from frequency ──
    delta_f = f_grid - f0;
    if abs(delta_f) < db
        P_freq = 0;
    else
        P_freq = -(1/R) * (delta_f/f0) * sum(P_rated);
        P_freq = max(min(P_freq, sum(P_rated)), -sum(P_rated));
    end

    % ── 2. Net power that fleet must supply ─────────────
    P_net = P_freq + (P_load - P_ren);   % +ve = need discharge

    % ── 3. Priority: protect low-SOC & respect departure ─
    avail_dis = (SOC > 0.25) & ~dep_flag;   % can discharge
    avail_chg = (SOC < 0.90);               % can charge

    if P_net > 0          % need V2G (discharge)
        weights = SOC .* avail_dis;         % prefer high SOC
        if sum(weights) > 0
            share = weights / sum(weights);
            P_alloc = share * P_net;
        else
            P_alloc = zeros(1,N);
        end
        for i = 1:N
            if avail_dis(i)
                mode(i) = 2;
                Iref(i) = -min(P_alloc(i)/V_nom, P_rated(i)/V_nom);
            end
        end
    else                 % need G2V (charge) or absorb surplus
        weights = (1-SOC) .* avail_chg;     % prefer low SOC
        if sum(weights) > 0
            share = weights / sum(weights);
            P_alloc = share * (-P_net);
        else
            P_alloc = zeros(1,N);
        end
        for i = 1:N
            if avail_chg(i)
                mode(i) = 1;
                Iref(i) = min(P_alloc(i)/V_nom, P_rated(i)/V_nom);
            end
        end
    end

    % ── 4. Force G2V for imminent departure ─────────────
    for i = 1:N
        if dep_flag(i) && SOC(i) < 0.60
            mode(i) = 1;
            Iref(i) = 0.5 * P_rated(i)/V_nom;  % half-rate charge
        end
    end
end

%% Latency_Sensitivity.m
% Sweeps ISO 15118 communication delay and records frequency nadir
% after a 0.3 Hz step disturbance. Quantifies control degradation.

SmartGrid_V2G_Parameters;   % load base params

delay_ms   = [0 50 100 150 200 300 500];
nadir_Hz   = zeros(size(delay_ms));
settle_s   = zeros(size(delay_ms));

for k = 1:length(delay_ms)
    % Set the Transport Delay block parameter via set_param
    set_param('SmartGrid_V2G_Model/CommDelay', 'DelayTime', ...
              num2str(delay_ms(k)/1000));
    
    % Run simulation (assumes model already configured for step Δf)
    simOut = sim('SmartGrid_V2G_Model', 'StopTime', '10');
    
    f_sig  = simOut.f_grid;           % timeseries from To Workspace
    nadir_Hz(k) = min(f_sig.Data);
    % settling time: first instant after which |f-50| < 0.05 Hz for 1 s
    err = abs(f_sig.Data - 50);
    idx = find(err < 0.05, 1, 'first');
    if ~isempty(idx)
        settle_s(k) = f_sig.Time(idx);
    else
        settle_s(k) = NaN;
    end
end

%% Plot
figure('Name','Latency Sensitivity');
subplot(2,1,1);
plot(delay_ms, nadir_Hz, 'o-', 'LineWidth', 2, 'Color', [0 .6 .6]);
ylabel('Frequency Nadir (Hz)'); grid on;
title('Effect of ISO 15118 Latency on Frequency Nadir');
subplot(2,1,2);
plot(delay_ms, settle_s, 's-', 'LineWidth', 2, 'Color', [.8 .3 0]);
xlabel('Communication Delay (ms)'); ylabel('Settling Time (s)');
grid on;
fprintf('Latency sweep complete. Results stored in nadir_Hz, settle_s\n');

%% Renewable_SelfConsumption.m
% Post-processing: calculates renewable self-consumption ratio
% and curtailment energy for a completed multi-EV simulation

% Assume To Workspace signals: P_pv, P_wind, P_fleet, P_grid
t      = P_pv.time;
P_ren  = P_pv.signals.values + P_wind.signals.values;
P_ev   = P_fleet.signals.values;   % positive = charging (absorbing)
P_exp  = max(P_grid.signals.values, 0);  % export to utility

% Energy calculations (trapezoidal integration)
E_ren  = trapz(t, P_ren) / 3.6e6;   % kWh
E_abs  = trapz(t, max(P_ev,0)) / 3.6e6;
E_curt = trapz(t, P_exp) / 3.6e6;

SCR = 100 * E_abs / max(E_ren, eps);   % Self-Consumption Ratio %

fprintf('\n── Renewable Integration Metrics ──────────────────\n');
fprintf('Total renewable energy : %8.2f kWh\n', E_ren);
fprintf('Absorbed by EV fleet  : %8.2f kWh\n', E_abs);
fprintf('Exported / curtailed  : %8.2f kWh\n', E_curt);
fprintf('Self-Consumption Ratio: %8.1f %%\n', SCR);
fprintf('───────────────────────────────────────────────────\n');

figure('Name','Renewable vs Fleet Power');
plot(t/3600, P_ren/1e3, 'g', 'LineWidth', 1.5); hold on;
plot(t/3600, P_ev/1e3,  'b', 'LineWidth', 1.5);
plot(t/3600, P_exp/1e3, 'r--','LineWidth', 1.2);
legend('P_{renewable}', 'P_{fleet} (charge)', 'P_{export}');
xlabel('Time (h)'); ylabel('Power (kW)');
title(['Renewable Self-Consumption = ' num2str(SCR,'%.1f') '%']);
grid on;

%% Islanded_Frequency_Restore.m
% Demonstrates secondary frequency restoration after islanding
% Uses a simple discrete PI on the aggregator side

% Parameters
f0      = 50;
Kp_sec  = 200;      % W/Hz  (secondary proportional)
Ki_sec  = 50;       % W/(Hz·s)
Ts      = 0.1;      % secondary sample time
Tsim    = 60;       % s

% Simulate a step load increase at t=5 s while islanded
t       = 0:Ts:Tsim;
P_load  = 30e3 * ones(size(t));   % 30 kW base
P_load(t>=5) = 45e3;             % +15 kW step

P_ev    = zeros(size(t));
f       = f0 * ones(size(t));
int_err = 0;
H_eq    = 3;                     % equivalent inertia (s)
D       = 0.02;                  % damping

for k = 2:length(t)
    % Primary droop response (already included in P_ev via R)
    delta_f = f(k-1) - f0;
    P_prim  = -(1/0.04) * (delta_f/f0) * 80e3;  % fleet rating 80 kW
    
    % Secondary PI
    int_err = int_err + delta_f * Ts;
    P_sec   = -Kp_sec * delta_f - Ki_sec * int_err;
    
    P_ev(k) = P_prim + P_sec;
    P_ev(k) = max(min(P_ev(k), 80e3), -80e3);
    
    % Simple swing equation
    dP = P_ev(k) - P_load(k);
    df = (dP/(2*H_eq*1e5) - D*delta_f) * Ts;  % scaled
    f(k) = f(k-1) + df;
end

figure('Name','Islanded Frequency Restoration');
subplot(2,1,1);
plot(t, f, 'b', 'LineWidth', 1.8);
yline(50, 'k--'); ylabel('Frequency (Hz)');
title('Secondary Control Restores Frequency after Load Step'); grid on;
subplot(2,1,2);
plot(t, P_ev/1e3, 'Color', [0 .6 .5], 'LineWidth', 1.8);
xlabel('Time (s)'); ylabel('Fleet Power (kW)'); grid on;